Program Splinpr;
{Demonstruje nestabilitu na prikladu porovnani kubickeho splinu
s nespravne zadanymi podminkami (obe na 1 hranici) ve srovnani
se spravnym splinem 1 podminka na obou hranicich}
{Zkusma funkce sin(x) od 0 do 2*pi}
{Vysledky do file SPL.dat v poradi x,y_spravny_spline,y_spatny_spline}
USES SPLIN;
const
  pi = 3.141592;
  nbod = 51;
  ndel = 501;
var
   X,Y, Y2,YPOM2 : RealArrayNP;
   xmax,yp1,ypn,dx,y21,xh,ygs,yfs,dx1: Real;
   i : integer;
   f1 : text;
BEGIN
   ASSIGN(f1,'spl.dat');
   REWRITE(f1);
   xmax := 2.*pi;
   dx := xmax/(nbod - 1);
   FOR i := 1 to nbod DO BEGIN
      x[i] := (i-1)*dx;
      y[i] := sin(x[i])
   END;
   Writeln (' Zadej derivaci y pro x = 0; x=2 pi');
   Read(yp1,ypn);
   Spline(X,Y,nbod,yp1,ypn,Y2);
   writeln('spoctena druha derivace na levem orkraji = ',y2[1]);
   {Writeln (' Pro chybny typ spline zadej 2 deriv. y pro x=0');
   Read(ypom2[1]);}
   Writeln (' Pro chybny typ spline zadame na levem okraji',
           'spoctenou 2 deriv. z predchoz.pripadu');
   Writeln ('Bez zaokrouhlovacich chyb by vysly naprosto stejne vysledky');
   ypom2[1] := y2[1];
   ypom2[2] := 2*(3*((y[2]-y[1])/dx - yp1)/dx-ypom2[1]);
   writeln ('x[1]=', x[1]:7:4,'  ypom2[1] =',ypom2[1]:13:-4,
            '  y2[1] =',y2[1]:13:-4);
   writeln ('x[2]=', x[2]:7:4,'  ypom2[2] =',ypom2[2]:13:-4,
            '  y2[2] =',y2[2]:13:-4);
   FOR i := 3 TO nbod DO BEGIN
       ypom2[i] := 6*(y[i]+y[i-2]-2*y[i-1])/dx/dx -
                   4*ypom2[i-1] - ypom2[i-2];
       writeln(' x[',i:2,']=',x[i]:7:4,'   ypom2[i]= ',ypom2[i]:13:-4,
               '  y2[i]= ',y2[i]:13:-4)
       END;
   dx1 := xmax/(ndel - 1);
   FOR i := 1 to ndel DO BEGIN
       xh := (i-1)*dx1;
       splint(x,y,y2,nbod,xh,ygs);
       splint(x,y,ypom2,nbod,xh,yfs);
       writeln(xh,' ',ygs,' ',yfs);
       writeln(f1,xh:10:6,' ',ygs:13:-4,' ',yfs:13:-4);
       END;
   CLOSE(f1)
END.