program IMM_Example_5;
{$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 four-step multistep methods of
// Milne or Nystrom types.
// 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;
     interval_function = function (const T : interval;
                                   const Y : array2interval) : array2interval;

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 PSI1 (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 {PSI1};

function PSI2 (const T : interval;
               const Y : array2interval) : array2interval;
var ys, yt : interval;
begin
  yt:=PSI1(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 {PSI2};

function PSI3 (const T : interval;
               const Y : array2interval) : array2interval;
var ys, yt, yu : interval;
begin
  yt:=PSI2(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)*PSI1(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 {PSI3};

function PSI4 (const T : interval;
               const Y : array2interval) : array2interval;
var ys, yt, yu : interval;
begin
  yt:=PSI3(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)*PSI2(T,Y)[2];
  yt:=yt-ys;
  yu:=FTY(T,Y)[2];
  ys:=Y[1]*yu+Y[2]*Y[2];
  ys:=100*ys*PSI1(T,Y)[2];
  yt:=yt-ys;
  ys:=(150*Y[2])*(yu*yu);
  Result[2]:=yt-ys
end {PSI4};

procedure M4 (const T0                     : interval;
              const Y0, Y1, Y2, Y3, Fdelta : array2interval;
              const h                      : Extended;
              const FTY, PSI               : interval_function;
              out T1, T2, T3, T4           : interval;
              out Y4                       : array2interval);
// explicit interval four-step method of Milne type
var s                          : Integer;
    i1, i2, ih, hh, T          : interval;
    FTY1, FTY2, FTY3, PSITY, Y : array2interval;
    strg                       : string;
begin
  Str (h:26, strg);
  ih.a:=left_read(strg);
  ih.b:=right_read(strg);
  Str (-3*h:26, strg);
  hh.a:=left_read(strg);
  hh.b:=ih.b;
  T1:=T0+ih;
  T2:=T1+ih;
  T3:=T2+ih;
  T4:=T3+ih;
  T:=T3+hh;
  for s:=1 to 2 do
    Y[s]:=Y3[s]+hh*Fdelta[s];
  PSITY:=PSI(T,Y);
  FTY3:=FTY(T3,Y3);
  FTY2:=FTY(T2,Y2);
  FTY1:=FTY(T1,Y1);
  for s:=1 to 2 do
    begin
      i1:=27*PSITY[s]-251*PSITY[s];
      i2:=ih*ih;
      i2:=i2*i2;
      i2:=i2*ih*i1;
      i1:=i2/720;
      i2:=2*FTY3[s]-FTY2[s]+2*FTY1[s];
      i2:=4*ih*i2/3;
      Y4[s]:=Y0[s]+(i2-i1)
    end;
end {M4};

procedure N4 (const T0                     : interval;
              const Y0, Y1, Y2, Y3, Fdelta : array2interval;
              const h                      : Extended;
              const FTY, PSI               : interval_function;
              out T1, T2, T3, T4           : interval;
              out Y4                       : array2interval);
// explicit interval four-step method of Nystrom type
var s                                : Integer;
    i1, i2, ih, hh, T                : interval;
    FTY0, FTY1, FTY2, FTY3, PSITY, Y : array2interval;
    strg                             : string;
begin
  Str (h:26, strg);
  ih.a:=left_read(strg);
  ih.b:=right_read(strg);
  Str (-3*h:26, strg);
  hh.a:=left_read(strg);
  hh.b:=ih.b;
  T1:=T0+ih;
  T2:=T1+ih;
  T3:=T2+ih;
  T4:=T3+ih;
  T:=T3+hh;
  for s:=1 to 2 do
    Y[s]:=Y3[s]+hh*Fdelta[s];
  PSITY:=PSI(T,Y);
  FTY3:=FTY(T3,Y3);
  FTY2:=FTY(T2,Y2);
  FTY1:=FTY(T1,Y1);
  FTY0:=FTY(T0,Y0);
  for s:=1 to 2 do
    begin
      i1:=PSI(T,Y)[s];
      i1:=251*PSITY[s]-19*PSITY[s];
      i2:=ih*ih;
      i2:=i2*i2;
      i2:=i2*ih*i1;
      i1:=i2/720;
      i2:=8*FTY3[s]-5*FTY2[s];
      i2:=i2+4*FTY1[s]-FTY0[s];
      i2:=ih*i2/3;
      Y4[s]:=Y2[s]+(i2+i1)
    end;
end {N4};

var k, n, m, method                     : Integer;
    amh, real_h, w                      : Extended;
    delta_t, h, T0, T1, T2, T3, T4      : interval;
    delta_y, Fdelta, Y0, Y1, Y2, Y3, Y4 : array2interval;
    left, right, strg                   : string;
    s                                   : AnsiString;
    find                                : Boolean;
    z                                   : Char;
    plik                                : TextFile;

procedure Results (const T    : interval;
                   const Y    : array2interval);
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 Writeln (plik, '     width =', w:11)
        else Writeln ('     width =', w:11)
    end
end {Results};

begin
  Writeln;
  Writeln ('Solving the initial value problem y''[1] = y[2], y[1](0) = 2');
  Writeln ('                                  y''[2] = 5(1-y[1]^2)y[2]-y[1], '
           +'y[2](0) = 0');
  Writeln;
  Writeln ('for Delta_t = {t in R : 0 <= t <= 1}');
  Writeln ('    Delta_y[1] = {y[1] in R :  1.8 <= y[1] <= 2.1}');
  Writeln ('    Delta_y[2] = {y[2] in R : -0.2 <= y[2] <= 0.1}');
  Writeln;
  Writeln ('with h = 0.0001');
  Writeln;
  Writeln ('by an interval multistep method.');
  Writeln;
  Writeln ('Choose an interval method:');
  Writeln (' 1 - explicit four-step fourth order of Milne type');
  Writeln (' 2 - explicit four-step fourth order of Nystrom type');
  Writeln;
  repeat
    Write ('Your choice of method (1 or 2): ');
    Readln (method)
  until (method=1) or (method=2);
  Writeln;
  Writeln ('You have choosen');
  Write ('the ');
  case method of
     1 : Write ('explicit four-step fourth order method of Milne');
     2 : Write ('explicit four-step fourth order method of Nystrom')
  end;
  Writeln (' type.');
  Writeln;
  delta_t.a:=0;
  delta_t.b:=1;
  delta_y[1].a:=left_read('1.8');
  delta_y[1].b:=right_read('2.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);
  amh:=1;
  real_h:=0.0001;
  delta_t.a:=0;
  delta_t.b:=1;
  w:=1.8;
  Str (w:26, s);
  delta_y[1].a:=left_read(s);
  w:=2.1;
  Str (w:26, s);
  delta_y[1].b:=right_read(s);
  w:=-0.2;
  Str (w:26, s);
  delta_y[2].a:=left_read(s);
  w:=0.1;
  Str (w:26, s);
  delta_y[2].b:=right_read(s);
  Fdelta:=FTY(delta_t,delta_y);
  T0.a:=0;
  T0.b:=0;
  Y0[1].a:=2;
  Y0[1].b:=2;
  Y0[2].a:=0;
  Y0[2].b:=0;
  T1.a:=left_read('0.0001');
  T1.b:=right_read('0.0001');
  Y1[1].a:=left_read('1.9999999900049981E+0000');
  Y1[1].b:=right_read('1.9999999900049982E+0000');
  Y1[2].a:=left_read('-1.9985007463979870E-0004');
  Y1[2].b:=right_read('-1.9985007463979869E-0004');
  T2.a:=left_read('0.0002');
  T2.b:=right_read('0.0002');
  Y2[1].a:=left_read('1.9999999600399701E+0000');
  Y2[1].b:=right_read('1.9999999600399702E+0000');
  Y2[2].a:=left_read('-3.9940059690355807E-0004');
  Y2[2].b:=right_read('-3.9940059690355806E-0004');
  T3.a:=left_read('0.0003');
  T3.b:=right_read('0.0003');
  Y3[1].a:=left_read('1.9999999101348489E+0000');
  Y3[1].b:=right_read('1.9999999101348490E+0000');
  Y3[2].a:=left_read('-5.9865201382483130E-0004');
  Y3[2].b:=right_read('-5.9865201382483129E-0004');
  find:=False;
  repeat
    Write ('maximum nunmber of steps (m*h <= a) m = ');
    Readln (m);
    if (m*real_h<=amh) and (m>=4)
      then find:=True
  until find;
  repeat
    Write ('output every k steps (k <= m)       k = ');
    Readln (k)
  until k<=m;
  Writeln;
  repeat
    Write ('Output (s - screen, f - file) : ');
    Readln (z)
  until (z='s') or (z='f');
  if z='f'
    then begin
           case method of
             1 : strg:='M4';
             2 : strg:='N4'
           end;
           strg:='IMM5'+strg+'RESULTS.TXT';
           AssignFile (plik, strg);
           Rewrite (plik);
           Writeln ('Output to '+strg+' in current directory');
           Writeln (plik, ' ');
           Writeln (plik, 'Solving the initial value problem y''[1] = y[2], '
                    +'y[1](0) = 2');
           Writeln (plik, '                                  y''[2] = '
                    +'5(1-y[1]^2)y[2]-y[1], y[2](0) = 0');
           Writeln (plik, ' ');
           Writeln (plik, 'for 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, 'with h = 0.0001');
           Writeln (plik, ' ');
           Write (plik, 'by an explicit four-step fourth order method ');
           case method of
             1 : Writeln (plik, 'of Milne type');
             2 : Writeln (plik, 'of Nystrom type')
           end;
           Writeln (plik, ' ')
         end;
  Writeln;
  Writeln ('Be patient ... Calculations may take a lot of minutes.');
  Writeln;
  for n:=4 to m do
    begin
      if method=1
        then M4 (T0, Y0, Y1, Y2, Y3, Fdelta, real_h, FTY, PSI4, T1, T2, T3, T4,
                 Y4)
        else N4 (T0, Y0, Y1, Y2, Y3, Fdelta, real_h, FTY, PSI4, T1, T2, T3, T4,
                 Y4);
      if (n mod k)=0
        then Results (T4, Y4);
      if n<m
        then begin
               T0:=T1;
               Y0:=Y1;
               T1:=T2;
               Y1:=Y2;
               T2:=T3;
               Y2:=Y3;
               T3:=T4;
               Y3:=Y4
             end;
    end;
  if z='f'
    then CloseFile (plik);
  Writeln;
  Writeln ('The end. Press Enter ...');
  Readln
end.