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

Пермский Национальный Исследовательский Политехнический Университет

Электротехнический факультет

Кафедра автоматики и телемеханики

 

Вычислительные методы.

Отчет по лабораторной работе №1.

Интерполирование функций.

Вариант №3.

 

 

Выполнил: студент гр. ТК-11

Ватутин.Н.В.

Проверил: Леготкина Т.С.

 

Пермь 2012 г.

Цель работы:

1. Научиться применять численные методы интерполяции функций, заданных таблично.

2. Изучить метод интерполирования функций с помощью кубических сплайн – функций.

 

Задано:

Задано множество вещественных абсцисс Х0, Х1, …, Хn и соответствующие координаты Y0, Y1, …, Yn. Здесь Х0< Х1< …< Хn, а каждое Yi – некоторое вещественное число, отвечающее Xi, определенное математически или являющееся результатом какого – либо наблюдения.

Задача одномерной интерполяции состоит в построении функции f(x), такой, что f(x_i)=y_i для всех i, т.е. нужно найти многочлен P_n(x) степени не выше n, значения которого в точках X_i (i=0,n) совпадают со значениями данной функции, т.е. P_n(x_i)=y_i. Многочлен P_n(x) является интерполяционным многочленом (полиномом), а точки x_i являются узлами интерполяции.

Задана таблица:

	0	1	2	3	4	5	6	7	8
xi	-4	-3	-2	-1	0	1	2	3	4
yi	0	2	0	0	0	0	0	2	0

Необходимо составить интерполяционный полином Лагранжа и Кубическую сплайн-функцию для нахождения аналитического выражения функции f(x), при этом на экран должны выводиться графики полинома Лагранжа и кубической сплайн - функции.

 

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

Пусть x_i (i=0,n) произвольные узлы, а y_i=f(x_i) – значение функции в узлах. Многочленом степени n, принимающем в точке x_i значения y_i, является интерполяционный многочлен Лагранжа

 

Сплайны:

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

На практике наиболее широкое применение получили сплайны третьей степени, имеющие на заданном отрезке непрерывную, по крайней мере, первую производную. Эти сплайны называются кубическими.

На участке [x₀, x_n] требуется найти функцию ф(х), удовлетворяющую следующим условиям:

 

1.Функция ф(х) на участке [x₀, x_n] должна быть непрерывна вместе со своими первой и второй производными.

 

2. На каждом участке [x_i-1, x_i] функция ф_i(х) должна описываться кубическим полиномом

 

a_oi+a_1i(x-x_i-1)+a_2i(x-x_i-1)²+a_3i(x-x_i-1)³ x с[x_i-1, x_i];

ф_i(x)=

0, x ¢[ x_i-1, x_i] i=1,n

 

Общий сплайн

 

при х с[x₁, x_n].

 

3. Значения сплайна в узлах интерполяции равны значениям заданной функции, т.е. ф(х_i)=y_i, i=0,n.

 

Текст программы на языке Pascal:

Program lab1;

Uses Crt,Graph;

Const

k=8;

xi:Array[0..k]Of Integer=(-4,-3,-2,-1,0,1,2,3,4);

yi:Array[0..k]Of Integer=(0,2,0,0,-3,0,0,2,0);

 

Var Gd,Gm,i,x0,y0,dx,dy:Integer;

x,y,p:Real;

 

Var

a,b,c,d,mz,f,g:Array[0..k]of Real;

Function Lg(x:Real):Real;

Var

i,j:Byte;

P1,P2,f:Real;

 

Begin

f:=0;

For i:=0 to k Do

Begin

P1:=1;

P2:=1;

For j:=0 To k Do

If j<>i Then

Begin

P1:=P1*(X-Xi[j]);

P2:=P2*(Xi[i]-Xi[j]);

End;

f:=f+(P1/P2)*yi[i];

End;

Lg:=F;

End;

 

Procedure Osi;

Var

i:Integer;

s:String[4];

Begin

SetColor(5);

Line(x0,0,x0,GetMaxY);

Line(0,y0,GetMaxX,y0);

i:=x0;

Repeat

Line(i,y0+3,i,y0-3);

Line(GetMaxX-i,y0+3,GetMaxX-i,y0-3);

str((i-x0) div dx,s);

OutTextXY(i-2,y0-16,s);

str(-(i-x0) div dx,s);

OutTextXY(GetMaxX-i-3,y0-16,s);

i:=i+DX;

Until (i-x0) div dx >4;

 

i:=y0;

Repeat

Line(x0-3,i,x0+3,i);

Line(x0-3,GetMaxY-i,x0+3,GetMaxY-i);

If i<>y0 Then

Begin

str(-(i-y0) div dy,s);

OutTextXY(x0+5,i,s);

str((i-y0) div dy,s);

OutTextXY(x0+5,GetMaxY-i,s);

End;

i:=i+Dy;

Until i>=GetMaxY;

End;

 

 

Procedure MakeSplain;

Var

i,j:Byte;

h,m:Array[0..k]of Real;

Begin

 

For i:=1 To k Do

h[i]:=xi[i]-xi[i-1];

g[1]:=2*(h[1]+h[2]);

f[1]:=6*((yi[2]-yi[1])/h[2]-(yi[1]-yi[0])/h[1]);

For i:=2 To k-1 Do

Begin

g[i]:=2*(h[i]+h[i+1])-sqr(h[i])/g[i-1];

f[i]:=6*((yi[i+1]-yi[i])/h[i+1]

-(yi[i]-yi[i-1])/h[i])

-f[i-1]*h[i-1]/g[i-1];

End;

 

 

m[0]:=0;m[k]:=0;

For i:=k-1 DownTo 1 Do

begin

m[i]:=(f[i]-h[i+1]*m[i+1])/g[i];

mz[i]:=m[i];

end;

For i:=1 To k Do

Begin

a[i]:=m[i-1]/(6*h[i]);

b[i]:=m[i]/(6*h[i]);

c[i]:=(yi[i-1]-(m[i-1]*sqr(h[i]))/6)/h[i];

d[i]:=(yi[i]-(m[i]*sqr(h[i]))/6)/h[i];

End;

 

End;

 

Function Splain(x:Real):Real;

Var

N:byte;

s:String[1];

Begin

 

n:=1;

For i:=0 To k-1 Do

If x>=xi[i] Then N:=i+1;

Splain:=a[N]*Sqr(xi[N]-x)*(xi[N]-x)

+b[N]*Sqr(x-xi[N-1])*(x-xi[N-1])

+c[N]*(xi[N]-x)

+d[N]*(x-xi[N-1]);

End;

 

Begin

DetectGraph(Gd,Gm);

InitGraph(Gd,Gm,'D:\bp\bgi');

dx:=GetMaxX Div k;dy:=GetMaxY Div 10;

p:=0.001;

x0:=GetMaxX Div 2;

y0:=GetMaxY Div 2;

Osi;

MakeSplain;

readkey;

x:=-4;

Repeat

y:=Lg(x);

PutPixel(x0+Round(dx*x),y0-Round(dy*y),7);

y:=splain(x);

PutPixel(x0+Round(dx*x),y0-Round(dy*y),9);

x:=x+p;

Until x>4;

ReadKey;

CloseGraph;

 

writeln;

for i:=0 to 8 do

begin

write(' x=',xi[i],', m[',i,']=');

writeln(mz[i]:2:2);

end;

readkey;

writeln;

For i:=1 To k Do

Begin

writeln('splain[',i,']: (',a[i]:2:2,')(x[i]-x)^3+(',b[i]:2:2,')(x-x[i-1])^3+('

,c[i]:2:2,')(x[i]-x)+(',d[i]:2:2,')(x-x[i-1]');

writeln;

End;

readkey;

 

End.

 

Графики полученных функций: