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

// The program solves the initial value problem
//    y'[1] = y[2],
//    y'[2] = 5(1 - y[1]^2)y[2] - y[1],     (*)
//    y[1](0) = 2, y[2](0) = 0,
// for 0 <= t <= tmax by the interval three-stage semi-implicit Runge-Kutta
// method of the fourth order (the Butcher method).
// The value of tmax is determine within the program and you have to enter an
// initial step size h0.
// When tmax is calculated, you will be asked to enter some data for
// calculating the step size h.
// 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 array2Extended = array [1..2] of Extended;
     array2interval = array [1..2] of interval;

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

function F (const T : interval;
            const Y : array2interval) : array2interval;
var s : interval;
begin
  Result[1]:=Y[2];
  s:=1-Y[1]*Y[1];
  s:=5*s*Y[2];
  Result[2]:=s-Y[1]
end {F};

function DF2 (const T : interval;
              const Y : array2interval) : interval;
var s, u : interval;
begin
  s:=F(T,Y)[2];
  u:=1-Y[1]*Y[1];
  s:=5*u*s;
  u:=10*(Y[1]*Y[2]);
  u:=(u+1)*Y[2];
  Result:=s-u
end {DF2};

function D2F2 (const T : interval;
               const Y : array2interval) : interval;
var s, u, v : interval;
begin
  s:=DF2(T,Y);
  u:=1-Y[1]*Y[1];
  s:=5*u*s;
  u:=30*(Y[1]*Y[2]);
  v:=F(T,Y)[2];
  u:=(u+1)*v;
  s:=s-u;
  u:=(10*Y[2])*(Y[2]*Y[2]);
  Result:=s-u
end {D2F2};

function D3F2 (const T : interval;
               const Y : array2interval) : interval;
var s, u, v : interval;
begin
  s:=D2F2(T,Y);
  u:=1-(Y[1]*Y[1]);
  s:=5*u*s;
  u:=40*(Y[1]*Y[2]);
  v:=DF2(T,Y);
  u:=(u+1)*v;
  s:=s-u;
  v:=F(T,Y)[2];
  u:=2*(Y[2]*Y[2]);
  u:=Y[1]*v+u;
  u:=30*(u*v);
  Result:=s-u
end {D3F2};

function D4F2 (const T : interval;
               const Y : array2interval) : interval;
var s, u, v, w : interval;
begin
  s:=D3F2(T,Y);
  u:=1-Y[1]*Y[1];
  s:=5*u*s;
  u:=50*(Y[1]*Y[2]);
  v:=D2F2(T,Y);
  u:=(u+1)*v;
  s:=s-u;
  v:=F(T,Y)[2];
  u:=Y[1]*v+Y[2]*Y[2];
  w:=DF2(T,Y);
  u:=100*(u*w);
  s:=s-u;
  u:=(150*Y[2])*(v*v);
  Result:=s-u
end {D4F2};

function PSI (const T : interval;
              const Y : array2interval) : array2interval;
var fty, df, d2f, d3f, d4f, s, u, v, w, z : interval;
begin
  fty:=F(T,Y)[2];
  df:=DF2(T,Y);
  d2f:=D2F2(T,Y);
  d3f:=D3F2(T,Y);
  d4f:=D4F2(T,Y);
  u:=10*(Y[1]*Y[2])+1;
  v:=5*(Y[1]*Y[1]-1);
  s:=d3f/720;
  w:=u*df+v*d2f;
  Result[1]:=s+w/288;
  s:=19*d4f/2880;
  w:=u*d2f+v*d3f;
  s:=s+w/192;
  w:=Y[1]*Y[2];
  z:=Y[2]*Y[2]-Y[1]*Y[1];
  z:=z-v*w;
  z:=z*df+w*d2f;
  s:=s+5*z/144;
  w:=u-v*v;
  u:=v*u;
  w:=w*d2f-u*df;
  s:=s+w/288;
  u:=2*(Y[1]*df);
  u:=Y[2]*fty+u;
  u:=25*(u*fty);
  Result[2]:=s+u/32
end {PSI};

procedure Find_tmax(const aa                 : Extended;
                    const b_left, b_right, M : array2Extended;
                    const interval_h0        : interval;
                    const interval_y0        : array2interval;
                    out interval_alpha       : array2interval;
                    out tmax                 : Extended;
                    out tmax_OK              : Boolean);
var i, s                              : Integer;
    alpha, ci_left, ci_right          : Extended;
    a, a1, b, b1, c, c1, d, d1, gamma,
    interval_delta_t, interval_h4     : interval;
    interval_delta_y, u               : array2interval;
    ci                                : array [1..4] of interval;
begin
  interval_delta_t.a:=0;
  interval_delta_t.b:=aa;
  for s:=1 to 2 do
    begin
      interval_delta_y[s].a:=b_left[s];
      interval_delta_y[s].b:=b_right[s];
      alpha:=M[s]*(interval_h0.a+interval_h0.b)/2;
      interval_alpha[s].a:=-alpha;
      interval_alpha[s].b:=alpha
    end;
  interval_h4:=interval_h0*interval_h0;
  interval_h4:=interval_h4*interval_h4;
  tmax:=aa;
  for s:=1 to 2 do
    u[s]:=F(interval_delta_t,interval_delta_y)[s];
  for s:=1 to 2 do
    begin
      a:=interval_y0[s];
      b:=interval_delta_y[s];
      for i:=1 to 3 do
        begin
          case i of
            1,3 : c:=0;
            2   : c:=u[s]/2
          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 3 do
        case i of
          1,3 : ci[i]:=u[s]/6;
          2   : ci[i]:=2*u[s]/3
        end;
      a:=interval_y0[s];
      d:=PSI(interval_delta_t,interval_delta_y)[s];
      d:=d+interval_alpha[s];
      d:=d*interval_h4;
      ci_left:=0;
      ci_right:=0;
      for i:=1 to 3 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
    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};

var i, k, l, max_l, n, s                    : Integer;
    at, max_f2, max_df2, max_d2f2, max_d3f2,
    max_d4f2, max_d5f2, real_h0, real_h,
    tmax, tmax1, w, w1, max_y1, max_y2      : Extended;
    h, h0, h05, h025, h5, T0, T1, v         : interval;
    alpha, K1, K2, K2prim, K3, Y0, Y1, YY   : array2interval;
    b_left, b_right, M                      : array2Extended;
    left, right                             : string;
    find, tmax_OK                           : Boolean;
    z                                       : Char;
    plik                                    : TextFile;

procedure Results (const T : interval;
                   const Y : array2interval);
var s           : Integer;
    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+']');
  for s:=1 to 2 do
    begin
      iends_to_strings (Y[s], left, right);
      if z='f'
        then Writeln (plik, 'Y[', s, ']('+hrealstring+') = ['+left+','
                      +right+']')
        else Writeln ('Y[', s, ']('+hrealstring+') = ['+left+','+right+']');
      w:=int_width(Y[s]);
      if z='f'
        then if s=1
               then Writeln (plik, '     width =', w:11)
               else Writeln (plik, '     width =', w:11,
                             '               max number of iterations = ',
                             max_l)
        else if s=1
               then Writeln ('     width =', w:11)
               else Writeln ('     width =', w:11,
                             '               max number of iterations = ',
                             max_l)
    end
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 ('Three-stage semi-implicit method of Butcher');
  Writeln;
  Writeln ('   Delta_t = {t in R: 0 <= t <= 1}');
  at:=1;
  Writeln ('Delta_y[1] = {y[1] in R:  1.8 <= y[1] <= 2.1');
  b_left[1]:=left_read('1.8');
  b_right[1]:=right_read('2.1');
  Writeln ('Delta_y[2] = {y[2] in R: -0.2 <= y[2] <= 0.1');
  b_left[2]:=left_read('-0.2');
  b_right[2]:=right_read('0.1');
  Writeln;
  repeat
    Write ('Enter 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 conditions:              y[1](0) = 2');
  Y0[1]:=2;
  Writeln ('                                 y[2](0) = 0');
  Y0[2]:=0;

// estimations of M[s],
// where M[s] >= Abs((r^(6)_k+1[s](theta*h))/6!)
//             = Abs(y^(6)[s](theta*h))/5760 = 7*Abs(f^(5)[s](theta*h))/5760,
// where s = 0, 1 and 0 < theta < 1, h <= h0;
// these estimations can be calculated:
  max_y1:=b_right[1];
  max_y2:=Abs(b_left[2]);
  w:=5*(1+Sqr(max_y1));
  max_f2:=w*max_y2+max_y1;
  w1:=max_y1*max_y2;
  max_df2:=w*max_f2+(10*w1+1)*max_y2;
  max_d2f2:=w*max_df2+(30*w1+1)*max_f2+10*max_y2*Sqr(max_y2);
  max_d3f2:=w*max_d2f2+(40*w1+1)*max_df2
            +30*(max_y1*max_f2+2*Sqr(max_y2))*max_f2;
  max_d4f2:=w*max_d3f2+(50*w1+1)*max_d2f2
            +100*(max_y1*max_f2+Sqr(max_y2))*max_df2+150*max_y2*Sqr(max_f2);
  M[1]:=7*max_d4f2/5760;
  max_d5f2:=w*max_d4f2+(60*w1+1)*max_d3f2
            +150*(max_y1*max_f2+Sqr(max_y2))*max_d2f2
            +100*(6*max_y2*max_f2+max_y1*max_df2)*max_df2
            +150*max_f2*Sqr(max_f2);
  M[2]:=7*max_d5f2/5760;
  Writeln;
  Writeln ('Constants M[s] calculated for RK method:  M[1] =', M[1]:11);
  Writeln ('                                          M[2] =', M[2]: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, 'IPCM8RKRESULTS.TXT');
           Rewrite (plik);
           Writeln ('Output to IPCM8RKRESULTS.TXT in current directory');
           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, 'Three-stage semi-implicit method of Butcher');
           Writeln (plik, ' ');
           Writeln (plik, '            delta_T = [', 0.0E0, ',', at, ']');
           for s:=1 to 2 do
             Writeln (plik, '         delta_Y[', s, '] = [', b_left[s], ',',
                      b_right[s], ']');
           Writeln (plik, '                 h0 =  ', real_h0);
           Writeln (plik, '            y[1](0) =  ', 2.0E0);
           Writeln (plik, '            y[2](0) =  ', 0.0E0);
           for s:=1 to 2 do
             Writeln (plik, '    calculated M[', s, '] =  ', M[s]);
           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;
  h05:=h/2;
  h025:=h/4;
  h5:=h*h;
  h5:=h5*h5;
  h5:=h*h5;
  i:=0;
  max_l:=0;
  repeat
    i:=i+1;
    K1:=F(T0,Y0);
    T1:=T0+h05;
    K2:=F(T1,Y0);
    l:=0;
    repeat
      l:=l+1;
      find:=True;
      for s:=1 to 2 do
        begin
          v:=h025*(K1[s]+K2[s]);
          Y1[s]:=Y0[s]+v
        end;
      K2prim:=F(T1,Y1);
      if (K2[1].a=0) or (K2[1].b=0) or (K2[2].a=0) or (K2[2].b=0)
        then if (Abs(K2prim[1].a-K2[1].a)>=1e-18)
                 or (Abs(K2prim[1].b-K2[1].b)>=1e-18)
                 or (Abs(K2prim[2].a-K2[2].a)>=1e-18)
                 or (Abs(K2prim[2].b-K2[2].b)>=1e-18)
               then find:=False  // 1e-18 - accuracy
               else
        else if (Abs((K2prim[1].a-K2[1].a)/K2[1].a)>=1e-18)
                 or (Abs((K2prim[1].b-K2[1].b)/K2[1].b)>=1e-18)
                 or (Abs((K2prim[2].a-K2[2].a)/K2[2].a)>=1e-18)
                 or (Abs((K2prim[2].b-K2[2].b)/K2[2].b)>=1e-18)
               then find:=False;
      K2:=K2prim
    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 begin
             for s:=1 to 2 do
               Y1[s]:=Y0[s]+h*K2[s];
             K3:=F(T0,Y1);
             for s:=1 to 2 do
               begin
                 v:=K1[s]+4*K2[s];
                 v:=v+K3[s];
                 v:=h*v/6;
                 Y1[s]:=Y0[s]+v;
                 v:=PSI(T0,Y0)[s];
                 v:=(v+alpha[s])*h5;
                 Y1[s]:=Y1[s]+v
               end;
             T1:=T0+h;
             if l>max_l
               then max_l:=l;
             if (i mod k)=0
               then Results (T1, Y1);
             if i<n
               then begin
                      T0:=T1;
                      Y0:=Y1
                    end
           end;
  until (i=n) or not find;
  if z='f'
    then CloseFile (plik);
  Writeln;
  Writeln ('The end. Press Enter ...');
  Readln
end.