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

// The program solves the initial value problem y' = 0.5y, y(0) = 1 for
// 0 <= t <=1 by the tree-step interval predictor-corrector method of
// Nystom-Milne-Simpson type for a number of different step-sizes (0.1, 0.01,
// 0.001, 0.0001, 0.00001, 0.000001).
// You can direct output to the screen (s) or to the file (f).
// All starting intervals are defined within the program.
// 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 Psi1Explicit (const T, Y : interval) : interval;
begin
  Result:=Y/4
end {Psi1Explicit};

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

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

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

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

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

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

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

var i, l, m, n                        : Integer;
    h0, w, t_exact, y_exact           : Extended;
    delta_t, delta_y, i1, i2, i3, i4,
    i5, i6, i7, i8, i9, h, hh, hh1,
    Fdelta, T, Y, T0, Y0, T1, Y1, T2,
    Y2, T3, Y3, T4, Y4, T5, Y5, T6,
    Y6, T7, Y7, T8, Y8                : 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, 'IPCM4RESULTS.TXT');
           Rewrite (plik);
           Writeln ('Output to IPCM4RESULTS.TXT in current directory');
           Writeln ('Please, wait ...')
         end
    else Writeln;
  if z='s'
    then begin
           Writeln ('Exact solution:');
           for i:=1 to 10 do
             Writeln ('y(', 0.01*i:4:2, ') = ', Exp(0.05*i):26);
           Writeln;
           Readln
         end
    else begin
           Writeln (plik, 'Exact solution:');
           Writeln (plik, '===============');
           for i:=1 to 10 do
             Writeln (plik, 'y(', 0.01*i:4:2, ') = ', Exp(0.05*i):26)
         end;
  h0:=0.1;
  m:=10;
  for i:=1 to 6 do
    begin
      if i>1
        then begin
               h0:=h0/10;
               m:=10*m
             end;
      Str (h0:26, s);
      h.a:=left_read(s);
      h.b:=right_read(s);
      delta_t.a:=0;
      delta_t.b:=1;
      delta_y.a:=1;
      s:='1.65';
      delta_y.b:=right_read(s);
      Fdelta:=FTY(delta_t,delta_y);
      T0:=0;
      Y0:=1;
      T1:=T0+h;
      y_exact:=Exp(0.5*h0);
      Str (y_exact:26, s);
      Y1:=int_read(s);
      T2:=T1+h;
      y_exact:=Exp(h0);
      Str (y_exact:26, s);
      Y2:=int_read(s);
      y_exact:=-2*h0;
      Str (y_exact:26, s);
      hh.a:=left_read(s);
      y_exact:=h0;
      Str (y_exact:26, s);
      hh.b:=right_read(s);
      y_exact:=-3*h0;
      Str (y_exact:26, s);
      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, '======================================'
                        +'==========');
               Writeln (plik, 'h0 = ', h0:8:6)
             end
        else begin
               Writeln;
               Writeln ('k = 3, Nystrom-Milne-Simpson predictor-corrector');
               Writeln ('h0 = ', h0:8:6)
             end;
      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 (i=1) or ((i=2) and (n mod 10 =0))
                   or ((i=3) and (n mod 100 = 0))
                   or ((i=4) and (n mod 1000 = 0))
                   or ((i=5) and (n mod 10000 = 0))
                   or ((i=6) and (n mod 100000 = 0))
                 then Results (T3, Y4, 'i');
        T0:=T1;
        Y0:=Y1;
        T1:=T2;
        Y1:=Y2;
        T2:=T3;
        Y2:=Y3
      until (n=m) or not find
  end;
  if z='f'
    then CloseFile (plik);
  Writeln;
  Writeln ('The end. Press Enter ...');
  Readln
end.