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

// The program solves the initial value problem y' = -y, y(0) = 1 by
// tree-step interval preictor-corrector methods of Adams-Bashforth-Moulton type
// and Nystrom-Milne-Simpson type. It is assume h = 0.0005, and all starting
// intervals are defined within the program.
// You can direct output to the screen (s) or to the file (f).
// Directing output to the screen, press Enter after each result displayed.

// 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!

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

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

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

function FTY (const T, Y : interval) : interval;
begin
  Result:=-Y
end {FTY};

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

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

begin
  repeat
    Write ('Output (s - screen, f - file) : ');
    Readln (z)
  until (z='s') or (z='f');
  if z='f'
    then begin
           AssignFile (plik, 'IPCM4LTRESULTS.TXT');
           Rewrite (plik);
           Writeln ('Output to IPCM4LTRESULTS.TXT in current directory')
         end
    else Writeln;
  m:=40000;
  s:='0.0001';
  s:='0.0005';
  h.a:=left_read(s);
  h.b:=right_read(s);
  delta_t.a:=0;
  delta_t.b:=1;
  s:='0.0000000021';
  delta_y.a:=left_read(s);
  delta_y.b:=1;
  Fdelta:=FTY(delta_t,delta_y);
  if z='s'
    then begin
           Writeln ('Initial value problem: y'' = -y, y(0) = 1');
           Writeln;
           Writeln ('Exact solution:');
           for n:=0 to 20 do
             Writeln ('y(', n:2, ') = ', Exp(-n):26);
           Writeln;
           Readln
         end
    else begin
           Writeln (plik, 'Initial value problem: y'' = -y, y(0) = 1');
           Writeln (plik, ' ');
           Writeln (plik, 'Exact solution:');
           Writeln (plik, '===============');
           for n:=0 to 20 do
             Writeln (plik, 'y(', n:2, ') = ', Exp(-n):26)
         end;
// k=3, Adams-Bashforth-Moulton predictor-corrector
  T0:=0;
  Y0:=1;
  T1:=T0+h;
  s:='0.9995001249791692';
  Y1.a:=left_read(s);
  s:='0.9995001249791693';
  Y1.b:=right_read(s);
  T2:=T1+h;
  s:='0.9990004998333749';
  Y2.a:=left_read(s);
  s:='0.9990004998333750';
  Y2.b:=right_read(s);
  s:='-0.001';
  hh.a:=left_read(s);
  s:='0.0005';
  hh.b:=right_read(s);
  s:='-0.0015';
  hh1.a:=left_read(s);
  hh1.b:=0;
  if z='f'
    then begin
           Writeln (plik, '   ');
           Writeln (plik, 'k = 3, Adams-Bashforth-Moulton predictor-corrector');
           Writeln (plik, '==================================================')
         end
    else Writeln ('k = 3, Adams-Bashforth-Moulton predictor-corrector');
  n:=2;
  repeat
    n:=n+1;
    T3:=T2+h;
    T:=T2+hh;
    Y:=Y2+hh*Fdelta;
    i1:=Psi3Explicit(T,Y);
    i2:=(h*h)*(h*h);
    i2:=3*i2*i1;
    i1:=i2/8;
    i2:=FTY(T2,Y2);
    i3:=FTY(T1,Y1);
    i4:=FTY(T0,Y0);
    i2:=23*i2-16*i3;
    i2:=i2+5*i4;
    i2:=h*i2/12;
    Y3:=Y2+(i2+i1);
    T:=T3+hh1;
//    iterations
    l:=0;  // number of iterations
    repeat
      l:=l+1;
      find:=True;
      Y:=Y3+hh1*Fdelta;
      i1:=Psi3Implicit(T,Y);
      i2:=h*h*h;
      i2:=i2*(h*h);
      i2:=19*i2*i1;
      i1:=i2/720;
      i2:=FTY(T3,Y3);
      i3:=FTY(T2,Y2);
      i4:=FTY(T1,Y1);
      i5:=FTY(T0,Y0);
      i2:=9*i2+19*i3;
      i2:=i2-5*i4;
      i2:=i2+i5;
      i2:=h*i2/24;
      Y4:=Y2+(i2-i1);
      if (Y3.a=0) or (Y3.b=0)
        then if (Abs(Y4.a-Y3.a)>=1e-18) or (Abs(Y4.b-Y3.b)>=1e-18)
               then find:=False  // 1e-18 - accuracy
               else
        else if (Abs((Y4.a-Y3.a)/Y3.a)>=1e-18) or (Abs((Y4.b-Y3.b)/Y3.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 1000)=0
             then Results (T3, Y3, 'i');
    T0:=T1;
    Y0:=Y1;
    T1:=T2;
    Y1:=Y2;
    T2:=T3;
    Y2:=Y3
  until (n=m) or not find;
  if z='s'
    then Readln;

// k=3, Nystrom-Milne-simpson predictor-corrector
  T0:=0;
  Y0:=1;
  T1:=T0+h;
  s:='0.9995001249791692';
  Y1.a:=left_read(s);
  s:='0.9995001249791693';
  Y1.b:=right_read(s);
  T2:=T1+h;
  s:='0.9990004998333749';
  Y2.a:=left_read(s);
  s:='0.9990004998333750';
  Y2.b:=right_read(s);
  s:='-0.001';
  hh.a:=left_read(s);
  s:='0.0005';
  hh.b:=right_read(s);
  s:='-0.0015';
  hh1.a:=left_read(s);
  hh1.b:=0;
  if z='f'
    then begin
           Writeln (plik, '   ');
           Writeln (plik, 'k = 3, Nystrom-Milne-Simpson predictor-corrector');
           Writeln (plik, '================================================')
         end
    else Writeln ('k = 3, Nystrom-Milne-Simpson predictor-corrector');
  n:=2;
  repeat
    n:=n+1;
    T3:=T2+h;
    T:=T2+hh;
    Y:=Y2+hh*Fdelta;
    i1:=Psi3Explicit(T,Y);
    i1:=9*i1-i1;
    i2:=(h*h)*(h*h);
    i1:=i2*i1/24;
    i2:=FTY(T2,Y2);
    i3:=FTY(T1,Y1);
    i4:=FTY(T0,Y0);
    i2:=7*i2-2*i3;
    i2:=i2+i4;
    i2:=h*i2/3;
    Y3:=Y1+(i2+i1);
    T:=T3+hh1;
//    iterations
    l:=0;  // number of iterations
    repeat
      l:=l+1;
      find:=True;
      Y:=Y3+hh1*Fdelta;
      i1:=Psi3Implicit(T,Y);
      i1:=11*i1-19*i1;
      i2:=(h*h)*(h*h);
      i2:=i2*h;
      i1:=i2*i1/720;
      i2:=FTY(T3,Y3);
      i3:=FTY(T2,Y2);
      i2:=i2+4*i3;
      i3:=FTY(T1,Y1);
      i2:=i2+i3;
      i2:=h*i2/3;
      Y4:=Y1+(i2+i1);
      if (Y3.a=0) or (Y3.b=0)
        then if (Abs(Y4.a-Y3.a)>=1e-18) or (Abs(Y4.b-Y3.b)>=1e-18)
               then find:=False  // 1e-18 - accuracy
               else
        else if (Abs((Y4.a-Y3.a)/Y3.a)>=1e-18) or (Abs((Y4.b-Y3.b)/Y3.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 1000)=0
             then Results (T3, Y4, 'i');
    T0:=T1;
    Y0:=Y1;
    T1:=T2;
    Y1:=Y2;
    T2:=T3;
    Y2:=Y3
  until (n=m) or not find;
  if z='f'
    then CloseFile (plik);
  Writeln;
  Writeln ('The end. Press Enter ...');
  Readln
end.