[Delphi] n-Order Interpolation (just for self documentation)

Of course, it can be linear, quadratic, cubic or any order interpolation

unit Unit1;

interface

uses
  Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
  Dialogs, ExtCtrls, StdCtrls;

type
  TForm1 = class(TForm)
    Timer1: TTimer;
    Image1: TImage;
    Memo1: TMemo;
    procedure awal;
    procedure datamentah;
    function newton(v:real):real;
    function lagrange(v:real):real;
    procedure diferensiasi;
    procedure interpolasi;
    procedure ilagrange;
    function fx(x:real):real;
    procedure Timer1Timer(Sender: TObject);
  private
    { Private declarations }
  public
    { Public declarations }
  end;

const n=11;
      orde=n;
      skalax=17;
      skalay=5;
var
  Form1: TForm1;
  x,y:array[0..n]of real;
  a,f:array[0..orde]of real;
  x0,y0:integer;
  xi,yi:real;
implementation

{$R *.dfm}

function tform1.lagrange(v:real):real;
var j,m:integer;lag,suku:real;
begin{}
  lag:=0;
  for j:=0 to orde do begin
    suku:=1;
    for m:=0 to orde do begin
      if j<> m then suku:=suku*(v-x[m])/(x[j]-x[m]);
    end;
    lag:=lag+y[j]*suku;
  end;
  lagrange:=lag;
end;

procedure tform1.ilagrange;
var i:integer;v:real;
begin{}
  image1.Canvas.Pen.Color:=clred;
  image1.Canvas.
    MoveTo(x0+round(-n*skalax),y0-round(lagrange(-n)*skalay));
  for i:=-n*100 to n*100 do begin
    v:=i/100;
    image1.Canvas.
      LineTo(x0+round(v*skalax),y0-round(lagrange(v)*skalay));
  end;
end;

procedure tform1.diferensiasi;
var i,j,k:integer;prod:real;
begin
  a[0]:=y[0];
  memo1.Lines.Append('a[0] = '+floattostr(a[0]));
  for i:=1 to orde do begin
    a[i]:=0;
    for j:=0 to i do begin
      prod:=1;
      for k:=0 to i do begin
        if (k<>j) then prod:=prod*(x[j]-x[k]);
      end;
      a[i]:=a[i]+y[j]/prod;
    end;
    memo1.Lines.Append('a['+inttostr(i)+'] = '+floattostr(a[i]));
  end;
end;

function tform1.newton(v:real):real;
var i:integer;newt,suku:real;
begin
  newt:=a[0];
  suku:=1;
  for i:=1 to orde do begin
    suku:=suku*(v-x[i-1]);
    newt:=newt+a[i]*suku;
  end;
  newton:=newt;
end;

procedure tform1.interpolasi;
var i:integer;v:real;
begin
  diferensiasi;
  {plot}
  image1.Canvas.Pen.Color:=cllime;
  image1.Canvas.
    MoveTo(x0+round(-n*skalax),y0-round(newton(-n)*skalay));
  for i:=-n*100 to n*100 do begin
    v:=i/100;
    image1.Canvas.
      LineTo(x0+round(v*skalax),y0-round(newton(v)*skalay));
  end;
end;

procedure tform1.awal;
begin
  randomize;
  image1.Width:=400;                              image1.Height:=400;
  x0:=round(image1.Width/2);                      y0:=round(image1.Height/2);
  image1.Canvas.                                  Pen.Color:=clblue;
  image1.Canvas.MoveTo(0,y0);                     image1.Canvas.LineTo(image1.Width,y0);
  image1.Canvas.MoveTo(x0,0);                     image1.Canvas.LineTo(x0,image1.Height);
  xi:=1.11;                                       memo1.Text:='';
  datamentah;                                     interpolasi;
  //ilagrange;
end;

procedure tform1.datamentah;
var i:integer;
begin
  for i:=0 to n do begin
    x[i]:=i-n/2;                                  y[i]:=fx(x[i]);
    with image1.Canvas                            do begin
      pen.Color:=clblue;                          brush.Color:=clgreen;
      pixels[x0+round(x[i]*skalax),               y0-round(y[i]*skalay)]:=clred;
      Ellipse(x0+round(x[i]*skalax)-7,            y0-round(y[i]*skalay)-7,
              x0+round(x[i]*skalax)+7,            y0-round(y[i]*skalay)+7);
    end;
  end;
end;

function tform1.fx(x:real):real;
begin
  fx:=x*x+random(20)-10;
end;

procedure TForm1.Timer1Timer(Sender: TObject);
begin
  timer1.Enabled:=false;
  awal;
end;

end.