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

V.ИНТЕРПОЛИРОВАНИЕФУНКЦИЙ

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

5.1. Постановказадачи

Пусть известные значения некоторой функции f(х) образуют следующую таблицу:

Исходные данные в задаче интерполяции

Требуется получить значение функции f(x) для значения аргумента xє [ x0, xn ], несовпадающего ни с одним из значений хi (i = 0, 1, ..., n).

Решение задачи сводится к поиску некоторой приближающей функции F(х), близкой в определенном смысле к функции f(x), для которой известно аналитическое выражение.

Классический подход к решению задачи построения приближающей функции основан на требовании строгого совпадения значений функций f(x) и F(x) в точках x,(i = 0, 1,...,n)

108

В данном случае нахождение приближенной функции называется интерполированием, а точки х0, x1 ,..., хn называются узлами интерполяции.

k 0

Условия (V.1), наложенные на многочлен, позволяют однозначно определить его коэффициенты. Действительно, требуя для Рn(х) выполнения условий (V.1), получаем линейную систему, состоящую из (n + 1) уравнения:

k 0

Решив систему (V.3) относительно неизвестных а0, а1, … , аn , находим значения этих неизвестных и, подставив в (V.2), находим аналитическое выражение аппроксимирующей функции.

Система (V.4) всегда имеет единственное решение, т. к. ее определитель

известный в алгебре как определитель Вандермонда, отличен от нуля. Следовательно, интерполяционный многочлен Рn(х) существует и единственен.

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

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

»N = 8;

»1=1: N;

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

»y=sin(x);

2. Визуализируйте табличную зависимость и истинные значения функции (рис. V.1).

Рис. V.1. График функции f (x) = sin (x) и табличных значений, используемых для решения задачи интерполяции

110

»М=1000;

»j=l:M;

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

»Y=sin(X) ,

»plot(x,y, 'о');

»hold on

»plot (x,y)

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

Листинг V.1. Файл Pol. M function z=Pol (x,a)

N = length(a); M= length(x); for j=l:M s=0;

for i=l:N

s= s+a(i)*x(j).^(N-i); end;

z(j)=s;

end;

4. Создайте файл Vandermortd.m (листинг V.2), содержащий описание функции, возвращающей значения элементов матрицы Вандермонда.

Листинг V.2. Файл Vandermond. m function z = Vandermond (x)

N = length(x);

z = ones (N,N) ; for i = 1:N

z(i,j) = x(i). ^(N-j); end;

end;

111

5.Вычислите значения элементов матрицы Вандермонда. » M=Vandennond(x) ;

6.Вычислите значения коэффициентов полинома.

» а = М^ - 1*у';

7.Вычислите значения полинома в заданных промежуточных точках

» Y1=Pol (X,a);

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

Рис. V.2. Погрешность аппроксимации функции f (x) = sin (x) полиномом 8-й степени.

112

5.2. Интерполяционный полином Лагранжа

Для функции, заданной табл. V.1, построим интерполяционный многочлен Ln(x), степень которого не выше п и выполнены условия (V.1).

Будем искать Ln (x) в виде

0, если i ≠ к

Очевидно, что требование (V.6) с учетом (V.5) обеспечивают выполнение условий (V.1).

Многочлены li ,(x) составим следующим образом:

li (x)=ci (x –x0)(x - xi) · …· (x – xi-1)(x – xi+1)· …·(x - xn), (V.7)

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

Формула (V.9) решает поставленную задачу.

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

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

»N=8;

»i = 1:N;

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

»y=sin(x);

2. Создайте файл Lagrange.m (листинг V.3), содержащий описание функции возвращающей значение многочлена li ,(x).

function z=Lagrange (x, i , X, Y)

%x - абсцисса точки интерполяции

%I - номер полинома Лагранжа

%X - вектор, содержащий абсциссы узлов интерполяции

%Y - вектор, содержащий ординаты точек интерполяции

N=length(X);

L=l;

for j=l:N

if not(j==i)

L=L*(x-X(j)) / (X(i)-X(j)); end;

end; z=L*Y(i);

3. Создайте файл Pol_Lagr.m (листинг V.4), содержащий описание функции, возвращающей значения полинома Лагранжа.

Листинг V.4 . Файл Lagr. m function z=Pol_I<agr(x,X,Y)

%x - абсцисса точки интерполяции

%i - номер полинома Лагранжа

%X - вектор, содержащий абсциссы узлов интерполяции

%Y - вектор, содержащий ординаты точек интерполяции

N=length(X); s=0,-

for i=l:N

s=s+Lagrange (x , i , X, Y ); end;

z=s;

114

4.Задайте число промежуточных точек, вычислите их координаты и точные значения интерполируемой функции. » М=1000;

» j=1:M;

» X(j) =2*pi / (M-1) * (j -1); » Y = sin (X)

5.Вычислите значения полинома Лагранжа в промежуточных точках.

» for j=l:M

Y2 ( j ) =Pol_Lagr (X ( j ) , x , у) ; end;

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

Рис. V.3. Погрешность аппроксимации функции f (х) = sin (х) полиномом Лагранжа

115

5.3. Интерполяционный полином Ньютона для равноотстоящих узлов

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

5.3.1 . Конечныеразности

Для функции, заданной табл. V.1 с постоянным шагом (V.10), определим разности между значениями функции в соседних узлах интерполяции:

yi = yi +1 – yi . (V.11)

Такое выражение называют конечными разностями первого порядка. Из конечных разностей первого порядка можно образовать конечные разности второго порядка:

2yi = yi +1– yi=(yi +2–yi +1)–(yi +1–yi)=yi +2– 2yi +1+yi. (V.12)

Аналогично получают выражение для конечных разностей третьего порядка:

3yi = 2yi +1 – 2yi = (yi +3 – 2yi +2 + yi +1) – (yi +2 – 2yi +1 + yi) =

5.3.2. Перваяинтерполяционнаяформула Ньютона

Будем искать интерполяционный полином в виде: 116

Значения коэффициентов a0, a1, … , an найдем из условия совпадения значений исходной функции и многочлена в узлах. Полагая х = х0 , из (V.15) найдем у0 = P0(x0) = a0 , откуда а0 = у0. Далее, последовательно придавая х значения x1, и х2, получаем: