program IMM_Example_3;
{$APPTYPE CONSOLE}
uses System.SysUtils, Winapi.Windows, IntervalArithmetic32and64;

// The program solves the initial value problem y' = (y - t)/(y + t), y(0) = 4,
// for delta_t = {t in R : 0 <= t <= 1} and delta_y = {y in R : 4 <= y <= 4.9}
// with h = 0.001 by interval explicit four-step fourth order methods of Milne
// and Nystrom types.
// 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 interval_function = function (const T, Y : interval) : interval;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

function D1FTY (const T, Y : interval) : interval;
begin
  Result:=D1FDT(T,Y)+D1FDY(T,Y)*FTY(T,Y)
end {D1FTY};

function D2FTY (const T, Y : interval) : interval;
var F, F2 : interval;
begin
  F:=FTY(T,Y);
  F2:=(D2FDY2(T,Y)*F+2*D2FDTDY(T,Y))*F;
  Result:=D2FDT2(T,Y)+F2+D1FDY(T,Y)*D1FTY(T,Y)
end {D2FTY};

function D3FTY (const T, Y : interval) : interval;
var F, F2, F3 : interval;
begin
  F:=FTY(T,Y);
  F3:=(D3FDY3(T,Y)*F+3*D3FDTDY2(T,Y))*F;
  F3:=(F3+3*D3FDT2DY(T,Y))*F;
  F3:=D3FDT3(T,Y)+F3;
  F2:=(D2FDTDY(T,Y)+D2FDY2(T,Y)*F)*D1FTY(T,Y);
  Result:=F3+3*F2+D1FDY(T,Y)*D2FTY(T,Y);
end {D3FTY};

function PSI (const T, Y : interval) : interval;
var DF, F, F2, F3, F4 : interval;
begin
  F:=FTY(T,Y);
  F4:=(D4FDY4(T,Y)*F+4*D4FDTDY3(T,Y))*F;
  F4:=(F4+6*D4FDT2DY2(T,Y))*F;
  F4:=(F4+4*D4FDT3DY(T,Y))*F;
  F4:=D4FDT4(T,Y)+F4;
  F3:=(D3FDY3(T,Y)*F+2*D3FDTDY2(T,Y))*F;
  F3:=6*(D3FDT2DY(T,Y)+F3);
  DF:=D1FTY(T,Y);
  F2:=D2FDY2(T,Y);
  F3:=(F3+3*F2*DF)*DF;
  F2:=4*(D2FDTDY(T,Y)+F2*F)*D2FTY(T,Y);
  F3:=F3+F2+D1FDY(T,Y)*D3FTY(T,Y);
  Result:=F4+F3
end {PSI};

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

procedure N4 (const T0, Y0, Y1, Y2, Y3, Fdelta : interval;
              const h                          : Extended;
              const FTY, PSI                   : interval_function;
              out T1, T2, T3, T4, Y4           : interval);
// explicit interval four-step method of Nystrom type
var i1, i2, i3, i4, i5, ih, hh, T, Y : interval;
    s                                : string;
begin
  Str (h:26, s);
  ih.a:=left_read(s);
  ih.b:=right_read(s);
  Str (-3*h:26, s);
  hh.a:=left_read(s);
  hh.b:=ih.b;
  T1:=T0+ih;
  T2:=T1+ih;
  T3:=T2+ih;
  T4:=T3+ih;
  T:=T3+hh;
  Y:=Y3+hh*Fdelta;
  i1:=PSI(T,Y);
  i1:=251*i1-19*i1;
  i2:=ih*ih;
  i2:=i2*i2;
  i2:=i2*ih*i1;
  i1:=i2/720;
  i2:=FTY(T3,Y3);
  i3:=FTY(T2,Y2);
  i4:=FTY(T1,Y1);
  i5:=FTY(T0,Y0);
  i2:=8*i2-5*i3;
  i2:=i2+4*i4-i5;
  i2:=ih*i2/3;
  Y4:=Y2+(i2+i1)
end {N4};

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

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

begin
  Writeln;
  Writeln ('Solving the initial value problem y'' = (y - t)/(y + t), y(0) = 4');
  Writeln;
  Writeln ('for delta_t = {t in R : 0 <= t <= 1}');
  Writeln ('and delta_y = {y in R : 4 <= y <= 4.9}');
  Writeln;
  Writeln ('with h = 0.001');
  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;
  amh:=1;
  real_h:=0.001;
  delta_t.a:=0;
  delta_t.b:=1;
  delta_y.a:=4;
  w:=4.9;
  Str (w:26, s);
  delta_y.b:=right_read(s);
  Fdelta:=FTY(delta_t,delta_y);
  T0.a:=0;
  T0.b:=0;
  Y0.a:=4;
  Y0.b:=4;
  T1.a:=left_read('0.001');
  T1.b:=right_read('0.001');
  Y1.a:=left_read('4.0009997500832942');
  Y1.b:=right_read('4.0009997500832943');
  T2.a:=left_read('0.002');
  T2.b:=right_read('0.002');
  Y2.a:=left_read('4.0019990006660423');
  Y2.b:=right_read('4.0019990006660424');
  T3.a:=left_read('0.003');
  T3.b:=right_read('0.003');
  Y3.a:=left_read('4.0029977522468409');
  Y3.b:=right_read('4.0029977522468410');
  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)           output_step = ');
    Readln (output_step)
  until output_step<=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:='IMM3'+strg+'RESULTS.TXT';
           AssignFile (plik, strg);
           Rewrite (plik);
           Writeln ('Output to '+strg+' in current directory')
         end
    else Writeln;
  for n:=4 to m do
    begin
      if method=1
        then M4 (T0, Y0, Y1, Y2, Y3, Fdelta, real_h, FTY, PSI, T1, T2, T3, T4,
                 Y4)
        else N4 (T0, Y0, Y1, Y2, Y3, Fdelta, real_h, FTY, PSI, T1, T2, T3, T4,
                 Y4);
      if (n mod output_step)=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.