x0
|
x1
|
x2
|
...
|
Xn-1
|
xn
|
y0
|
y1
|
y2
|
...
|
yn-1
|
yn
|
При этом требуется получить значение функции f в точке x,
принадлежащей
отрезку [x0..xn]
но не совпадающей ни с одним значением xi.Часто при этом не известно аналитическое
выражение функции f(x), или оно не пригодно для вычислений.
В этих случаях используется прием построения
приближающей функции F(x), которая очень близка к f(x) и совпадает
с ней в точках x0, x1, x2,...
xn. При этом нахождение приближенной функции называется
интерполяцией, а точки x0,x1,x2,...xn
- узлами интерполяции. Обычно интерполирующую ищут в виде полинома n степени:
Pn(x)=a0xn+a1xn-1+a2xn-2+...+an-1x+an
Для каждого набора точек имеется только один
интерполяционный многочлен, степени не больше n. Однозначно определенный многочлен может быть
представлен в различных видах. Рассмотрим интерполяционный многочлен Ньютона и
Лагранжа.
Интерполяционная формула Лагранжа.
Формула
Лагранжа является наиболее общей, может применяться к таким узлам интерполяции,
что расстояние между соседними узлами не постоянная величина.
Построим
интерполяционный полином Ln(x) степени не больше n, и для
которого выполняются условия Ln(xi)=yi . Запишем его в виде суммы:
Ln(x)=l0(x)+
l1(x)+ l2(x)+...+ ln(x),
(1)
Тогда многочлен lk(x) имеет следующий вид:
lk(x)=
(2)
Подставим
(2) в (1) и перепишем Ln(x)
в виде:
Если
функция f(x), подлежащая интерполяции, дифференцируема
больше чем n+1 раз, то погрешность интерполяции оценивается следующим
образом:
где0<θ<1 (3)
Интерполяционная
формула Ньютона.
Построение
интерполяционного многочлена в форме Ньютона применяется главным образом когда
разность xi+1-xi=h
постоянна для всех значений x=0..n-1.
Конечная
разность k-го порядка:
Δyi=yi+1-yi
Δ2yi=
Δyi+1-
Δyi=yi+2-2yi+1+yi
………………………………
Δkyi=yi+k-kyi+1-k+k(k-1)/2!*yi+k-2+...+(-1)kyi
Будем искать интерполяционный
многочлен в виде:
Pn(x)=a0+a1(x-x0)+a2(x-x0)(x-x1)+...+an(x-x0)(x-x1)...(x-xn-1)
Найдем значения
коэффициентов a0, a1, a2,
...,an:
Полагая x=x0, находим a0=P(x0)=y0;
Далее подставляя
значения x1, x2, ...,xn получаем:
a1=Δy0/h
a2=Δ2y0/2!h2
a3=Δ3y0/3!h3
....................
an=Δny0/n!hn
Таким образом:
Pn(x)=y0+ Δy0/h*(x-x0)+ Δ2y0/2!h2*(x-x0)(x-x1)+...+
Δny0/n!hn*(x-x0)(x-x1)...(x-xn-1) (1)
Практически формула (1)
применяется в несколько ином виде:
Возьмем: t=(x-x0)/h, тогда x=x0+th и формула
(1) переписывается как:
Pn(x)=y0+tΔy0+t(t-1)/2! Δ2y0+...+t(t-1)...(t-n+1)/n!Δny0
(2)
Формула (2) называется
интерполяционной формулой Ньютона.
Погрешность метода
Ньютона оценивается следующим образом:
(3)
Интерполяция
сплайнами.
При большом количестве
узлов интерполяции сильно возрастает степень интерполяционных многочленов, что
делает их неудобными для проведения вычислений.
Высокой степени
многочленов можно избежать, разбив отрезок интерполирования на несколько
частей, с построением в каждой части своего интерполяционного полинома. Такой
метод называется интерполяцией сплайнами. Наиболее распространенным является
построение на каждом отрезке [xi, xi+1],
i=0..n-1 кубической
функции. При этом сплайн – кусочная функция, на каждом отрезке заданная
кубической функцией, является кусочно-непрерывной, вместе со своими первой и
второй производной.
Будем искать кубический
сплайн на каждом из частичных отрезков [xi, xi+1]
в виде:
,
где ai,bi,ci,di – неизвестные.
Из того что Si(xi)=yi
получим:
В силу непрерывности
потребуем совпадения значений в узлах, т.е.:
,i=0..n-1;
(1)
Также потребуем
совпадения значений первой и второй производной:
,i=0..n-2; (2)
,i=0..n-2;
(3)
Из (1) получим n линейных
уравнений с 3n неизвестными
,i=0..n-1; (1*)
Из (2) и (3) получим 2(n-1) линейных
уравнений с теми же неизвестными:
,i=0..n-1;
(2*)
,i=1..n-1;
(3*)
Недостающие два уравнения определим следующим
образом. Предположим, что в точках х0 и хn производная равна нулю и получим еще два уравнения. Получим систему из
3*n линейных уравнений с 3*n
неизвестными. Решим ее любым из методов и построим интерполяционную функцию,
такую что на отрезке [xi, xi+1] она равна Si.
Метод Лагранжа
procedure
TForm1.Button1Click(Sender: TObject);
type tip=array of real;
var x,y:tip;
i,j,n:byte;
p,s,xx:real;
begin
n:=edt.Count;
setlength(x,n);
setlength(y,n);
for i:=0 to n-1 do
x[i]:=edt.massiv[i];edt.Lines.Delete(0);
for i:=0 to n-1 do
y[i]:=edt.massiv[i];edt.Lines.Delete(0);
xx:=strtofloat(edt.Text);
edt.Lines.Delete(0);
s:=0;
for i:=0 to n-1 do
begin
p:=1;
for j:=0 to n-1 do if
i<>j then p:=p*(xx-x[j])/(x[i]-x[j]);
p:=p*y[i];
s:=s+p;
end;
edt.writer('',s,1);
end;
Сплайн –
интерполяция (программа
составляет систему линейных уравнений, решая которую находим коэффициенты
кубических сплайнов).
procedure
TForm1.Button1Click(Sender: TObject);
var b,c,d,x,y:array of real;
urm:array of array of real;
i,j,k,n :byte;
begin
n:=edt.Count;
setlength(x,n);setlength(y,n);
for i:=0 to n-1 do
x[i]:=edt.massiv[i];edt.Lines.Delete(0);
for i:=0 to n-1 do
y[i]:=edt.massiv[i];edt.Lines.Delete(0);
setlength(b,n-1);setlength(c,n-1);setlength(d,n-1);
setlength(urm,3*(n-1),3*(n-1)+1);
for i:=0 to 3*(n-1)-1 do
for j:=0 to 3*(n-1) do
urm[i,j]:=0;
for i:=0 to n-1 do edt.writer('
',y[i],0);
for i:=0 to n-2 do
begin
urm[i,3*i+0]:=x[i+1]-x[i];
urm[i,3*i+1]:=(x[i+1]-x[i])*(x[i+1]-x[i]);
urm[i,3*i+2]:=(x[i+1]-x[i])*(x[i+1]-x[i])*(x[i+1]-x[i]);
urm[i,3*(n-1)]:=y[i+1]-y[i];
end;
for i:=0 to n-3 do
begin
urm[i+n-1,3*i+0]:=1;
urm[i+n-1,3*i+1]:=2*(x[i+1]-x[i]);
urm[i+n-1,3*i+2]:=3*(x[i+1]-x[i])*(x[i+1]-x[i]);
urm[i+n-1,3*i+3]:=-1;
end;
for i:=0 to n-3 do
begin
urm[i+2*n-3,3*i+1]:=1;
urm[i+2*n-3,3*i+2]:=3*(x[i+1]-x[i]);
urm[i+2*n-3,3*i+4]:=-1;
end;
urm[3*n-5,0]:=1;
urm[3*n-5,3*(n-1)]:=0;
urm[3*n-4,3*(n-1)-3]:=1;urm[i+2*n-3,3*(n-1)-2]:=2*(y[n-1]-y[n-2])]
urm[3*n-4,3*(n-1)-1]:=3*(y[n-1]-y[n-2])
*(y[n-1]-y[n-2]);
urm[i+2*n-3,3*(n-1)]:=0
for i:=0 to 3*(n-1)-1 do
begin
edt.writer('',1);
for j:=0 to 3*(n-1) do
edt.writer(' ',urm[i,j],0);
end;
end;
Выполнить
интерполяцию сплайнами третьей степени. Построить график и отметить на нем узлы
интерполяции.
Решение.
Будем
искать кубический сплайн на каждом из частичных отрезков [xi, xi+1],
i=0..2 в виде:
, где ai,bi,ci,di – неизвестные.
Из
того что Si(xi)=yi получим:
В
соответствии с теоретическим положениями изложенными выше, составим систему
линейных уравнений, матрица которой будет иметь вид:
При
этом мы потребовали равенства производной нулю.
Решая
систему уравнений получим вектор решений [b1,c1,d1,b2,c2,d2]:
Подставляя
в уравнение значения b1,c1,d1,
получим на отрезке [7..9]:
Если
выражение упростить то:
или
График имеет вид:
Метод Ньютона
procedure
TForm1.Button1Click(Sender: TObject);
type tip=array of real;
var x,y:tip;
i,j,n:byte;
p,s,xx,t,h:real;
kp:array of array of real;
begin
n:=edt.Count;
setlength(x,n);
setlength(y,n);
for i:=0 to n-1 do
x[i]:=edt.massiv[i];edt.Lines.Delete(0);
for i:=0 to n-1 do
y[i]:=edt.massiv[i];edt.Lines.Delete(0);
xx:=strtofloat(edt.Text);
edt.Lines.Delete(0);
setlength(kp,n,n);
for i:=0 to n-1 do for j:=0 to
n-1 do kp[i,j]:=0;
for i:=0 to n-1 do
kp[0,i]:=y[i];
for i:= 1 to n-1 do
for j:=0 to n-i-1 do
kp[i,j]:=kp[i-1,j+1]-kp[i-1,j];
for i:= 0 to n-1 do
begin
for j:=0 to n-1 do
edt.writer(' ',kp[i,j],0);
edt.writer('',1);
end;
edt.writer('',1);
h:=0.5;
t:=(xx-x[0])/h;
s:=y[0];
for i:=1 to n-1 do
begin
p:=1;
for j:=0 to i-1 do
p:=p*(t-j)/(j+1);
s:=s+p*kp[i,0];
end;
edt.writer('',s,1);;
end;
Построить
интерполяционный многочлен Ньютона. Начертить график и отметить на нем узлы
интерполяции. Вычислить значение функции в точке х=1.25.
xi
|
1
|
1.5
|
2
|
2.5
|
3
|
3.5
|
yi
|
0.5
|
2.2
|
2
|
1.8
|
0.5
|
2.25
|
Решение.
Построим таблицу
конечных разностей в виде матрицы:
Воспользуемся
интерполяционной формулой Ньютона:
Pn(x)=y0+tΔy0+t(t-1)/2! Δ2y0+...+t(t-1)...(t-n+1)/n!Δny0
Подставив значения
получим многочлен пятой степени, упростив который получим:
P5(x)=2.2x5-24x4+101.783x3-20.2x2+211.417x-80.7
Вычислим значение
функции в точке x=1.25; P(1.25)=2.0488;
График функции имеет
вид:
Построить
интерполяционный многочлен Лагранжа. Начертить график и отметить на нем узлы
интерполяции. Вычислить значение в точке х=1.2.
xi
|
0
|
1.25
|
2.125
|
3.25
|
yi
|
5.0
|
4.6
|
5.7
|
5.017
|
4.333
|
Решение.
Построим
интерполяционный многочлен Лагранжа L4(x), подставив значения из таблицы в формулу:
Напишем
программу и вычислим значение многочлена в точке х=1.2:
L4(1.2)=5.657;
Полученный многочлен
имеет четвертую степень. Упростим его и получим:
Построим график
полученного полинома: