Материал: Спецглавы высшей математики. численные методы. Пантелеев И.Н
»N=8;
»i = 1:N;
»x (i) = 2*pi / (N-1)*(i-1);
»y = sin (x) ;
2.Задайте значения абсцисс точек, в которых вычисляется значение интерполяционного полинома.
» М=1000; » j = 1:M;
» x (j) = 2*pi / (N-1)*(i-1)
» Y=sin (X) ; % вычисление точных значений интерполируемой функции
3.Вычислите интерполируемые значения функции в узлах координатной сетки.
» yy=spline(x,y,X);
4.Визуализируйте результаты сплайн-интерполяции и разности между точными и интерполированными значениями
(рис. V.6, V.7).
» plot(x,y, 'o'.X ,yy); » plot(X.Y-yy);
5.Вычислите и визуализируйте значения первых производных сплайна (рис. V.8).
» m=1:М-1
» yyl(m) = (yy(m+l)-yy(m))/(2*pi/(M-l)); » plot(X(m),yyl(m))
6.Вычислите и визуализируйте значения вторых производных сплайна (рис. V.9).
» m=1:М-2;
» уу2(m)=(yyl(m+1)-yyl(m))/(2*pi/(M-1)); » plot (X (m) ,уу2(m))
Рис. V.6. Визуализация исходных данных и результатов сплайн-интерполяции
129
128
Рис. V.7. Разность между точными и интерполированными значениями функции f (х) = sin (x)
Рис. V.8. Первая производная, вычисленная по численным значениям, полученным сплайн-интерполяцией
130
Рис. V.9. Вторая производная, вычисленная по численным значениям, полученным сплайн-интерполяцией
Рис. V.10. Третья производная, вычисленная по численным значениям, полученным сплайн-интерполяцией
131
7. Вычислите и визуализируйте значения третьих производных |
» N=8; |
|
сплайнов (рис. V.10). |
» i=l:N; |
|
» ууЗ (т) = (уу2 (m+1)-уу2 (m)) / (2*pi/ (M-1)) ; » plot(X(т),ууЗ); |
» x(i)=2*pi/(N-l)*(i-1); |
|
» m=1:М-3; |
» y=sin(x); |
|
Как видно из рис. V.7—V.10, первая и вторая производные |
|
|
сплайнов являются непрерывными функциями, третья и более |
2. Задайте значения абсцисс точек, в которых вычисляется |
|
высокого порядка производные — разрывными функциями. |
значение интерполяционного полинома. |
|
5.6. Решение задачи одномерной интерполяции средствами |
» М=1000; |
|
» j=l:M, |
||
|
пакета MATLAB |
» X(j)=2*pi/(M-l)*(j-1); |
Для решения задачи одномерной сплайн-интерполяции в |
» Y=sin(X); % вычисление точных значений интерполируемой |
|
пакете MATLAB используется функция interp 1 (), имеющая |
функции |
|
следующий синтаксис: |
|
|
yi = interpi (х,у,xi) — линейная интерполяция табличных |
3. Вычислите интерполируемые значения функции в узлах |
|
значений х, у в точках, абсциссы которых находятся в векторе |
координатной сетки и визуализируйте точные исходные |
|
xi. |
|
данные и точные интерполированные значения (рис. V.11) |
yi = interpi (у,xi) — линейная интерполяция табличных |
» yi=interpl(x,y,X); |
|
значений у в точках, абсциссы которых находятся в векторе xi |
» plot(x,y,'о1,X,Y,X,yi) |
|
в предположении, что |
|
|
x=l:length(у). |
|
|
yi = interpi (х,у,xi,method) — интерполяция линейных значений |
|
|
с использованием выбранного метода интерполяции. |
|
|
Возможные значения переменной method: |
|
|
'nearest' — интерполяция с использованием ближайших |
|
|
узлов; |
|
|
'linear' — линейная интерполяция (по умолчанию); |
|
|
'spline' — интерполяция кубическим сплайном; |
|
|
'pchip' — интерполяция полиномами Эрмита третьей степени; |
|
|
' cubic' — аналогично pchip. |
|
|
Пример 5.6. Решить задачу интерполяции для функции f (х) = |
|
|
sin |
(х), заданной таблично в восьми точках на интервале |
|
[0, 2π], используя встроенную функцию пакета MATLAB |
|
|
interpi (). |
|
|
1. |
Задайте табличные значения интерполируемой функции. |
Рис. V.11. Исходные данные и результат линейной интерполяции |
|
|
|
132 |
133 |
|
VI. ЧИСЛЕННОЕ ДИФФЕРЕНЦИРОВАНИЕ И ИНТЕГРИРОВАНИЕ
В данном разделе рассмотрим формулы численного дифференцирования и интегрирования (формулу прямоугольников, формулу трапеций, формулу Симпсона), погрешности данных формул, изложим основные идеи вычисления определенных интегралов методом Монте-Карло, обсудим использование соответствующих функций пакета
MATLAB.
6.1. Численное дифференцирование функций, заданных аналитически
По определению производная функции f(x) равна:
f ' (x) lim |
f (x x) f (x) |
(VI.1) |
|
x |
|||
0 |
|
Переходя в (VI.1) от бесконечно малых разностей к конечным,
получаем |
приближенную |
формулу |
численного |
||
дифференцирования: |
f (x x) f ()x |
|
|
||
|
f ' (x) |
. |
(VI.2) |
||
|
|
||||
|
|
x |
|
|
|
Формула (VI.2) позволяет построить простой вычислительный алгоритм:
1.Задать значение точки, в которой вычисляется производная.
2.Задать значение приращения Ах.
3.Вычислить производную в соответствие с формулой (VI.2). Замена бесконечно малых приращений конечными является причиной возникновения ошибки. Для оценки ее величины
разложим функцию f(x) в точке x x в ряд Тэйлора:
f (x x) f (x) |
f ' (x) |
x |
f '' (x) |
( x)2 |
|
f (2) |
( x)n (VI.3) |
|
1! |
2! |
n! |
||||||
|
|
134 |
|
|
||||
|
|
|
|
|
|
|
Подставив (VI.3) в (VI.2) и приведя подобные члены, получим:
f ' (x) f ' (x) |
f '' (x) |
x |
(VI.4) |
|
2! |
||||
|
|
|
Из (VI.4) видно, что все члены начиная со второго, определяют отличие численного значения производной от ее точного значения. Основной член погрешности равен
f '' (x) x / 2!, т. к. данный член ~ x , говорят, что формула
(VI.2) имеет первый порядок точности по x .
Можно вычислять производную, используя симметричную разностную схему:
f ' (x) |
f (x x) f (x x) |
. |
(VI.5) |
|
|||
|
2 x |
|
|
Для оценки точности данной формулы необходимо удержать первые четыре члена в разложении в ряд Тэйлора:
|
|
|
1 |
|
|
|
x |
|
|
|
( x) |
2 |
|
|
( x) |
3 |
|
|
f |
' (x) |
f (x) f ' |
(x) |
f '' (x) |
|
f ''' (x) |
|
|
|
|||||||||
|
1! |
2! |
3! |
|
|
|||||||||||||
|
|
|
2 x |
|
|
|
|
|
|
|
|
|
|
|
||||
|
1 |
|
|
x |
|
|
|
( x) |
2 |
|
|
( x) |
3 |
|
|
|||
|
f (x) f ' (x) |
f '' (x) |
|
f ''' |
(x) |
|
. |
(VI.6) |
||||||||||
|
1! |
2! |
|
3! |
|
|||||||||||||
|
2 x |
|
|
|
|
|
|
|
|
|
|
|
||||||
Раскрыв в (VI.6) скобки и приведя подобные члены, получаем:
f ' (x) |
f ' (x) f ''' (x) |
( x)2 |
|
(VI.7) |
|
3! |
|||||
|
|
|
|
Из (VI.7) видно, что основной член погрешности равен f ''' ( x )( x )2 / 3! , т. к.
данный член ~ ( x)2 , говорят, что формула (VI.5) имеет
второй порядок точности по x .
Пример VI.1. Вычислить значения производной функции f(x) = sin(0.01x2) на интервале [0, 10π], используя пакет MATLAB.
135
Решение:
» f = inline ('sin(0.01*х.^2)'); % задание дифференцируемой функции
» dx=0.01; % шаг изменения координатной сетки » x=0:dx:10*pi; % вычисление координат узлов
»yf=feval (f ,x); % вычисление значений функции в узлах % выполнение процедуры численного дифференцирования
»N=length(x);
»m=l:N-l;
»df(m)=(yf(m+1)-yf(m))/dx;
»plot(df) % визуализация производной функции (рис. VI.1)
»fl=inline('0.02*x.*cos(0.01*x.^2)'); % задание функции,
%описывающей первую производную
»ya=feval(£l,x); % вычисление значений первой производной
%по аналитической формуле
»plot(x(m),abs(yf(m)-уа(m)); % визуализация разности между
%численными и аналитическими
%значениями производной (рис. VI.2)
Рис. VI.1. График производной функции f (x) = sin (0.01x2)
Аналогичным образомпоступаютпри вычислении производных высшихпорядков.
Отметим, что для аппроксимации производных конечными разностями в пакете MATLAB имеется функция diff(), синтаксис которой имеет следующийвид:
diff(x) — возвращает конечные разности, вычисленные по смежным элементам вектора х. Длина вектора, возвращаемого функцией diff(x), на единицу меньше длины вектора х. Если х
— матрица, то возвращает матрицу, содержащую конечные разности, вычисленные по каждому столбцу;
diff (x,n) — возвращает конечные разности л-го порядка;
diff (x,n,dim) — возвращает конечные разности n-ro порядка по столбцам (dim=i) или строкам матрицы.
Рис. VI.2. Разность между точными и численными значениями производных функции f (х) = sin (0.01x2)
137
136