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

// The program solves the initial value problem y' = (y - t)/(y + t) (*),
// y(0) = 4 for 0 <= t <= 1 by the four order tree-step interval
// predictor-corrector method of Nystom-Milne-Simpson with the step size
// h = 0.001.
// You can direct output to the screen (s) or to the file (f).
// All starting intervals are defined within the program on the basis of the
// interval Runge-Kutta of order four with step size 0.0001.

// 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 F (const T, Y : interval) : interval;
begin
  Result:=(Y-T)/(Y+T)
end {F};

function F1T (const T, Y : interval) : interval;
var s : interval;
begin
  s:=Y+T;
  s:=2*Y/(s*s);
  Result:=-s
end {F1T};

function F1Y (const T, Y : interval) : interval;
var s : interval;
begin
  s:=Y+T;
  Result:=2*T/(s*s)
end {F1Y};

function F2T2 (const T, Y : interval) : interval;
var s : interval;
begin
  s:=Y+T;
  s:=s*s*s;
  Result:=4*Y/s
end {F2T2};

function F2TY (const T, Y : interval) : interval;
var s : interval;
begin
  s:=Y+T;
  s:=s*s*s;
  Result:=2*(Y-T)/s
end {F2TY};

function F2Y2 (const T, Y : interval) : interval;
var s : interval;
begin
  s:=Y+T;
  s:=s*s*s;
  s:=4*T/s;
  Result:=-s
end {F2Y2};

function F3T3 (const T, Y : interval) : interval;
var s : interval;
begin
  s:=Y+T;
  s:=s*s;
  s:=s*s;
  s:=12*Y/s;
  Result:=-s
end {F3T3};

function F3T2Y (const T, Y : interval) : interval;
var s : interval;
begin
  s:=Y+T;
  s:=s*s;
  s:=s*s;
  s:=(T-2*Y)/s;
  Result:=4*s
end {F3T2Y};

function F3TY2 (const T, Y : interval) : interval;
var s : interval;
begin
  s:=Y+T;
  s:=s*s;
  s:=s*s;
  s:=(2*T-Y)/s;
  Result:=4*s
end {F3TY2};

function F3Y3 (const T, Y : interval) : interval;
var s : interval;
begin
  s:=Y+T;
  s:=s*s;
  s:=s*s;
  Result:=12*T/s
end {F3Y3};

function F4T4 (const T, Y : interval) : interval;
var s, s1 : interval;
begin
  s1:=Y+T;
  s:=s1*s1;
  s:=s1*(s*s);
  Result:=48*Y/s
end {F4T4};

function F4T3Y (const T, Y : interval) : interval;
var s, s1 : interval;
begin
  s1:=Y+T;
  s:=s1*s1;
  s:=s1*(s*s);
  s1:=3*Y-T;
  Result:=12*s1/s
end {F4T3Y};

function F4T2Y2 (const T, Y : interval) : interval;
var s, s1 : interval;
begin
  s1:=Y+T;
  s:=s1*s1;
  s:=s1*(s*s);
  Result:=24*(Y-T)/s
end {F4T2Y2};

function F4TY3 (const T, Y : interval) : interval;
var s, s1 : interval;
begin
  s1:=Y+T;
  s:=s1*s1;
  s:=s1*(s*s);
  s1:=Y-3*T;
  Result:=12*s1/s
end {F4TY3};

function F4Y4 (const T, Y : interval) : interval;
var s, s1 : interval;
begin
  s1:=Y+T;
  s:=s1*s1;
  s:=s1*(s*s);
  s:=48*T/s;
  Result:=-s
end {F4Y4};

function DY2 (const T, Y : interval) : interval;
var S : interval;
begin
  S:=F1Y(T,Y)*F(T,Y);
  Result:=F1T(T,Y)+S
end {DY2};

function DY3 (const T, Y : interval) : interval;
var FTY, S : interval;
begin
  FTY:=F(T,Y);
  S:=F2Y2(T,Y)*FTY;
  S:=(2*F2TY(T,Y)+S)*FTY;
  S:=F2T2(T,Y)+S;
  Result:=S+F1Y(T,Y)*DY2(T,Y)
end {DY3};

function DY4 (const T, Y : interval) : interval;
var FTY, S : interval;
begin
  FTY:=F(T,Y);
  S:=F3Y3(T,Y)*FTY;
  S:=(3*F3TY2(T,Y)+S)*FTY;
  S:=(3*F3T2Y(T,Y)+S)*FTY;
  S:=F3T3(T,Y)+S;
  FTY:=F2TY(T,Y)+F2Y2(T,Y)*FTY;
  S:=S+3*FTY*DY2(T,Y);
  Result:=S+F1Y(T,Y)*DY3(T,Y)
end {DY4};

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

function Psi4Explicit (const T, Y : interval) : interval;
var FTY, F2Y2TY, S, S1, Y2TY : interval;
begin
  FTY:=F(T,Y);
  S:=F4Y4(T,Y)*FTY;
  S:=(4*F4TY3(T,Y)+S)*FTY;
  S:=(6*F4T2Y2(T,Y)+S)*FTY;
  S:=(4*F4T3Y(T,Y)+S)*FTY;
  S:=F4T4(T,Y)+S;
  S1:=F3Y3(T,Y)*FTY;
  S1:=(2*F3TY2(T,Y)+S1)*FTY;
  S1:=F3T2Y(T,Y)+S1;
  Y2TY:=DY2(T,Y);
  S:=S+6*S1*Y2TY;
  F2Y2TY:=F2Y2(T,Y);
  S1:=F2TY(T,Y)+F2Y2TY*FTY;
  S:=S+4*S1*DY3(T,Y);
  S1:=Y2TY*Y2TY;
  S:=S+3*F2Y2TY*S1;
  Result:=S+F1Y(T,Y)*DY4(T,Y)
end {Psi4Explicit};

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

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

procedure Results (const T, Y : interval;
                   const kind : Char);
var 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+']');
  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
  Writeln;
  Writeln ('I N T E R V A L   N Y S T R O M - M I L N E - S I M P S O N   ');
  Writeln;
  Writeln ('M E T H O D   O F   O R D E R   4');
  Writeln;
  repeat
    Write ('Output (s - screen, f - file) : ');
    Readln (z)
  until (z='s') or (z='f');
  if z='f'
    then begin
           AssignFile (plik, 'IPCM5NMSRESULTS.TXT');
           Rewrite (plik);
           Writeln ('Output to IPCM5NMSRESULTS.TXT in current directory');
           Writeln ('Please, wait ...')
         end
    else Writeln;
  m:=1000;
  s:='0.001';
  h.a:=left_read(s);
  h.b:=right_read(s);
  T0:=0;
  Y0:=4;
  T1:=T0+h;
  s:='4.0009997500832942';
  Y1.a:=left_read(s);
  s:='4.0009997500832943';
  Y1.b:=right_read(s);
  T2:=T1+h;
  s:='4.0019990006660423';
  Y2.a:=left_read(s);
  s:='4.0019990006660424';
  Y2.b:=right_read(s);
  if z='f'
    then begin
           Writeln (plik, '   ');
           Writeln (plik, 'I N T E R V A L   N Y S T R O M - M I L N E - '
                    +'S I M P S O M   M E T H O D');
           Writeln (plik, '   ');
           Writeln (plik, 'O F   O R D E R   4');
           Writeln (plik, '   ');
           Writeln (plik, 'h = 0.001');
           Writeln (plik, '   ');
           Writeln (plik, 'Starting intervals:');
           iends_to_strings (Y0, left, right);
           Writeln (plik, 'Y0 = ['+left+','+right+']');
           iends_to_strings (Y1, left, right);
           Writeln (plik, 'Y1 = ['+left+','+right+']');
           iends_to_strings (Y2, left, right);
           Writeln (plik, 'Y2 = ['+left+','+right+']');
           Writeln (plik, '   ');
         end
    else Writeln ('h = 0.001');
  present1:=Now;
  delta_t.a:=0;
  delta_t.b:=1;
  delta_y.a:=4;
  s:='4.9';
  delta_y.b:=right_read(s);
  Fdelta:=F(delta_t,delta_y);
  s:='-0.002';
  hh.a:=left_read(s);
  s:='0.001';
  hh.b:=right_read(s);
  s:='-0.003';
  hh1.a:=left_read(s);
  hh1.b:=0;
  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:=F(T2,Y2);
    i3:=F(T1,Y1);
    i4:=F(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:=F(T3,Y3);
      i3:=F(T2,Y2);
      i2:=i2+4*i3;
      i3:=F(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 100)=0
             then Results (T3, Y4, 'i');
    T0:=T1;
    Y0:=Y1;
    T1:=T2;
    Y1:=Y2;
    T2:=T3;
    Y2:=Y3
  until (n=m) or not find;
  present2:=Now;
  DecodeTime(present2-present1, hour, min, sec, ms);
  Writeln;
  Writeln ('Time of computations: ', min, ' min ', sec, '.', ms, ' sec');
  if z='f'
    then begin
           Writeln (plik, ' ');
           Writeln (plik, 'Time of computations: ', min, ' min ', sec, '.', ms,
                    ' sec');
           CloseFile (plik)
         end;
  Writeln;
  Writeln ('The end. Press Enter ...');
  Readln
end.