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

// The program solves the initial value problem y' = 0.5y, y(0) = 1 by
// the four-step explicit interval methods of Milne type.
// 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 intervalfunction = function (const T, Y : interval) : interval;

function Milne (const Tk1, Yk1, Tk2, Yk2,
                      Tk3, Yk3, Yk4,
                      DeltaT, DeltaY, H : interval;
                const F, PSI            : intervalfunction) : interval;
var H3H, H5, Tk1H3H, Yk1H3H, iv1, iv2 : interval;
begin
  iv1:=2*F(Tk1,Yk1)-F(Tk2,Yk2);
  iv1:=iv1+2*F(Tk3,Yk3);
  iv1:=4*H*iv1/3; // the term with Fk-1, Fk-2 and Fk-3
  H3H.a:=-3*H.b;
  H3H.b:=H.b;
  Tk1H3H:=Tk1+H3H; // the first argument for PSI
  Yk1H3H:=Yk1+H3H*F(DeltaT,DeltaY); // the second argument for PSI
  iv2:=PSI(Tk1H3H,Yk1H3H);
  iv2:=27*iv2-251*iv2;
  H5:=H*H;
  H5:=H5*H5;
  H5:=H*H5; // h^5
  iv2:=H5*iv2/720; // the term with PSI
  Result:=Yk4+iv1-iv2
end {Milne};

function F (const T, Y : interval) : interval;
begin
  Result:=Y/2
end {F};

function PSI (const T, Y : interval) : interval;
begin
  Result:=Y/32
end {PSI};

var i, k, m                                : Integer;
    hour, min, sec, msec                   : Word;
    exact_y, h_value, width                : Extended;
    Tk, Yk, Tk1, Yk1, Tk2, Yk2,
    Tk3, Yk3, Tk4, Yk4, H, DeltaT, DeltaY  : interval;
    hh                                     : AnsiString;
    file_name, left, right                 : string;
    agree, sf                              : Char;
    ftxt                                   : TextFile;
    time1, time2                           : TDateTime;

begin
  Write ('SOLVING THE INITIAL VALUE PROBLEM BY AN INTERVAL VERSION OF ');
  Writeln ('MILNE''S METHOD');
  Writeln;
  Writeln ('Example: y'' = 0.5y, y(0) = 1');
  Writeln;
  repeat
    Write ('Output (s - screen, f - file): ');
    Readln (sf)
  until (sf='s') or (sf='f');
  Writeln;
  if sf='s'
    then Writeln ('Output to the screen')
    else begin
           repeat
             repeat
               Write ('Enter file''s name (.txt will be added '
                      +'automatically): ');
               Readln (file_name);
               file_name:=file_name+'.txt';
               Writeln ('Your file''s name is: '+file_name);
               repeat
                 Write ('Do you agree (y - yes, n - no): ');
                 Readln (agree)
               until (agree='y') or (agree='n');
             until agree='y';
             AssignFile (ftxt, file_name);
             Rewrite (ftxt);
           until agree='y';
           Write (ftxt, 'SOLVING THE INITIAL VALUE PROBLEM BY AN INTERVAL ');
           Writeln (ftxt, 'VERSION OF MILNE''S METHOD');
           Writeln (ftxt, ' ');
           Writeln (ftxt, 'Example: y'' = 0.5y, y(0) = 1');
           Writeln (ftxt, ' ')
         end;
  Writeln;
  repeat
    Write ('Step size (>0 & <1):  h = ');
    Readln (hh);
    Val (hh, h_value, m)
  until (m=0) and (h_value>0) and (h_value<1);
  if sf='f'
    then Writeln (ftxt, 'Step size:            h = ', hh);
  repeat
    Write ('Number of steps (>3): m = ');
    Readln (m)
  until m>3;
  if sf='f'
    then Writeln (ftxt, 'Number of steps:      m = ', m);
  repeat
    Write ('Output every k steps: k = ');
    Readln (k)
  until k>0;
  Writeln;
  if sf='f'
    then begin
           Writeln (ftxt, ' ');
           Writeln ('Solving ...')
         end;
  H:=int_read(hh);
  Tk4:=0;     // t0 = tk-4
  Tk3:=Tk4+H; // t1 = tk-3
  Tk2:=Tk3+H; // t2 = tk-2
  Tk1:=Tk2+H; // t3 = tk-1
  DeltaT.a:=0;
  DeltaT.b:=1;
  DeltaY.a:=1;
  DeltaY.b:=right_read('1.65');
  Yk4:=1; // y0 = y(0)
  if hh='0.0001'
    then begin
           Yk3.a:=left_read('1.0000500012500208');
           Yk3.b:=right_read('1.0000500012500209'); // y1 = y(0.0001)
           Yk2.a:=left_read('1.0001000050001666');
           Yk2.b:=right_read('1.0001000050001667'); // y2 = y(0.0002)
           Yk1.a:=left_read('1.0001500112505625');
           Yk1.b:=right_read('1.0001500112505626')  // y3 = y(0.0003)
         end
    else if hh='0.0005'
           then begin
                  Yk3.a:=left_read('1.0002500312526043');
                  Yk3.b:=right_read('1.0002500312526044'); // y1 = y(0.0005)
                  Yk2.a:=left_read('1.0005001250208359');
                  Yk2.b:=right_read('1.0005001250208360'); // y2 = y(0.0010)
                  Yk1.a:=left_read('1.0007502813203256');
                  Yk1.b:=right_read('1.0007502813203257')  // y3 = y(0.0015)
                end
           else if hh='0.001'
                  then begin
                         Yk3.a:=left_read('1.0005001250208359');
                         Yk3.b:=right_read('1.0005001250208360');
                                                           // y1 = y(0.001)
                         Yk2.a:=left_read('1.0010005001667083');
                         Yk2.b:=right_read('1.0010005001667084');
                                                           // y2 = y(0.002)
                         Yk1.a:=left_read('1.0015011255627110');
                         Yk1.b:=right_read('1.0015011255627111')
                                                           // y3 = y(0.003)
                       end
                  else if hh='0.002'
                         then begin
                                Yk3.a:=left_read('1.0010005001667083');
                                Yk3.b:=right_read('1.0010005001667084');
                                                           // y1 = y(0.002)
                                Yk2.a:=left_read('1.0020020013340002');
                                Yk2.b:=right_read('1.0020020013340003');
                                                           // y2 = y(0.004)
                                Yk1.a:=left_read('1.0030045045033770');
                                Yk1.b:=right_read('1.0030045045033771')
                                                           // y3 = y(0.006)
                              end
                         else begin
                                Writeln ('Add aditional starting points in '
                                         +'source code - see lines 174-178.');
                                Readln;
                                Exit
// insert Yk3, Yk2 and Yk1 for another value of h
                                     end;
// writing initial data
  if sf='s'
    then begin
           Writeln ('Initial data:');
           Writeln ('-------------');
           Writeln;
           Write ('     T = [')
         end
    else begin
           Writeln (ftxt, 'Initial data:');
           Writeln (ftxt, '-------------');
           Writeln (ftxt, ' ');
           Write (ftxt, '     T = [')
         end;
  iends_to_strings (Tk4, left, right);
  if sf='s'
    then begin
           Writeln (left, ',', right, ']');
           Write ('     Y = [')
         end
    else begin
           Writeln (ftxt, left, ',', right, ']');
           Write (ftxt, '     Y = [')
         end;
  iends_to_strings (Yk4, left, right);
  if sf='s'
    then begin
           Writeln (left, ',', right, ']');
           Writeln;
           Write ('     T = [')
         end
    else begin
           Writeln (ftxt, left, ',', right, ']');
           Writeln (ftxt, ' ');
           Write (ftxt, '     T = [')
         end;
  iends_to_strings (Tk3, left, right);
  if sf='s'
    then begin
           Writeln (left, ',', right, ']');
           Write ('     Y = [')
         end
    else begin
           Writeln (ftxt, left, ',', right, ']');
           Write (ftxt, '     Y = [')
         end;
  iends_to_strings (Yk3, left, right);
  if sf='s'
    then begin
           Writeln (left, ',', right, ']');
           Writeln;
           Write ('     T = [')
         end
    else begin
           Writeln (ftxt, left, ',', right, ']');
           Writeln (ftxt, ' ');
           Write (ftxt, '     T = [')
         end;
  iends_to_strings (Tk2, left, right);
  if sf='s'
    then begin
           Writeln (left, ',', right, ']');
           Write ('     Y = [')
         end
    else begin
           Writeln (ftxt, left, ',', right, ']');
           Write (ftxt, '     Y = [')
         end;
  iends_to_strings (Yk2, left, right);
  if sf='s'
    then begin
           Writeln (left, ',', right, ']');
           Writeln;
           Write ('     T = [')
         end
    else begin
           Writeln (ftxt, left, ',', right, ']');
           Writeln (ftxt, ' ');
           Write (ftxt, '     T = [')
         end;
  iends_to_strings (Tk1, left, right);
  if sf='s'
    then begin
           Writeln (left, ',', right, ']');
           Write ('     Y = [')
         end
    else begin
           Writeln (ftxt, left, ',', right, ']');
           Write (ftxt, '     Y = [')
         end;
  iends_to_strings (Yk1, left, right);
  if sf='s'
    then begin
           Writeln (left, ',', right, ']');
           Writeln;
           Write ('DeltaT = [')
         end
    else begin
           Writeln (ftxt, left, ',', right, ']');
           Writeln (ftxt, ' ');
           Write (ftxt, 'DeltaT = [')
         end;
  iends_to_strings (DeltaT, left, right);
  if sf='s'
    then begin
           Writeln (left, ',', right, ']');
           Write ('DeltaY = [')
         end
    else begin
           Writeln (ftxt, left, ',', right, ']');
           Write (ftxt, 'DeltaY = [')
         end;
  iends_to_strings (DeltaY, left, right);
  if sf='s'
    then begin
           Writeln (left, ',', right, ']');
           Writeln;
           Writeln ('Solution:');
           Writeln ('---------');
           Writeln
         end
    else begin
           Writeln (ftxt, left, ',', right, ']');
           Writeln (ftxt, ' ');
           Writeln (ftxt, 'Solution:');
           Writeln (ftxt, '---------');
           Writeln (ftxt, ' ')
         end;
  time1:=Time;
  for i:=4 to m do
    begin
      Tk:=Tk1+H;
      Tk2:=Tk1-H;
      Tk3:=Tk2-H;
      Yk:=Milne(Tk1,Yk1,Tk2,Yk2,Tk3,Yk3,Yk4,DeltaT,DeltaY,H,F,PSI);
      if i mod k = 0
        then begin
               iends_to_strings (Tk, left, right);
               if sf='s'
                 then begin
                        Write ('     T = [');
                        Writeln (left, ',', right, ']');
                        Write ('     Y = [')
                      end
                 else begin
                        Write (ftxt, '     T = [');
                        Writeln (ftxt, left, ',', right, ']');
                        Write (ftxt, '     Y = [')
                      end;
               iends_to_strings (Yk, left, right);
               if sf='s'
                 then Writeln (left, ',', right, ']')
                 else Writeln (ftxt, left, ',', right, ']');
               width:=int_width(Yk);
               if sf='s'
                 then Writeln ('         width =', width:11)
                 else Writeln (ftxt, '         width =', width:11);
               exact_y:=Exp(0.5*h_value*i);
               if sf='s'
                 then Write (' Exact =  ', exact_y:25)
                 else Write (ftxt, ' Exact =  ', exact_y:25);
               if (Yk.a<=exact_y) and (exact_y<=Yk.b)
                 then if sf='s'
                        then Writeln ('  Exact solution within the interval Y')
                        else Writeln (ftxt, '  Exact solution within the '
                                      +'interval Y')
                 else if sf='s'
                        then Writeln ('  Exact solution outside the interval Y')
                        else Writeln (ftxt, '  Exact solution outside the '
                                      +'interval Y');
               if sf='s'
                 then Writeln
                 else Writeln (ftxt, ' ')
             end;
      if i<m
        then begin
               Tk4:=Tk3;
               Yk4:=Yk3;
               Tk3:=Tk2;
               Yk3:=Yk2;
               Tk2:=Tk1;
               Yk2:=Yk1;
               Tk1:=Tk;
               Yk1:=Yk
             end
    end;
  time2:=Time;
  DecodeTime (time2-time1, hour, min, sec, msec);
  if sf='s'
    then begin
           Writeln;
           Write ('Time of calculations: ', hour, ':');
           if min<10
             then Write ('0', min, ':')
             else Write (min, ':');
           if sec<10
             then Writeln ('0' , sec, '.', msec)
             else Writeln (sec, '.', msec)
         end
    else begin
           Writeln (ftxt, ' ');
           Write (ftxt, 'Time of calculations: ', hour, ':');
           if min<10
             then Write (ftxt, '0', min, ':')
             else Write (ftxt, min, ':');
           if sec<10
             then Writeln (ftxt, '0' , sec, '.', msec)
             else Writeln (ftxt, sec, '.', msec);
           CloseFile (ftxt);
         end;
  Writeln;
  Writeln ('End of calculations. Press Enter');
  Readln
end.