program IPCM_Example_8_NMS;
{$APPTYPE CONSOLE}
uses System.SysUtils, IntervalArithmetic32and64;

// The program solves the initial value problem
//    y'[1] = y[2],
//    y'[2] = 5(1 - y[1]^2)y[2] - y[1],     (*)
//    y[1](0) = 2, y[2](0) = 0,
// for 0 <= t <= tmax <= 1 by the interval three-step Nystrom-Milne-Simpson
// predictor-corrector method.
// The value of tmax is determine on the basis of step-size h = 0.0001
// and number of steps m, which you have to enter. All starting intervals
// are determined within the program.
// You can direct output to the screen (s) or to the file (f).

// Important note: Delphi's 64 bit compiler do not support 80-bit Extended
// floating point values on Win64 (Extended = Double on Win64).
// The uTExtendedX87 unit provides for Win64 a replacement FPU-backed 80-bit
// Extended floating point type called TExtendedX87 (see the comments in
// IntervalArithmetic32and64 unit for details).

// Note: Do not compile this program in any Delphi environment older than XE!

type array2interval = array [1..2] of interval;

function FTY (const T : interval;
              const Y : array2interval) : array2interval;
var ys : interval;
begin
  Result[1]:=Y[2];
  ys:=1-Y[1]*Y[1];
  ys:=5*ys*Y[2];
  Result[2]:=ys-Y[1]
end {FTY};

function Psi1Explicit (const T : interval;
                       const Y : array2interval) : array2interval;
var ys, yt : interval;
begin
  yt:=FTY(T,Y)[2];
  Result[1]:=yt;
  ys:=1-Y[1]*Y[1];
  yt:=5*ys*yt;
  ys:=10*(Y[1]*Y[2]);
  ys:=(ys+1)*Y[2];
  Result[2]:=yt-ys
end {Psi1Explicit};

function Psi2Explicit (const T : interval;
                       const Y : array2interval) : array2interval;
var ys, yt : interval;
begin
  yt:=Psi1Explicit(T,Y)[2];
  Result[1]:=yt;
  ys:=1-Y[1]*Y[1];
  yt:=5*ys*yt;
  ys:=30*(Y[1]*Y[2]);
  ys:=(ys+1)*FTY(T,Y)[2];
  yt:=yt-ys;
  ys:=(10*Y[2])*(Y[2]*Y[2]);
  Result[2]:=yt-ys
end {Psi2Explicit};


function Psi3Explicit (const T : interval;
                       const Y : array2interval) : array2interval;
var ys, yt, yu : interval;
begin
  yt:=Psi2Explicit(T,Y)[2];
  Result[1]:=yt;
  ys:=1-Y[1]*Y[1];
  yt:=5*ys*yt;
  ys:=40*(Y[1]*Y[2]);
  ys:=(ys+1)*Psi1Explicit(T,Y)[2];
  yt:=yt-ys;
  yu:=FTY(T,Y)[2];
  ys:=2*(Y[2]*Y[2]);
  ys:=Y[1]*yu+ys;
  ys:=30*ys*yu;
  Result[2]:=yt-ys
end {Psi3Explicit};

function Psi4Explicit (const T : interval;
                       const Y : array2interval) : array2interval;
var ys, yt, yu : interval;
begin
  yt:=Psi3Explicit(T,Y)[2];
  Result[1]:=yt;
  ys:=1-Y[1]*Y[1];
  yt:=5*ys*yt;
  ys:=50*(Y[1]*Y[2]);
  ys:=(ys+1)*Psi2Explicit(T,Y)[2];
  yt:=yt-ys;
  yu:=FTY(T,Y)[2];
  ys:=Y[1]*yu+Y[2]*Y[2];
  ys:=100*ys*Psi1Explicit(T,Y)[2];
  yt:=yt-ys;
  ys:=(150*Y[2])*(yu*yu);
  Result[2]:=yt-ys
end {Psi4Explicit};

function Psi3Implicit (const T : interval;
                       const Y : array2interval) : array2interval;
begin
  Result:=Psi4Explicit(T,Y)
end {Psi3Implicit};

var i, k, l, n, m, max_l            : Integer;
    w                               : Extended;
    delta_t, h, hh, hh1, T, T0, T1,
    T2, T3, T4                      : interval;
    delta_y, i1, i2, i3, i4, Fdelta,
    Y, Y0, Y1, Y2, Y3, Y4           : array2interval;
    left, right                     : string;
    s                               : AnsiString;
    find                            : Boolean;
    kind, z                         : Char;
    plik                            : TextFile;

procedure Results (const T    : interval;
                   const Y    : array2interval;
                   const kind : Char);
var i           : Integer;
    hreal       : Extended;
    hrealstring : string;
begin
  hreal:=(T.a+T.b)/2;
  Str (hreal:4:2, hrealstring);
  iends_to_strings (T, left, right);
  if z='f'
    then Writeln (plik, '   T('+hrealstring+') = ['+left+','+right+']')
    else Writeln ('   T('+hrealstring+') = ['+left+','+right+']');
  for i:=1 to 2 do
    begin
      iends_to_strings (Y[i], left, right);
      if z='f'
        then Writeln (plik, 'Y[', i, ']('+hrealstring+') = ['+left+','
                      +right+']')
        else Writeln ('Y[', i, ']('+hrealstring+') = ['+left+','+right+']');
      w:=int_width(Y[i]);
      if z='f'
        then if (kind='i') and (i=2)
               then Writeln (plik, '     width =', w:11,
                             '               max number of iterations : ',
                             max_l)
               else Writeln (plik, '     width =', w:11)
        else if (kind='i') and (i=2)
               then Writeln ('     width =', w:11,
                             '               max number of iterations : ',
                             max_l)
               else Writeln ('     width =', w:11)
    end
end {Results};

begin
  Writeln;
  Writeln ('THREE-STEP INTERVAL NYSTROM-MILNE-SIMPSON PREDICTOR-CORRECTOR'
           +' METHOD');
  Writeln;
  Writeln ('   Delta_t = {t in R: 0 <= t <= 1}');
  delta_t.a:=0;
  delta_t.b:=1;
  Writeln ('Delta_y[1] = {y[1] in R:  1.8 <= y[1] <= 2.1');
  delta_y[1].a:=left_read('1.8');
  delta_y[1].b:=right_read('2.1');
  Writeln ('Delta_y[2] = {y[2] in R: -0.2 <= y[2] <= 0.1');
  delta_y[2].a:=left_read('-0.2');
  delta_y[2].b:=right_read('0.1');
  Fdelta:=FTY(delta_t,delta_y);
  Writeln;
  Writeln ('Data:');
  Writeln ('            step size h = 0.0001');
  s:='0.0001';
  h.a:=left_read(s);
  h.b:=right_read(s);
  repeat
    Write ('      number of steps m = ');
    Readln (m)
  until m*(h.a+h.b)/2<=1;
  repeat
    Write ('output every k steps  k = ');
    Readln (k)
  until (k>=1) and (k<=m);
  Writeln;
  repeat
    Write ('Output (s - screen, f - file) : ');
    Readln (z)
  until (z='s') or (z='f');
  if z='f'
    then begin
           AssignFile (plik, 'IPCM8NMSRESULTS.TXT');
           Rewrite (plik);
           Writeln ('Output to IPCM8NMSRESULTS.TXT in current directory');
           Writeln (plik, ' ');
           Writeln (plik, 'THREE-STEP INTERVAL NYSTROM-MILNE-SIMPSON PREDICTOR-'
                    +'CORRECTOR METHOD');
           Writeln (plik, '   Delta_t = {t in R: 0 <= t <= 1}');
           Writeln (plik, 'Delta_y[1] = {y[1] in R:  1.8 <= y[1] <= 2.1');
           Writeln (plik, 'Delta_y[2] = {y[2] in R: -0.2 <= y[2] <= 0.1');
           Writeln (plik, ' ');
           Writeln (plik, 'Data:');
           Writeln (plik, '      step size h = 0.0001');
           Writeln (plik, 'number of steps m = ', m);
           Writeln (plik, ' ')
         end
    else Writeln;
  Writeln ('Be patient ... Calculations may take a lot of minutes.');
  Writeln;
  T0:=0;
  Y0[1]:=2;
  Y0[2]:=0;
  T1:=T0+h;
  s:='1.9999999900049981';
  Y1[1].a:=left_read(s);
  s:='1.9999999900049982';
  Y1[1].b:=right_read(s);
  s:='-0.0001998500746398';
  Y1[2].a:=left_read(s);
  s:='-0.0001998500746397';
  Y1[2].b:=right_read(s);
  T2:=T1+h;
  s:='1.9999999600399701';
  Y2[1].a:=left_read(s);
  s:='1.9999999600399702';
  Y2[1].b:=right_read(s);
  s:='-0.0003994005969036';
  Y2[2].a:=left_read(s);
  s:='-0.0003994005969035';
  Y2[2].b:=right_read(s);
  s:='-0.0002';
  hh.a:=left_read(s);
  s:='0.0001';
  hh.b:=right_read(s);
  s:='-0.0003';
  hh1.a:=left_read(s);
  hh1.b:=0;
  n:=2;
  max_l:=0;
  repeat
    n:=n+1;
    T3:=T2+h;
    T:=T2+hh;
    for i:=1 to 2 do
      Y[i]:=Y2[i]+hh*Fdelta[i];
    i1:=Psi3Explicit(T,Y);
    for i:=1 to 2 do
      begin
        i1[i]:=9*i1[i]-i1[i];
        i2[i]:=(h*h)*(h*h);
        i1[i]:=i2[i]*i1[i]/24
      end;
    i2:=FTY(T2,Y2);
    i3:=FTY(T1,Y1);
    i4:=FTY(T0,Y0);
    for i:=1 to 2 do
      begin
        i2[i]:=7*i2[i]-2*i3[i];
        i2[i]:=i2[i]+i4[i];
        i2[i]:=h*i2[i]/3;
        Y3[i]:=Y1[i]+(i2[i]+i1[i])
      end;
    T:=T3+hh1;
//    iterations
    l:=0;  // number of iterations
    repeat
      l:=l+1;
      find:=True;
      for i:=1 to 2 do
        Y[i]:=Y3[i]+hh1*Fdelta[i];
      i1:=Psi3Implicit(T,Y);
      for i:=1 to 2 do
        begin
          i1[i]:=11*i1[i]-19*i1[i];
          i2[i]:=(h*h)*(h*h);
          i2[i]:=i2[i]*h;
          i1[i]:=i2[i]*i1[i]/720
        end;
      i2:=FTY(T3,Y3);
      i3:=FTY(T2,Y2);
      for i:=1 to 2 do
        i2[i]:=i2[i]+4*i3[i];
      i3:=FTY(T1,Y1);
      for i:=1 to 2 do
        begin
          i2[i]:=i2[i]+i3[i];
          i2[i]:=h*i2[i]/3;
          Y4[i]:=Y1[i]+(i2[i]+i1[i])
        end;
      for i:=1 to 2 do
        if (Y3[i].a=0) or (Y3[i].b=0)
          then if (Abs(Y4[i].a-Y3[i].a)>=1e-18) or (Abs(Y4[i].b-Y3[i].b)>=1e-18)
                 then find:=False  // 1e-18 - accuracy
                 else
          else if (Abs((Y4[i].a-Y3[i].a)/Y3[i].a)>=1e-18)
                   or (Abs((Y4[i].b-Y3[i].b)/Y3[i].b)>=1e-18)
                 then find:=False;
      Y3:=Y4;
    until (l=20) or find;  // 20 - max number of iterations
    if not find
      then Writeln ('For n = ', n, ' accuracy not achieved after 20 '
                    +'iterations')
      else if (n mod k)=0
             then Results (T3, Y4, 'i');
    T0:=T1;
    Y0:=Y1;
    T1:=T2;
    Y1:=Y2;
    T2:=T3;
    Y2:=Y3;
    if l>max_l
      then max_l:=l
  until (n=m) or not find;
  if z='f'
    then CloseFile (plik);
  Writeln;
  Writeln ('The end. Press Enter ...');
  Readln
end.