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

// The program solves the simple pendulum motion problem given by the initial
// value problem y'[1] = -u^2*y[2], y'[2] = y[1], y[1](0) = 0, y[2] = phi0,
// where u = Sqrt(g/L), g = 9.80665 , L = 1, phi0 = Pi/6, by two-step
// interval preictor-corrector method of Adams-Bashforth-Moulton type
// and Nystrom-Milne-Simpson type. It is assume h = 0.001, and all starting
// intervals are defined within the program.
// You can direct output to the screen (s) or to the file (f).
// 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!

type array2interval = array [1..2] of interval;

function Psi2Explicit (const u2, T : interval;
                       const Y     : array2interval) : array2interval;
var u4 : interval;
begin
  u4:=u2*u2;
  Result[1]:=u4*Y[2];
  Result[2]:=-u2*Y[1]
end {Psi2Explicit};

function Psi3Explicit (const u2, T : interval;
                       const Y     : array2interval) : array2interval;
var u4 : interval;
begin
  u4:=u2*u2;
  Result[1]:=u4*Y[1];
  Result[2]:=u4*Y[2]
end {Psi3Explicit};

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

function FTY (const u2, T : interval;
              const Y     : array2interval) : array2interval;
begin
  Result[1]:=-u2*Y[2];
  Result[2]:=Y[1]
end {FTY};

var i,l, n, m                       : Integer;
    w                               : Extended;
    delta_t, h, hh, hh1, T, T0, T1,
    T2, T3, T4, u2                  : Interval;
    i1, i2, i3, i4, delta_y, Fdelta,
    Y, Y0, Y1, Y2, Y3, Y4           : array2interval;
    left, right, s                  : string;
    find                            : Boolean;
    z                               : Char;
    plik                            : TextFile;

procedure Results (const T    : interval;
                   const Y    : array2interval;
                   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+']');
  if z='f'
    then begin
           iends_to_strings (Y[1], left, right);
           Writeln (plik, 'Y[1]('+hrealstring+') = ['+left+','+right+']');
           w:=int_width(Y[1]);
           Writeln (plik, 'width = ', w:11);
           iends_to_strings (Y[2], left, right);
           Writeln (plik, 'Y[2]('+hrealstring+') = ['+left+','+right+']');
           w:=int_width(Y[2]);
           Writeln (plik, 'width = ', w:11);
           Writeln (plik, 'number of iterations : ', l)
         end
    else begin
           iends_to_strings (Y[1], left, right);
           Writeln ('Y[1]('+hrealstring+') = ['+left+','+right+']');
           w:=int_width(Y[1]);
           Writeln ('width = ', w:11);
           iends_to_strings (Y[2], left, right);
           Writeln ('Y[2]('+hrealstring+') = ['+left+','+right+']');
           w:=int_width(Y[2]);
           Writeln ('width = ', w:11);
           Writeln ('number of iterations : ', l)
         end
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, 'IPCM6RESULTS.TXT');
           Rewrite (plik);
           Writeln ('Output to IPCM6RESULTS.TXT in current directory')
         end
    else Writeln;
  m:=2000;
  s:='9.80665';
  u2.a:=left_read(s);
  u2.b:=right_read(s);
  s:='0.001';
  h.a:=left_read(s);
  h.b:=right_read(s);
  delta_t.a:=0;
  delta_t.b:=2;
  s:='-1.8';
  delta_y[1].a:=left_read(s);
  s:='1.8';
  delta_y[1].b:=right_read(s);
  s:='-0.6';
  delta_y[2].a:=left_read(s);
  s:='0.6';
  delta_y[2].b:=right_read(s);
  Fdelta:=FTY(u2,delta_t,delta_y);
  w:=Sqrt(9.80665);
  if z='s'
    then begin
           Writeln ('Initial value problem: y''[1] = -u^2*y[2], y''[2] = '
                    +'y[1],');
           Writeln ('                       y[1](0) = 0,       y[2](0) = '
                    +'phi0,');
           Writeln ('where u = Sqrt(g/L), g = 9.80665, L = 1, phi0 = Pi/6');
           Writeln;
           Writeln ('Exact solution:');
           Writeln ('y[1](0.5) = ', -Pi*w*sin(0.5*w)/6:26,
                    '  y[2](0.5) = ', Pi*cos(0.5*w)/6:26);
           Writeln ('y[1](1.0) = ', -Pi*w*sin(w)/6:26,
                    '  y[2](1.0) = ', Pi*cos(w)/6:26);
           Writeln ('y[1](1.5) = ', -Pi*w*sin(1.5*w)/6:26,
                    '  y[2](1.5) = ', Pi*cos(1.5*w)/6:26);
           Writeln ('y[1](2.0) = ', -Pi*w*sin(2*w)/6:26,
                    '  y[2](2.0) = ', Pi*cos(2*w)/6:26);
           Writeln;
           Readln
         end
    else begin
           Writeln (plik, 'Initial value problem: y''[1] = -u^2*y[2], y''[2] = '
                    +'y[1],');
           Writeln (plik, '                       y[1](0) = 0,       y[2](0) = '
                    +'phi0,');
           Writeln (plik, 'where u = Sqrt(g/L), g = 9.80665, L = 1, '
                    +'phi0 = Pi/6');
           Writeln (plik, ' ');
           Writeln (plik, 'Exact solution:');
           Writeln (plik, '===============');
           Writeln (plik, 'y[1](0.5) = ', -Pi*sin(0.5*w)/6:26,
                    '  y[2](0.5) = ', Pi*cos(0.5*w)/6:26);
           Writeln (plik, 'y[1](1.0) = ', -Pi*sin(w)/6:26,
                    '  y[2](1.0) = ', Pi*cos(w)/6:26);
           Writeln (plik, 'y[1](1.5) = ', -Pi*sin(1.5*w)/6:26,
                    '  y[2](1.5) = ', Pi*cos(1.5*w)/6:26);
           Writeln (plik, 'y[1](2.0) = ', -Pi*sin(2*w)/6:26,
                    '  y[2](2.0) = ', Pi*cos(2*w)/6:26)
         end;

// k=2, Adams-Bashforth-Moulton predictor-corrector
  T0:=0;
  Y0[1]:=0;
  s:='0.52359877559829887';
  Y0[2].a:=left_read(s);
  s:='0.52359877559829888';
  Y0[2].b:=right_read(s);
  T1:=T0+h;
  s:='-0.00513474154027124';
  Y1[1].a:=left_read(s);
  s:='-0.00513474154027123';
  Y1[1].b:=right_read(s);
  s:='0.52359620822533077';
  Y1[2].a:=left_read(s);
  s:='0.52359620822533078';
  Y1[2].b:=right_read(s);
  s:='-0.001';
  hh.a:=left_read(s);
  s:='0.001';
  hh.b:=right_read(s);
  s:='-0.002';
  hh1.a:=left_read(s);
  hh1.b:=0;
  if z='f'
    then begin
           Writeln (plik, '   ');
           Writeln (plik, 'k = 2, Adams-Bashforth-Moulton predictor-corrector');
           Writeln (plik, '==================================================')
         end
    else Writeln ('k = 2, Adams-Bashforth-Moulton predictor-corrector');
  n:=1;
  repeat
    n:=n+1;
    T2:=T1+h;
    T:=T1+hh;
    for i:=1 to 2 do
      Y[i]:=Y1[i]+hh*Fdelta[i];
    i1:=Psi2Explicit(u2,T,Y);
    for i:=1 to 2 do
      begin
        i2[i]:=h*h*h;
        i2[i]:=5*i2[i]*i1[i];
        i1[i]:=i2[i]/12
      end;
    i2:=FTY(u2,T1,Y1);
    i3:=FTY(u2,T0,Y0);
    for i:=1 to 2 do
      begin
        i2[i]:=3*i2[i]-i3[i];
        i2[i]:=h*i2[i]/2;
        Y2[i]:=Y1[i]+(i2[i]+i1[i])
      end;
    T:=T2+hh1;
//    iterations
    l:=0;  // number of iterations
    repeat
      l:=l+1;
      find:=True;
      for i:=1 to 2 do
        Y[i]:=Y2[i]+hh1*Fdelta[i];
      i1:=Psi2Implicit(u2,T,Y);
      for i:=1 to 2 do
        begin
          i2[i]:=(h*h)*(h*h);
          i1[i]:=i2[i]*i1[i]/24
        end;
      i2:=FTY(u2,T2,Y2);
      i3:=FTY(u2,T1,Y1);
      i4:=FTY(u2,T0,Y0);
      for i:=1 to 2 do
        begin
          i2[i]:=5*i2[i]+8*i3[i];
          i2[i]:=i2[i]-i4[i];
          i2[i]:=h*i2[i]/12;
          Y3[i]:=Y1[i]+(i2[i]-i1[i]);
          if (Y2[i].a=0) or (Y2[i].b=0)
            then if (Abs(Y3[i].a-Y2[i].a)>=1e-18)
                     or (Abs(Y3[i].b-Y2[i].b)>=1e-18)
                   then find:=False  // 1e-18 - accuracy
                   else
            else if (Abs((Y3[i].a-Y2[i].a)/Y2[i].a)>=1e-18)
                     or (Abs((Y3[i].b-Y2[i].b)/Y2[i].b)>=1e-18)
                   then find:=False
        end;
      Y2:=Y3
    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=500) or (n=1000) or (n=1500) or (n=2000)
             then Results (T2, Y2, 'i');
    T0:=T1;
    Y0:=Y1;
    T1:=T2;
    Y1:=Y2
  until (n=m) or not find;
  if z='s'
    then Readln;

// k=2, Nystrom-Milne-Simpson predictor-corrector
  T0:=0;
  Y0[1]:=0;
  s:='0.52359877559829887';
  Y0[2].a:=left_read(s);
  s:='0.52359877559829888';
  Y0[2].b:=right_read(s);
  T1:=T0+h;
  s:='-0.00513474154027124';
  Y1[1].a:=left_read(s);
  s:='-0.00513474154027123';
  Y1[1].b:=right_read(s);
  s:='0.52359620822533077';
  Y1[2].a:=left_read(s);
  s:='0.52359620822533078';
  Y1[2].b:=right_read(s);
  s:='-0.001';
  hh.a:=left_read(s);
  s:='0.001';
  hh.b:=right_read(s);
  s:='-0.002';
  hh1.a:=left_read(s);
  hh1.b:=0;
  if z='f'
    then begin
           Writeln (plik, '   ');
           Writeln (plik, 'k = 2, Nystrom-Milne-Simpson predictor-corrector');
           Writeln (plik, '================================================')
         end
    else Writeln ('k = 2, Nystrom-Milne-Simpson predictor-corrector');
  n:=1;
  repeat
    n:=n+1;
    T2:=T1+h;
    T:=T1+hh;
    for i:=1 to 2 do
      Y[i]:=Y1[i]+hh*Fdelta[i];
    i1:=Psi2Explicit(u2,T,Y);
    for i:=1 to 2 do
      begin
        i1[i]:=5*i1[i]-i1[i];
        i2[i]:=h*h*h;
        i2[i]:=i2[i]*i1[i]/12
      end;
    i1:=FTY(u2,T1,Y1);
    for i:=1 to 2 do
      begin
        i1[i]:=2*h*i1[i];
        Y2[i]:=Y0[i]+(i1[i]+i2[i])
      end;
    T:=T2+hh1;
//    iterations
    l:=0;  // number of iterations
    repeat
      l:=l+1;
      find:=True;
      for i:=1 to 2 do
        Y[i]:=Y2[i]+hh1*Fdelta[i];
      i1:=Psi2Implicit(u2,T,Y);
      for i:=1 to 2 do
        begin
          i1[i]:=i1[i]-i1[i];
          i2[i]:=(h*h)*(h*h);
          i1[i]:=i2[i]*i1[i]/24
        end;
      i2:=FTY(u2,T2,Y2);
      i3:=FTY(u2,T1,Y1);
      for i:=1 to 2 do
        i2[i]:=i2[i]+4*i3[i];
      i3:=FTY(u2,T0,Y0);
      for i:=1 to 2 do
        begin
          i2[i]:=i2[i]+i3[i];
          i2[i]:=h*i2[i]/3;
          Y3[i]:=Y0[i]+(i2[i]+i1[i]);
          if (Y2[i].a=0) or (Y2[i].b=0)
            then if (Abs(Y3[i].a-Y2[i].a)>=1e-18)
                     or (Abs(Y3[i].b-Y2[i].b)>=1e-18)
                   then find:=False  // 1e-18 - accuracy
                   else
            else if (Abs((Y3[i].a-Y2[i].a)/Y2[i].a)>=1e-18)
                     or (Abs((Y3[i].b-Y2[i].b)/Y2[i].b)>=1e-18)
                   then find:=False
        end;
      Y2:=Y3
    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=500) or (n=1000) or (n=1500) or (n=2000)
             then Results (T2, Y2, 'i');
    T0:=T1;
    Y0:=Y1;
    T1:=T2;
    Y1:=Y2
  until (n=m) or not find;
  if z='f'
    then CloseFile (plik);
  Writeln;
  Writeln ('The end. Press Enter ...');
  Readln
end.