Материал: Спецглавы высшей математики. численные методы. Пантелеев И.Н

Подставляя (V.16) в выражение для многочлена (V.15), получаем

Формула (V.17) называется первым интерполяционным полиномом Ньютона.

Пример V.3. Решить, используя пакет MATLAB, задачу интерполяции с помощью первого интерполяционного полинома Ньютона для функции f(х) = sin (x), заданной таблично в восьми точках на интервале [0, 2π].

1. Создайте файл Newtonl.m (листинг V.5), содержащий описание функции, возвращающей локальное значение интерполяционного полинома Ньютона.

Листинг V.5. Файл Newton1. m

f unction z=Newtonl ( t , х, у)

%t - абсцисса точки, в которой вычисляется

%значение интерполяционногополинома

%х, у - координаты точек, заданных таблично

N=length(x) ; for i=l:N

f(i,1) = y(i)

end;

for k=2:N

for i=l:N-k+l

f(i,k) = (f(i+1,k-1)-f(i,k-1)) / (x(i+k-1) –x(i)); end;

end;

s=y(l);

for k=2:N

118

r=l;

for i=l:k-l

r=r.*(t-x(i)); end; s=f(l,k)*r+s;

end;

z=S;

2. Задайте табличные значения интерполируемой функции.

»N=8;

»i=1:N

»x(i)=2*pi / (N-1)*(i-1);

»y=sin(x)

3. Задайте значения абсцисс точек, в которых вычисляется значение интерполяционного полинома.

»М=1000;

»j=1:M

»x(j)=2*pi / (M-1)*(j-1);

»Y=sin(X); % вычисление точных значений интерполируемой функции

4.Вычислите значения полинома Ньютона в узлах заданной координатной сетки.

» for k=l:M Z(j)=Newtonl(X(j),x,y); end;

5.Постройте разности между точными и интерполированными значениями функции (рис. V.4).

» plot(X,Y-Z)

Отдавая дань традициям преподавания численных методов в прошлом веке, приведем описание модификации формулы (V.17), применявшейся при ручных вычислениях. Положим

119

x hx0 t , т.е. x x0 ht , тогда:

Подставляя данные выражения в (V.17), окончательно получаем

120

Рис. V.4. Разность между точными и интерполированными значениями функции с помощью первого полинома Ньютона

Формула (V.18) называется первой интерполяционной формулой Ньютона. Данная формула применяется для интегрирования в начале отрезка, когда t мало по абсолютной величине.

5.3.3. Втораяинтерполяционнаяформула Ньютона

Когда значение аргумента находится ближе к концу отрезка интерполяции, используется вторая интерполяционная формула Ньютона, которая получается, если искать интерполяционный полином в виде:

Коэффициенты полинома (V.19), находятся из условия совпадения значений функции и интерполяционного многочлена в узлах:

ной формулы Ньютона, используемой при ручных вычислениях:

Pn (x) Pn (xn th) yn t yn 1 t(t2!1) 2 yn 2

(V.21)

t(t 1) n!(t n 1) n y0

Пример V.4. Решить, используя пакет MATLAB, задачу интерполяции с помощью второго интерполяционного полинома Ньютона для функции f(x) = sin(x), заданной таблично в восьми точках на интервале [0, 2π].

1. Создайте файл Newton2.m (листинг V.6), содержащий описание функции, возвращающей локальное значение интерполяционного полинома Ньютона.

Листинг V.6. Файл Newton2. m function z=Newton2 ( t , х, у) N=length(x) ;

for i=l:N

f (i,1) = y(i); end;

for k=2:N

for i=l:N-k+l

f(i ,k ) =(f(i+1, k-1)-f(i, k-1)) / (x(i+k-1)-x(i)); end;

end;

s=y(l);

122

for k=l:N r=l;

for i=l:k-l

r = r*(t-x(N-i) ) ; end;

s=f (N-k,k)*r+s end;

z=s;

2.Задайте табличные значения интерполируемой функции. » N=8;

» i =1:N;

» x (i) =2 * pi / (N-1) * (i-1); » y= sin (x);

3.Задайте значения абсцисс точек, в которых вычисляется значение интерполяционного полинома.

» М=1000; » j=l:M;

» X(j)=2*pi/(M-l)*(j-l);

» Y=sin(X); % вычисление точных значений

4.Вычислите значения полинома Ньютона в узлах заданной координатной сетки.

» for k=l:M Z(j)=Newtonl(X(j) ,x,y) ; end;

5.Постройте разности между точными интерполированными значениями

функции (рис. V.5). » plot(X,Y-Z)

123

Рис. V.5. Разность между точными и интерполированными значениями функции с помощью второго полинома Ньютона

5.4. Погрешность интерполяции

Погрешность интерполяции полиномом Лагранжа оценивается по формуле:

Погрешность интерполяции полиномом Ньютона оценивается

по формулам:

124

5.5. Сплайн-интерполяция

При большом количестве узлов интерполяции приходится использовать интерполяционные полиномы высокой степени, что создает определенные неудобства при вычислениях. Можно избежать высокой степени интерполяционного многочлена, разбив отрезок интерполяции на несколько частей и построив на каждой части самостоятельный интерполяционный многочлен. Однако такое интерполирование обладает существенным недостатком: в точках сшивки разных интерполяционных полиномов будет разрывной их первая производная, поэтому для решения задачи кусочно-линейной интерполяции используют особый вид кусочно-полиномиальной интерполяции — сплайнинтерполяцию.

Сплайн — это функция, которая на каждом частичном отрезке интерполяции является алгебраическим многочленом, а на всем заданном отрезке непрерывна вместе с несколькими своими производными.

Пусть интерполируемая функция f(x) задана своими

значениями у, в узлах xi , (i 0,1, ,n). Длину частичного

отрезка xi 1, xi обозначим hi :

hi xi xi 1,(i 1,2, ,n).

где ai,bi,ci,di— четверка неизвестных коэффициентов. Можно доказать, что задача нахождения кубического сплайна имеет единственное решение. Потребуем совпадения значений S(x) в узлах с табличными значениями функции f(x):

Число этих уравнений 2π в два раза меньше числа неизвестных коэффициентов. Для того чтобы получить дополнительные условия, потребуем также непрерывности первой и второй производных сплайна во всех точках, включая узлы. Для этого следует приравнять левые и правые производные S' (x 0), S' (x 0), S'' (x 0), S' (x 0), S'' (x 0) во

внутреннем узле xi .

Вычислив выражения для производных S'(x), S''(x) последовательным дифференцированием (V.27):

Приравняв левые и правые производные, получаем:

получения недостающих уравнений накладывают требования к поведению сплайна на концах отрезка интерполяции. Если потребовать нулевой кривизны сплайна на концах отрезка интерполяции (т. е. равенство нулю второй производной), то получим:

Исключив из уравнений (V.28)-(V.33) я неизвестных at , получаем систему уравнений:

c1 0 ,

cn 3dnhn 0 ,

где i 0, 1, ,n 1 .

Система (V.35) состоит из Зn уравнений. Решив систему (V.35), получаем значения неизвестныхbi , ,ci ,di ,

определяющих совокупность всех формул для искомого интерполяционного сплайна

Si (x) yi 1 bi( (x xi 1) ci (x xi 1)2 di (x xi 1)3 , (V.36)

где i 0, 1, ,n 1 .

Программа, реализующая метод сплайн-интерполяции, оказывается достаточно громоздкой, поэтому мы ограничимся обсуждением решения задачи об интерполяции синуса с помощью кубических сплайнов, используя функцию пакета MATLAB spline().

Пример V.5. Решить, используя пакет MATLAB, задачу интерполяции с помощью кубических сплайнов для функции f(x) = sin(x), заданной таблично в восьми точках на интервале

[0, 2π].

1. Задайте табличные значения интерполируемой функции.

127