Interpolacja Newtona, źle zaprogramowany wzór.

Mój program:

program newton111;
{$APPTYPE CONSOLE}
uses
  SysUtils;

const
max = 1000;
var
n, i, j :integer;
x,y:array[0..max] of double;
ilorazy:array[0..max,0..max] of double;
arg,omega,w:double;

begin
writeln ('podaj liczbe wezlow');
read(n);
n:=n-1;
  
for i:=0 to n do
    begin
    writeln('podaj wezel o indeksie ', i+1);
    readln(x[i]);
    writeln ('podaj wartosc w wezle');
    readln(ilorazy[i,0]);
  
  end;
writeln('podaj argument dla ktorego chcesz obliczyc wartosc wielomianu');
readln(arg);

  for j:=1 to n do
 
   for i:=0 to n-j do
    ilorazy[i,j]:=(ilorazy[i+1,j-1] - ilorazy[i,j-1])/(x[i+j] - x[i]);
    w:=ilorazy[0,0];
       omega:=1;
     
  for i:=1 to n do
        begin
         omega:=(arg-x[i-1]);
         w:=w+ilorazy[0,i]*omega;
         end;

writeln ('wartosc wielomainu w punkcie ', arg,' wynosi ', w);
readln;
end.

Byłem pewien, że liczy wszystko dobrze, gdyż jak sprawdziłem na przykładzie prostej funkcji liniowej dobrze przybliżył mi jej wartość.

Okazało się jednak, że nie liczy do końca poprawnie, gdyż jak się dowiedziałem zaprogramowany jest wzór:

Wn(x)=f(x0) + f(x0,x1) (x-x0) + f(x0,x1,x2) (x-x1) + ... + f(x0,x1,...,xn) (x-x{n-1})
zamiast poprawnego:

Wn(x)=f(x0) + f(x0,x1)(x-x0) + f(x0,x1,x2)(x-x0)(x-x1) + ... + f(x0,x1,...,xn) (x-x0)(x-x1)...(x-x{n-1})
No i utknąłem w tym momencie.

dodanie znacznika <code class="delphi"> i poprawienie tagu "pacsal"- fp

Jakie jest pytanie i dlaczego ten kod jest napisany w tak nieczytelny sposób?

Uczę się pascala od pół roku i to jest pewnie powodem nieczytelności kodu.

Pytaniem jest jak zmodyfikować program aby zaimplementować poprawny wzór czyli aby:
wartości funkcji w poszczególnych węzłach:

f(x0) ... f(x0,x1) ... f(x0,x1,x2) ...f(x0,x1,x2...xn-1)
były mnożone nie przez różnice węzłu dla którego obliczamy przybliżenie z węzłem o następnym indeksie,

(x-x0) ... (x-x1) ... (x-x{n-1})
ale przez tego rodzaju iloczyn:

(x-x0) ... (x-x0)(x-x1) ... (x-x0)(x-x1)(x-x2) ... [(x-x0)(x-x1)...(x-x{n-1})]
właśnie w ten sposób:

Wn(x)=f(x0) + f(x0,x1)(x-x0) + f(x0,x1,x2)(x-x0)(x-x1) + ... + f(x0,x1,...,xn) (x-x0)(x-x1)...(x-x{n-1})
I nad tym problemem siedzę od rana, kombinując na różne sposoby.

EDIT:

Dobra po kolejnych 90 minutach patrzenia się w kod, spłodziłem to:

program newton;

{$APPTYPE CONSOLE}

uses
  SysUtils;

const
   max=1000;

var
   i,j,n : word;
   omega,arg : double;
   x,w : array[0..max] of double;
   ilorazy : array [0..max,0..max] of double;

begin

write('Podaj liczbe wezlow: ');
readln(n);

n:=n-1;


for i:=0 to n do

    begin
    writeln('podaj wezel o indeksie ', i+1);
    readln(x[i]);
    writeln ('podaj wartosc w wezle');
    readln(ilorazy[i,0]);
    end;
   write('Podaj argument, dla ktorych mam obliczyc wartosci wielomianu interpolacyjnego: ');

     readln(arg);


   for j:=1 to n do
     for i:=0 to n-j do
       ilorazy[i,j]:= (ilorazy[i+1,j-1] - ilorazy[i,j-1])/(x[i+j] - x[i]);


   begin
     w[j]:=ilorazy[0,0];

     omega:=1;
     for i:=1 to n do
     begin
       omega:=omega*(arg-x[i-1]);
       w[j]:=w[j]+ilorazy[0,i]*omega;
     end;

     writeln('Wartosc wielomianu interpolacyjnego w punkcie ',arg ,' wynosi: ',w[j]);

end;

readln;
end.

Test funkcją kwadratową wypadł pomyślnie, więc chyba wszystko działa. Prosiłbym jedynie o rzucenie okiem na kod przez kogoś bardziej doświadczonego w programowaniu.

dodanie znacznika <code class="delphi"> - fp