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

// The program solves the initial value problem y' = (y - t)/(y + t) (*),
// y(0) = 4 for 0 <= t <= tmax by the interval Runge-Kutta method of order 4.
// The value of tmax is determine within the program and you have to enter:
// a such that delta_T = [0, a],
// b_left and b_right such that delta_Y = [b_left, b_right],
// initial step size h0.
// 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!

// The right-hand function in (*) and its derivatives:

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};

// Functions for interval Runge Kutta method of order 4:

function Ki1 (const i    : Integer;
              const T, Y : interval) : interval;
var S : interval;
begin
  if i=1
    then Result:=0
    else begin
           S:=F1Y(T,Y)*F(T,Y);
           case i of
             2,3 : Result:=(F1T(T,Y)+S)/2;
             4   : Result:=S
           end
         end
end {Ki1};

function lambda_i1 (const i    : Integer;
                    const T, Y : interval) : interval;
var S : interval;
begin
  if i=1
    then Result:=0
    else begin
           case i of
             2 : S:=Ki1(1,T,Y)/2;
             3 : S:=Ki1(2,T,Y)/2;
             4 : S:=Ki1(3,T,Y)
           end;
           Result:=S
         end
end {lamba_i1};

function Ki2 (const i    : Integer;
              const T, Y : interval) : interval;
var FTY, S : interval;
begin
  if i=0
    then Result:=0
    else begin
           FTY:=F(T,Y);
           S:=F2Y2(T,Y)*FTY;
           S:=(2*F2TY(T,Y)+S)*FTY;
           S:=F2T2(T,Y)+S;
           if (i=2) or (i=3)
             then S:=S/4;
           Result:=S+2*F1Y(T,Y)*lambda_i1(i,T,Y)
         end
end {Ki2};

function lambda_i2 (const i    : Integer;
                    const T, Y : interval) : interval;
var S : interval;
begin
  if i=1
    then Result:=0
    else begin
           case i of
             2 : S:=Ki2(1,T,Y)/2;
             3 : S:=Ki2(2,T,Y)/2;
             4 : S:=Ki2(3,T,Y)
           end;
           Result:=S
         end
end {lamba_i2};

function Ki3 (const i    : Integer;
              const T, Y : interval) : interval;
var FTY, S, S1 : interval;
begin
  if i=1
    then Result:=0
    else 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;
           if (i=2) or (i=3)
             then S:=S/4;
           S1:=F2TY(T,Y)+F2Y2(T,Y)*FTY;
           S1:=6*S1*lambda_i1(i,T,Y);
           S:=S+S1;
           if (i=2) or (i=3)
             then S:=S/2;
           Result:=S+3*F1Y(T,Y)*lambda_i2(i,T,Y)
         end
end {Ki3};

function lambda_i3 (const i    : Integer;
                    const T, Y : interval) : interval;
var S : interval;
begin
  if i=1
    then Result:=0
    else begin
           case i of
             2 : S:=Ki3(1,T,Y)/2;
             3 : S:=Ki3(2,T,Y)/2;
             4 : S:=Ki3(3,T,Y)
           end;
           Result:=S
         end
end {lamba_i3};

function Ki4 (const i    : Integer;
              const T, Y : interval) : interval;
var FTY, F2Y2TY, Li1, S, S1 : interval;
begin
  if i=1
    then Result:=0
    else 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;
           if (i=2) or (i=3)
             then S:=S/4;
           S1:=F3Y3(T,Y)*FTY;
           S1:=(2*F3TY2(T,Y)+S1)*FTY;
           S1:=F3T2Y(T,Y)+S1;
           Li1:=lambda_i1(i,T,Y);
           S1:=12*S1*Li1;
           S:=S+S1;
           if (i=2) or (i=3)
             then S:=S/4;
           F2Y2TY:=F2Y2(T,Y);
           S1:=F2TY(T,Y)+F2Y2TY*FTY;
           S1:=S1*lambda_i2(i,T,Y);
           case i of
             2,3 : S1:=6*S1;
             4   : S1:=12*S1
           end;
           S:=S+S1;
           S1:=(12*F2Y2TY)*(Li1*Li1);
           S:=S+S1;
           Result:=S+4*F1Y(T,Y)*lambda_i3(i,T,Y)
         end
end {Ki4};

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

function PSI (const T, Y : interval) : interval;
var s : interval;
begin
  s:=(Ki4(1,T,Y)+Ki4(4,T,Y))/6;
  s:=s+(Ki4(2,T,Y)+Ki4(3,T,Y))/3;
  s:=DY5(T,Y)-5*s;
  Result:=s/120
end {PSI};

procedure Find_tmax(const aa, b_left, b_right, M   : Extended;
                    const interval_h0, interval_y0 : interval;
                    out interval_alpha             : interval;
                    out tmax                       : Extended;
                    out tmax_OK                    : Boolean);
var i                                 : Integer;
    alpha, ci_left, ci_right          : Extended;
    a, a1, b, b1, c, c1, d, d1, gamma,
    interval_delta_t, interval_delta_y,
    interval_h4                       : interval;
    ci                                : array [1..4] of interval;
begin
  interval_delta_t.a:=0;
  interval_delta_t.b:=aa;
  interval_delta_y.a:=b_left;
  interval_delta_y.b:=b_right;
  interval_h4:=interval_h0*interval_h0;
  interval_h4:=interval_h4*interval_h4;
  alpha:=M*(interval_h0.a+interval_h0.b)/2;
  interval_alpha.a:=-alpha;
  interval_alpha.b:=alpha;
  tmax:=aa;
  a:=interval_y0;
  b:=interval_delta_y;
  for i:=1 to 4 do
    begin
      case i of
        1   : c:=0;
        2,3 : c:=F(interval_delta_t,interval_delta_y)/2;
        4   : c:=F(interval_delta_t,interval_delta_y)
      end;
      if (c.a>=0) and (c.b>0)
        then begin
               a1.a:=a.b;
               a1.b:=a1.a;
               b1.a:=b.b;
               b1.b:=b1.a;
               c1.a:=c.b;
               c1.b:=c1.a;
               gamma:=(b1-a1)/c1;
               if gamma.a<tmax
                 then tmax:=gamma.a
             end
        else if (c.a<0) and (c.b<=0)
               then begin
                      a1.a:=a.a;
                      a1.b:=a1.a;
                      b1.a:=b.a;
                      b1.b:=b1.a;
                      c1.a:=c.a;
                      c1.b:=c1.a;
                      gamma:=(b1-a1)/c1;
                      if gamma.a<tmax
                        then tmax:=gamma.a
                    end
               else if (c.a<0) and (c.b>0)
                      then begin
                             a1.a:=a.b;
                             a1.b:=a1.a;
                             b1.a:=b.b;
                             b1.b:=b1.a;
                             c1.a:=c.b;
                             c1.b:=c1.a;
                             gamma:=(b1-a1)/c1;
                             if gamma.a<tmax
                               then tmax:=gamma.a;
                             a1.a:=a.a;
                             a1.b:=a1.a;
                             b1.a:=b.a;
                             b1.b:=b1.a;
                             c1.a:=c.a;
                             c1.b:=c1.a;
                             gamma:=(b1-a1)/c1;
                             if gamma.a<tmax
                               then tmax:=gamma.a
                           end
    end;
  for i:=1 to 4 do
    case i of
      1,4 : ci[i]:=F(interval_delta_t,interval_delta_y)/6;
      2,3 : ci[i]:=F(interval_delta_t,interval_delta_y)/3
    end;
  a:=interval_y0;
  d:=PSI(interval_delta_t,interval_delta_y);
  d:=d+interval_alpha;
  d:=d*interval_h4;
  ci_left:=0;
  ci_right:=0;
  for i:=1 to 4 do
    begin
      ci_left:=ci_left+ci[i].a;
      ci_right:=ci_right+ci[i].b
    end;
  if (ci_left>=0) and (ci_right>0)
    then begin
           a1.a:=a.b;
           a1.b:=a1.a;
           b1.a:=b.b;
           b1.b:=b1.a;
           c1.a:=ci_right;
           c1.b:=c1.a;
           d1.a:=d.b;
           d1.b:=d1.a;
           gamma:=(b1-a1-d1)/c1;
           if gamma.a<tmax
             then tmax:=gamma.a
         end
    else if (ci_left<0) and (ci_right<=0)
           then begin
                  a1.a:=a.a;
                  a1.b:=a1.a;
                  b1.a:=b.a;
                  b1.b:=b1.a;
                  c1.a:=ci_left;
                  c1.b:=c1.a;
                  d1.a:=d.a;
                  d1.b:=d1.a;
                  gamma:=(b1-a1-d1)/c1;
                  if gamma.a<tmax
                    then tmax:=gamma.a
                end
           else if (ci_left<0) and (ci_right>0)
                  then begin
                         a1.a:=a.b;
                         a1.b:=a1.a;
                         b1.a:=b.b;
                         b1.b:=b1.a;
                         c1.a:=ci_right;
                         c1.b:=c1.a;
                         d1.a:=d.b;
                         d1.b:=d1.a;
                         gamma:=(b1-a1-d1)/c1;
                         if gamma.a<tmax
                           then tmax:=gamma.a;
                         a1.a:=a.a;
                         a1.b:=a1.a;
                         b1.a:=b.a;
                         b1.b:=b1.a;
                         c1.a:=ci_left;
                         c1.b:=c1.a;
                         d1.a:=d.a;
                         d1.b:=d1.a;
                         gamma:=(b1-a1-d1)/c1;
                         if gamma.a<tmax
                           then tmax:=gamma.a
                       end;
  if tmax<=0
    then begin
           Writeln ('It is impossible to calculate tmax.');
           Writeln ('You should consider to change your data.');
           tmax_OK:=False
         end
    else tmax_OK:=True
end {Find_tmax};

// The end of functions for interval Runge-Kutta method of order 4

var i, k, n                : Integer;
    hour, min, ms, s       : Word;
    at, b_left, b_right, M,
    real_h0, real_h, tmax,
    tmax1, w               : Extended;
    sec                    : Single;
    alpha, h, h0, h05, h5,
    K1, K2, K3, K4,
    T0, Y0, T1, Y1, YY     : interval;
    left, right            : string;
    tmax_OK                : Boolean;
    z                      : Char;
    plik                   : TextFile;
    present1, present2     : TDateTime;

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

begin
  Writeln;
  Writeln ('I N T E R V A L   R U N G E - K U T T A   M E T H O D   O F'
           +'   O R D E R   4');
  Writeln;
  Writeln ('Enter data for finding tmax:');
  Writeln;
  repeat
    Write ('Delta_t = {t in R: 0 <= t <= a} (a <= 4 recommended) a = ');
    Readln (right);
    at:=right_read(right)
  until at>0;
  repeat
    Write ('Delta_y = {y in R: b_left <= y <= b_right}      b_left = ');
    Readln (left);
    b_left:=left_read(left);
    Write ('                                               b_right = ');
    Readln (right);
    b_right:=right_read(right)
  until b_right>b_left;
  repeat
    Write ('Initial step size (0 < h0 < 1)                      h0 = ');
    Readln(left);
    h0:=int_read(left);
    Val (left, real_h0, i)
  until (real_h0>0) and (real_h0<1) and (i=0);
  T0:=0;
  Writeln ('Initial condition                                 y(0) = 4');
  Y0:=4;

// an estimation of M >= Abs((r^(6)_k+1(theta*h))/6!) = Abs(y^(6)(theta*h))/5760,
// where 0 < theta < 1, h <= h0;
// this estimation can be calculated;
// assuming that Abs(y(t)) <= 6.3 for t in [0, 20] we have
// M <= (2.22 + 2.18*h0 + 0.82*h0^2 + 0.13*h0^3 + 0.01*h0^4)*10^(-2),
// what for h0 = 1 yields M <= 5.37*10^(-2)
  M:=((((0.02*real_h0+0.13)*real_h0+0.82)*real_h0+2.18)*real_h0+2.22)*0.01;
  Writeln;
  Writeln ('Constant M calculated for interval RK method    M =', M:11);

// calling the Find_tmax procedure
  Writeln ('Calculating tmax ... Please, wait!');
  Find_tmax (at, b_left, b_right, M, H0, Y0, alpha, tmax1, tmax_OK);
  if tmax_OK
    then Writeln ('Calculated integration interval              tmax =', tmax1);

// selecting the step size h
  Writeln;
  Write ('Your choice of tmax                          tmax = ');
  if tmax_OK
    then repeat
           Readln (tmax);
         until tmax<tmax1
    else Readln (tmax);
  repeat
    Write ('Number of steps                                 n = ');
    Readln (n);
    real_h:=tmax/n
  until real_h<=real_h0;
  Writeln;
  Writeln ('Calculated step size                            h = ', real_h:6:4);
  Str (real_h:26, left);
  h:=int_read(left);
  Writeln;
  repeat
    Write ('Output every k steps                            k = ');
    Readln (k)
  until (k>=1) and (k<=n);
  Writeln;
  repeat
    Write ('Output (s - screen, f - file) : ');
    Readln (z)
  until (z='s') or (z='f');
  if z='f'
    then begin
           AssignFile (plik, 'IPCM5RKRESULTS.TXT');
           Rewrite (plik);
           Writeln ('Output to IPCM5RKRESULTS.TXT in current directory');
           Writeln ('Please, wait ...');
           Writeln (plik, 'I N T E R V A L   R U N G E - K U T T A   '
                    +'M E T H O D   O F   O R D E R   4');
           Writeln (plik, ' ');
           Writeln (plik, '            delta_T = [', 0.0E0, ',', at, ']');
           Writeln (plik, '            delta_Y = [', b_left, ',', b_right, ']');
           Writeln (plik, '                 h0 =  ', real_h0);
           Writeln (plik, '               y(0) =  ', 4.0E0);
           Writeln (plik, '       calculated M =  ', M);
           if tmax_OK
             then Writeln (plik, '    calculated tmax =  ', tmax1);
           Writeln (plik, 'your choice of tmax =  ', tmax);
           Writeln (plik, '                  h =  ', real_h);
           Writeln (plik, ' ')
         end
    else Writeln;
  Writeln ('Be patient ... Calculations may take a lot of minutes.');
  Writeln;
  present1:=Now;
  h05:=h/2;
  h5:=h*h;
  h5:=h5*h5;
  h5:=h*h5;
  for i:=1 to n do
    begin
      K1:=F(T0,Y0);
      T1:=T0+h05;
      Y1:=Y0+h05*K1;
      K2:=F(T1,Y1);
      Y1:=Y0+h05*K2;
      K3:=F(T1,Y1);
      T1:=T0+h;
      Y1:=Y0+h*K3;
      K4:=F(T1,Y1);
      YY:=K1+2*K2;
      YY:=YY+2*K3;
      YY:=h*(YY+K4)/6;
      Y1:=(PSI(T0,Y0)+alpha)*h5;
      Y1:=Y0+YY+Y1;
      if (i mod k)=0
        then Results (T1, Y1);
      if i<n
        then begin
               T0:=T1;
               Y0:=Y1
             end
    end;
  present2:=Now;
  DecodeTime(present2-present1, hour, min, s, ms);
  Writeln;
  Writeln ('Time of computations: ', min, ' min ', s, '.', ms, ' sec');
  if z='f'
    then begin
           Writeln (plik, ' ');
           Writeln (plik, 'Time of computations: ', min, ' min ', s, '.', ms,
                    ' sec');
           CloseFile (plik)
         end;
  Writeln;
  Writeln ('The end. Press Enter ...');
  Readln
end.