Показать сообщение отдельно
Старый 28.02.2013, 10:56   #6 (permalink)
Vladimir_S
Специалист
 
Регистрация: 27.08.2008
Адрес: Санкт-Петербург
Сообщений: 27,807
Сказал(а) спасибо: 340
Поблагодарили 583 раз(а) в 208 сообщениях
Репутация: 113184
По умолчанию

Ну вот, нарисовал. Паскаль (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.
Vladimir_S вне форума   Ответить с цитированием
Ads

Яндекс

Member
 
Регистрация: 31.10.2006
Сообщений: 40200
Записей в дневнике: 0
Сказал(а) спасибо: 0
Поблагодарили 0 раз(а) в 0 сообщениях
Репутация: 55070