Интерполяция функций (работа 1)

Интерполяция функций

Лабораторная работа по дисциплине «Вычислительные методы линейной алгебры».

Министерство образования Российской Федерации.

Хабаровский государственный Технический Университет.

Кафедра «Прикладная математика и информатика»

Хабаровск 2003

Задание.

1) Построить интерполяционный многочлен Ньютона. Начертить график и отметить на нем узлы интерполяции. Вычислить значения в точке х=1.25.

x>i>

1

1.5

2

2.5

3

3.5

y>i>

0.5

2.2

2

1.8

0.5

2.25

2) Построить интерполяционный многочлен Лагранжа. Начертить график и отметить на нем узлы интерполяции. Вычислить значение в точке х=1.2.

x>i>

0

0.25

1.25

2.125

3.25

y>i>

5.0

4.6

5.7

5.017

4.333

3) Выполнить интерполяцию сплайнами третьей степени. Построить график и отметить на нем узлы интерполяции.

x>i>

7

9

13

y>i>

2

-2

3

Постановка задачи интерполяция.

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

x>0>

x>1>

x>2>

...

X>n-1>

x>n>

y>0>

y>1>

y>2>

...

y>n-1>

y>n>

При этом требуется получить значение функции f в точке x, принадлежащей

отрезку [x>0>..x>n>] но не совпадающей ни с одним значением x>i>.Часто при этом не известно аналитическое выражение функции f(x), или оно не пригодно для вычислений.

В этих случаях используется прием построения приближающей функции F(x), которая очень близка к f(x) и совпадает с ней в точках x>0>, x>1>, x>2>,... x>n>. При этом нахождение приближенной функции называется интерполяцией, а точки x>0>,x>1>,x>2>,...x>n> - узлами интерполяции. Обычно интерполирующую ищут в виде полинома n степени:

P>n>(x)=a>0>xn+a>1>xn-1+a>2>xn-2+...+a>n-1>x+a>n>

Для каждого набора точек имеется только один интерполяционный многочлен, степени не больше n. Однозначно определенный многочлен может быть представлен в различных видах. Рассмотрим интерполяционный многочлен Ньютона и Лагранжа.

Интерполяционная формула Лагранжа.

Формула Лагранжа является наиболее общей, может применяться к таким узлам интерполяции, что расстояние между соседними узлами не постоянная величина.

Построим интерполяционный полином L>n>(x) степени не больше n, и для которого выполняются условия L>n>(x>i>)=y>i> . Запишем его в виде суммы:

L>n>(x)=l>0>(x)+ l>1>(x)+ l>2>(x)+...+ l>n>(x), (1)

где l>k>(x>i>)= y>i>, если i=k, и l>k>(x>i>)= 0, если i≠k;

Тогда многочлен l>k>(x) имеет следующий вид:

l>k>(x)= (2)

Подставим (2) в (1) и перепишем L>n>(x) в виде:

Если функция f(x), подлежащая интерполяции, дифференцируема больше чем n+1 раз, то погрешность интерполяции оценивается следующим образом:

где0<θ<1 (3)

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

Построение интерполяционного многочлена в форме Ньютона применяется главным образом когда разность x>i>>+1>-x>i>=h постоянна для всех значений x=0..n-1.

Конечная разность k-го порядка:

Δy>i>=y>i>>+1>-y>i>

Δ2y>i>= Δy>i>>+1>- Δy>i>=y>i>>+2>-2y>i>>+1>+y>i>

………………………………

Δky>i>=y>i+k>-ky>i+1-k>+k(k-1)/2!*y>i+k-2>+...+(-1)ky>i>

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

Pn(x)=a>0>+a>1>(x-x>0>)+a>2>(x-x>0>)(x-x>1>)+...+a>n>(x-x>0>)(x-x>1>)...(x-x>n-1>)

Найдем значения коэффициентов a>0>, a>1>, a>2>, ...,a>n>:

Полагая x=x>0>, находим a>0>=P(x>0>)=y>0>;

Далее подставляя значения x>1>, x>2>, ...,x>n> получаем:

a>1>=Δy>0>/h

a>2>=Δ2y>0>/2!h2

a>3>=Δ3y>0>/3!h3

....................

a>n>=Δny>0>/n!hn

Таким образом:

Pn(x)=y>0>+ Δy>0>/h*(x-x>0>)+ Δ2y>0>/2!h2*(x-x>0>)(x-x>1>)+...+ Δny>0>/n!hn*(x-x>0>)(x-x>1>)...(x-x>n-1>) (1)

Практически формула (1) применяется в несколько ином виде:

Возьмем: t=(x-x>0>)/h, тогда x=x>0>+th и формула (1) переписывается как:

P>n>(x)=y>0>+tΔy>0>+t(t-1)/2! Δ2y>0>+...+t(t-1)...(t-n+1)/n!Δny>0> (2)

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

Погрешность метода Ньютона оценивается следующим образом:

(3)

Интерполяция сплайнами.

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

Высокой степени многочленов можно избежать, разбив отрезок интерполирования на несколько частей, с построением в каждой части своего интерполяционного полинома. Такой метод называется интерполяцией сплайнами. Наиболее распространенным является построение на каждом отрезке [x>i>, x>i>>+1>], i=0..n-1 кубической функции. При этом сплайн – кусочная функция, на каждом отрезке заданная кубической функцией, является кусочно-непрерывной, вместе со своими первой и второй производной.

Будем искать кубический сплайн на каждом из частичных отрезков [x>i>, x>i>>+1>] в виде:

, где a>i>,b>i>,c>i>,d>i> – неизвестные.

Из того что S>i>(x>i>)=y>i> получим:

В силу непрерывности потребуем совпадения значений в узлах, т.е.:

,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 неизвестными. Решим ее любым из методов и построим интерполяционную функцию, такую что на отрезке [x>i>, x>i>>+1>] она равна S>i>.

Метод Лагранжа

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('',1);

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;

Выполнить интерполяцию сплайнами третьей степени. Построить график и отметить на нем узлы интерполяции.

x>i>

7

9

13

y>i>

2

-2

3

Решение.

Будем искать кубический сплайн на каждом из частичных отрезков [x>i>, x>i>>+1>], i=0..2 в виде:

, где a>i>,b>i>,c>i>,d>i> – неизвестные.

Из того что S>i>(x>i>)=y>i> получим:

В соответствии с теоретическим положениями изложенными выше, составим систему линейных уравнений, матрица которой будет иметь вид:

При этом мы потребовали равенства производной нулю.

Решая систему уравнений получим вектор решений [b>1>,c>1>,d>1>,b>2>,c>2>,d>2>]:

Подставляя в уравнение значения b>1>,c>1>,d>1>, получим на отрезке [7..9]:

Если выражение упростить то:

Аналогично подставляя в уравнение значения b>2>,c>2>,d>2>, получим на отрезке [9..13]:

или

График имеет вид:

Метод Ньютона

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.

x>i>

1

1.5

2

2.5

3

3.5

y>i>

0.5

2.2

2

1.8

0.5

2.25

Решение.

Построим таблицу конечных разностей в виде матрицы:

Воспользуемся интерполяционной формулой Ньютона:

P>n>(x)=y>0>+tΔy>0>+t(t-1)/2! Δ2y>0>+...+t(t-1)...(t-n+1)/n!Δny>0>

Подставив значения получим многочлен пятой степени, упростив который получим:

P>5>(x)=2.2x5-24x4+101.783x3-20.2x2+211.417x-80.7

Вычислим значение функции в точке x=1.25; P(1.25)=2.0488;

График функции имеет вид:

Построить интерполяционный многочлен Лагранжа. Начертить график и отметить на нем узлы интерполяции. Вычислить значение в точке х=1.2.

x>i>

0

0.25

1.25

2.125

3.25

y>i>

5.0

4.6

5.7

5.017

4.333

Решение.

Построим интерполяционный многочлен Лагранжа L>4>(x), подставив значения из таблицы в формулу:

Напишем программу и вычислим значение многочлена в точке х=1.2:

L>4>(1.2)=5.657;

Полученный многочлен имеет четвертую степень. Упростим его и получим:

Построим график полученного полинома: