Ну вот, нарисовал. Паскаль (Turbo, Free), естественно. Метод решения был взят отсюда (да простят меня коллеги за ссылки, но без них - никак):
Многочлен Ньютона интерполяционный
Лагранжа интерполяционная формула
Между прочим, в первой (Ньютон) лёгкий брёх: поскольку нумерация начинается с нуля, то считать надо не до
n+1, как у них, а до
n. Учтено.
В программе по запросам предусмотрены следующие ветвления: расчет интерполяции функции пользователя (в этом случае задаются начальное и конечное значения интервала, степень полинома
n и
(n+1) значений функции), работа с тестовой функцией (
Sin(x) от 0 до 6), расчет интерполяционного полинома каждым из методов отдельно либо обоими сразу.
С Днём рождения, Никита!
Код:
Var
x,y,b:Array[0..9] of Real;
n,m,i,j,k,N_L,T_R:Byte;
Xbeg,Xfin,Dx,DX1,P,P1,P2,S,xx,yy:Real;
Begin
Write('Your function - 1, Test function (Sin(x) from 0 to 6) - 2 ');
Readln(T_R);
If T_R=1 then
begin
Write('Chooze the method: Newton - 1, Lagrange - 2, Both - 3 ');
Readln(N_L);
end
else N_L:=3;
If T_R=1 then
begin
Write('Polinomial degree (n<10) = ');
Readln(n);
Write('Xbeg = ');
Readln(Xbeg);
Write('Xfin = ');
Readln(Xfin);
DX:=(Xfin-Xbeg)/n;
for i:=0 to n do x[i]:=Xbeg+DX*i;
for i:=0 to n do
begin
write('y[',i,']= ');
readln(y[i]);
end;
end
else
begin
n:=6;
Xbeg:=0;
Xfin:=6;
DX:=1;
b[0]:=y[0];
for i:=0 to n do
begin
x[i]:=Xbeg+DX*i;
y[i]:=Sin(x[i]);
end;
end;
If (N_L=1) or (N_L=3) then {Newton coefficients}
begin
for i:=1 to n do
begin
S:=b[0];
for j:=1 to i-1 do
begin
P:=1;
for k:=0 to j-1 do
P:=P*(x[i]-x[k]);
S:=S+b[j]*P;
end;
P:=1;
for k:=0 to i-1 do
P:=P*(x[i]-x[k]);
b[i]:=(y[i]-S)/P;
end;
end;
If T_R=1 then
begin
Write('New n: ');
Readln(m);
end else m:=n*2;
DX1:=(Xfin-Xbeg)/m;
Write(' x ');
If (N_L=1) or (N_L=3) then Write(' Newton ');
If (N_L=2) or (N_L=3) then Write(' Lagrange ');
If T_R=2 then Writeln(' Accurate ') else Writeln;
For i:=0 to m do
begin
xx:=Xbeg+DX1*i;
Write(xx:7:3);
if (N_L=1) or (N_L=3) then { Newton }
begin
yy:=b[0];
for j:=1 to n do
begin
P:=1;
for k:=0 to j-1 do P:=P*(xx-x[k]);
yy:=yy+b[j]*P;
end;
write(yy:13:3);
end;
if (N_L=2) or (N_L=3) then { Lagrange }
begin
yy:=0;
for j:=0 to n do
begin
P1:=1;
P2:=1;
for k:=0 to n do
if k<>j then
begin
P1:=P1*(xx-x[k]);
P2:=P2*(x[j]-x[k]);
end;
yy:=yy+y[j]*P1/P2;
end;
write(yy:13:3);
end;
if T_R=2 then writeln(Sin(xx):13:3) else writeln;
end;
Readln;
End.