program ABVSS_Example_1_ASMD;
{------------------------------------------------------------------------------}
{                                                                              }
{ The program solves the initial value problem                                 }
{                           y' = 0.5y,  y(0) = 1,                              }
{ with the exact solution y = exp(0.5t), by interval methods of                }
{ Adams-Bashforth type with variable step sizes.                               }
{ Note: In order to speed up calculations in this program we use iadd, isub,   }
{       imul and idiv functions from our IntervalArithmetic32and64 unit, but   }
{       one can use overloading operators also implemented in this unit        }
{       (see ABVSS_Example1 program).                                          }
{                                                                              }
{------------------------------------------------------------------------------}
{$APPTYPE CONSOLE}
{$R *.res}
uses System.SysUtils, System.DateUtils, IntervalArithmetic32and64;

const i0  : interval = (a:0; b:0);
      i1  : interval = (a:1; b:1);
      i2  : interval = (a:2; b:2);
      i3  : interval = (a:3; b:3);
      i4  : interval = (a:4; b:4);
      i5  : interval = (a:5; b:5);
      i6  : interval = (a:6; b:6);
      i8  : interval = (a:8; b:8);
      i16 : interval = (a:16; b:16);
      i24 : interval = (a:24; b:24);
      i32 : interval = (a:32; b:32);

function modify_time (const time1, time2 : TDateTime) : string;
var code, days, hour_numbers : Integer;
    hours, time, tstr        : string;
    h, m, s, ms              : Word;
begin
  DecodeTime (time1-time2, h, m, s, ms);
  Str (h, tstr);
  if h<10
    then time:='0'+tstr
    else time:=tstr;
  Str (m, tstr);
  if m<10
    then time:=time+':0'+tstr
    else time:=time+':'+tstr;
  Str (s, tstr);
  if s<10
    then time:=time+':0'+tstr
    else time:=time+':'+tstr;
  Str (ms, tstr);
  if ms<10
    then time:=time+'.00'+tstr
    else if ms<100
           then time:=time+'.0'+tstr
           else time:=time+'.'+tstr;
  days:=DaysBetween(time1,time2);
  if days>0
    then begin
           hours:=Copy(time,1,2);
           Val (hours, hour_numbers, code);
           hour_numbers:=24*days+hour_numbers;
           Str (hour_numbers, hours);
           Result:=hours+Copy(time,3,10)
         end
    else Result:=time
end {modify_time};

function F (const T, Y : interval) : interval;
begin
  Result:=idiv(Y,i2)
end {F};
function PSI1(const T, Y : interval) : interval;
begin
  Result:=idiv(Y,i4)
end {PSI1};
function PSI2(const T, Y : interval) : interval;
begin
  Result:=idiv(Y,i8)
end {PSI2};
function PSI3(const T, Y : interval) : interval;
begin
  Result:=idiv(Y,i16)
end {PSI3};
function PSI4(const T, Y : interval) : interval;
begin
  Result:=idiv(Y,i32)
end {PSI4};

var k, m                          : Integer;
    a, b_left, b_right, real_t, w : Extended;
    file_name, left, right, time  : string;
    answer, output                : Char;
    h                             : array of Extended;
    delta_t, delta_y, g, ih, ih1,
    ih2, ih3, PSI, T1, T2, T3, T4,
    Th, Y, Y1, Y2, Y3, Y4         : interval;
    results                       : TextFile;
    time1, time2                  : TDateTime;

begin
  Writeln ('Solving the initial value problem');
  Writeln ('y'' = 0.5y, y(0) = 1');
  Writeln ('by interval Adams-Bashforth methods of order p (p = 1, 2, 3, 4)');
  Writeln ('with variable step sizes for t in [0, a] and y in [b_left,'
           +' b_right]');
  Writeln;
  Writeln ('Enter data:');
  repeat
    Write ('a       = ');
    Readln (a);
    if a<=0
      then Writeln ('a must be greater then 0')
  until a>0;
  Write ('b_left  = ');
  Readln (b_left);
  repeat
    Write ('b_right = ');
    Readln (b_right);
    if b_right<=b_left
      then Writeln ('b_right must be greater then b_left')
  until b_right>b_left;
  Str (a:26, right);
  delta_t.a:=0;
  delta_t.b:=right_read(right);
  Str (b_left:26, left);
  delta_y.a:=left_read(left);
  Str (b_right:26, right);
  delta_y.b:=right_read(right);
  Write ('Number of step sizes m = ');
  Readln (m);
  SetLength (h, m+1);
  for k:=1 to m do
    repeat
      Write ('h[', k:2, '] = ');
      Readln (h[k]);
      if h[k]<=0
        then Writeln ('h[', k:2, '] must be greater then zero')
    until h[k]>0;
  repeat
    Write ('Select the output (s - screen, f - file)                 output'
           +' = ');
    Readln (output)
  until (output='s') or (output='f');
  if output='f'
    then begin
           repeat
             Writeln ('Enter the name of file');
             Write ('(.txt will be added automatically)     file_name = ');
             Readln (file_name);
             file_name:=file_name+'.txt';
             Writeln ('The program will create a resulting file:          ',
                      file_name);
             Write ('Do you accept this name (y - yes, n - no) answer = ');
             Readln (answer)
           until answer='y';
           Assign (results, file_name);
           Rewrite (results);
           Writeln (results, 'Solving the initial value problem');
           Writeln (results, 'y'' = 0.5y, y(0) = 1');
           Writeln (results, 'by interval Adams-Bashforth methods of order p '
                    +'(p = 1, 2, 3, 4)');
           Writeln (results, 'with variable step sizes for t in [0, a] and y '
                    +'in [b_left, b_right]');
           Writeln (results, ' ');
           Writeln (results, 'a       = ', a:4:2);
           Writeln (results, 'b_left  = ', b_left:4:2);
           Writeln (results, 'b_right = ', b_right:4:2);
           Writeln (results, 'Number of step sizes m = ', m);
           for k:=1 to m do
             Writeln (results, 'h[', k:2, '] = ', h[k]:4:2);
           Writeln (results, ' ')
         end;
  Writeln;
  Writeln ('Press Enter ...');
  Readln;
  if output='s'
    then begin
           Writeln ('Enclosures of the exact solution obtained by the'
                    +' Adams-Bashforth method of first order');
           Writeln
         end
    else begin
           Writeln (results, 'Enclosures of the exact solution obtained by the'
                    +' Adams-Bashforth method of first order');
           Writeln (results, ' ')
         end;
  k:=0;
  real_t:=0;
  T1:=i0;
  iends_to_strings (T1, left, right);
  if output='s'
    then begin
           Writeln ('Initial intervals:');
           Writeln ('T[0] = [', left, ',', right, ']')
         end
    else begin
           Writeln (results, 'Initial intervals:');
           Writeln (results, 'T[0] = [', left, ',', right, ']')
         end;
  Y1:=i1;
  iends_to_strings (Y1, left, right);
  if output='s'
    then begin
           Writeln ('Y[0] = [', left, ',', right, ']');
           Writeln
         end
    else begin
           Writeln (results, 'Y[0] = [', left, ',', right, ']');
           Writeln (results, ' ')
         end;
  time1:=Now;
  repeat
    k:=k+1;
    real_t:=real_t+h[k];
    Str (h[k]:26, left);
    ih:=int_read(left);
    Th.a:=0;
    Th.b:=ih.b;
    PSI:=iadd(Y1,imul(Th,F(delta_t,delta_y)));
    PSI:=PSI1(iadd(T1,Th),PSI);
    Y:=iadd(Y1,imul(ih,F(T1,Y1)));
    Y:=iadd(Y,imul(idiv(imul(ih,ih),i2),PSI));
    iends_to_strings (Y, left, right);
    w:=int_width(Y);
    if output='s'
      then begin
             Writeln ('t = ', real_t:5:2, '  Y[', k:2, '] = [', left, ', ',
                      right, ']  width = ', w:11);
             Writeln ('           exact =  ', Exp(0.5*real_t):25)
           end
      else begin
             Writeln (results, 't = ', real_t:5:2, '  Y[', k:2, '] = [', left,
                      ', ', right, ']  width = ', w:11);
             Writeln (results, '           exact =  ', Exp(0.5*real_t):25);
           end;
    T1:=iadd(T1,ih);
    Y1:=Y
  until k=m;
  time2:=Now;
  time:=modify_time(time1,time2);
  if output='s'
    then begin
           Writeln;
           Writeln ('Time of calculations: ', time)
         end
    else begin
           Writeln (results, ' ');
           Writeln (results, 'Time of calculations: ', time)
         end;
  if output='s'
    then begin
           Writeln;
           Writeln ('Press Enter ...');
           Readln;
           Writeln ('Enclosures of the exact solution obtained by the'
                    +' Adams-Bashforth method of second order');
           Writeln;
         end
    else begin
           Writeln (results, ' ');
           Writeln (results, 'Enclosures of the exact solution obtained by the'
                    +' Adams-Bashforth method of second order');
           Writeln (results, ' ')
         end;
  k:=1;
  real_t:=h[1];
  Str (real_t:26, left);
  T1:=int_read(left);
  w:=Exp(0.5*real_t);
  Str (w:26, left);
  Y1:=int_read(left);
  T2:=i0;
  Y2:=i1;
  iends_to_strings (T2, left, right);
  if output='s'
    then begin
           Writeln ('Initial intervals:');
           Writeln ('T[0] = [', left, ',', right, ']')
         end
    else begin
           Writeln (results, 'Initial intervals:');
           Writeln (results, 'T[0] = [', left, ',', right, ']')
         end;
  iends_to_strings (Y2, left, right);
  if output='s'
    then Writeln ('Y[0] = [', left, ',', right, ']')
    else Writeln (results, 'Y[0] = [', left, ',', right, ']');
  iends_to_strings (T1, left, right);
  if output='s'
    then Writeln ('T[1] = [', left, ',', right, ']')
    else Writeln (results, 'T[1] = [', left, ',', right, ']');
  iends_to_strings (Y1, left, right);
  if output='s'
    then begin
           Writeln ('Y[1] = [', left, ',', right, ']');
           Writeln
         end
    else begin
           Writeln (results, 'Y[1] = [', left, ',', right, ']');
           Writeln (results, ' ')
         end;
  time1:=Now;
  repeat
    k:=k+1;
    real_t:=real_t+h[k];
    Str (h[k-1]:26, left);
    ih1:=int_read(left);
    Str (h[k]:26, left);
    ih:=int_read(left);
    Th.a:=-ih1.b;
    Th.b:=ih.b;
    PSI:=iadd(Y1,imul(Th,F(delta_t,delta_y)));
    PSI:=PSI2(iadd(T1,Th),PSI);
    Y:=isub(F(T1,Y1),F(T2,Y2));
    Y:=idiv(imul(idiv(ih,ih1),Y),i2);
    Y:=iadd(F(T1,Y1),Y);
    Y:=iadd(Y1,imul(ih,Y));
    g:=idiv(idiv(ih1,ih),i2);
    g:=idiv(iadd(idiv(i1,i3),g),i2);
    g:=imul(imul(imul(ih,ih),ih),g);
    Y:=iadd(Y,imul(g,PSI));
    iends_to_strings (Y, left, right);
    w:=int_width(Y);
    if output='s'
      then begin
             Writeln ('t = ', real_t:5:2, '  Y[', k:2, '] = [', left, ', ',
                      right, ']  width = ', w:11);
             Writeln ('           exact =  ', Exp(0.5*real_t):25)
           end
      else begin
             Writeln (results, 't = ', real_t:5:2, '  Y[', k:2, '] = [', left,
                      ', ', right, ']  width = ', w:11);
             Writeln (results, '           exact =  ', Exp(0.5*real_t):25)
           end;

    T2:=T1;
    Y2:=Y1;
    T1:=iadd(T1,ih);
    Y1:=Y
  until k=m;
  time2:=Now;
  time:=modify_time(time1,time2);
  if output='s'
    then begin
           Writeln;
           Writeln ('Time of calculations: ', time)
         end
    else begin
           Writeln (results, ' ');
           Writeln (results, 'Time of calculations: ', time)
         end;
  if output='s'
    then begin
           Writeln;
           Writeln ('Press Enter ...');
           Readln;
           Writeln ('Enclosures of the exact solution obtained by the'
                    +' Adams-Bashforth method of third order');
           Writeln
         end
    else begin
           Writeln (results, ' ');
           Writeln (results, 'Enclosures of the exact solution obtained by the'
                    +' Adams-Bashforth method of third order');
           Writeln (results, ' ')
         end;
  k:=2;
  real_t:=h[1]+h[2];
  Str (real_t:26, left);
  T1:=int_read(left);
  w:=Exp(0.5*real_t);
  Str (w:26, left);
  Y1:=int_read(left);
  w:=h[1];
  Str (w:26, left);
  T2:=int_read(left);
  w:=Exp(0.5*w);
  Str (w:26, left);
  Y2:=int_read(left);
  T3:=i0;
  Y3:=i1;
  iends_to_strings (T3, left, right);
  if output='s'
    then begin
           Writeln ('Initial intervals:');
           Writeln ('T[0] = [', left, ',', right, ']')
         end
    else begin
           Writeln (results, 'Initial intervals:');
           Writeln (results, 'T[0] = [', left, ',', right, ']')
         end;
  iends_to_strings (Y3, left, right);
  if output='s'
    then Writeln ('Y[0] = [', left, ',', right, ']')
    else Writeln (results, 'Y[0] = [', left, ',', right, ']');
  iends_to_strings (T2, left, right);
  if output='s'
     then Writeln ('T[1] = [', left, ',', right, ']')
     else Writeln (results, 'T[1] = [', left, ',', right, ']');
  iends_to_strings (Y2, left, right);
  if output='s'
    then Writeln ('Y[1] = [', left, ',', right, ']')
    else Writeln (results, 'Y[1] = [', left, ',', right, ']');
  iends_to_strings (T1, left, right);
  if output='s'
    then Writeln ('T[2] = [', left, ',', right, ']')
    else Writeln (results, 'T[2] = [', left, ',', right, ']');
  iends_to_strings (Y1, left, right);
  if output='s'
    then begin
           Writeln ('Y[2] = [', left, ',', right, ']');
           Writeln
         end
    else begin
           Writeln (results, 'Y[2] = [', left, ',', right, ']');
           Writeln (results, ' ')
         end;
  time1:=Now;
  repeat
    k:=k+1;
    real_t:=real_t+h[k];
    Str (h[k-2]:26, left);
    ih2:=int_read(left);
    Str (h[k-1]:26, left);
    ih1:=int_read(left);
    Str (h[k]:26, left);
    ih:=int_read(left);
    Th:=iadd(ih2,ih1);
    Th.a:=-Th.b;
    Th.b:=ih.b;
    PSI:=iadd(Y1,imul(Th,F(delta_t,delta_y)));
    PSI:=PSI3(iadd(T1,Th),PSI);
    Y:=isub(F(T2,Y2),F(T3,Y3));
    Y:=imul(idiv(ih1,ih2),Y);
    Y:=isub(isub(F(T1,Y1),F(T2,Y2)),Y);
    g:=imul(idiv(ih,ih1),idiv(iadd(ih,ih1),iadd(ih1,ih2)));
    g:=imul(isub(i1,idiv(ih,imul(i3,iadd(ih,ih1)))),g);
    Y:=idiv(imul(g,Y),i2);
    g:=isub(F(T1,Y1),F(T2,Y2));
    g:=idiv(imul(idiv(ih,ih1),g),i2);
    Y:=iadd(iadd(F(T1,Y1),g),Y);
    Y:=iadd(Y1,imul(ih,Y));
    g:=idiv(imul(idiv(ih1,ih),idiv(iadd(ih1,ih2),ih)),i2);
    g:=iadd(idiv(idiv(iadd(imul(i2,ih1),ih2),ih),i3),g);
    g:=idiv(iadd(idiv(i1,i4),g),i6);
    Y:=iadd(Y,imul(imul(imul(ih,ih),imul(ih,ih)),imul(g,PSI)));
    iends_to_strings (Y, left, right);
    w:=int_width(Y);
    if output='s'
      then begin
             Writeln ('t = ', real_t:5:2, '  Y[', k:2, '] = [', left, ', ',
                      right, ']  width = ', w:11);
             Writeln ('           exact =  ', Exp(0.5*real_t):25)
           end
      else begin
             Writeln (results, 't = ', real_t:5:2, '  Y[', k:2, '] = [', left,
                      ', ', right, ']  width = ', w:11);
             Writeln (results, '           exact =  ', Exp(0.5*real_t):25)
           end;
    T3:=T2;
    Y3:=Y2;
    T2:=T1;
    Y2:=Y1;
    T1:=iadd(T1,ih);
    Y1:=Y
  until k=m;
  time2:=Now;
  time:=modify_time(time1,time2);
  if output='s'
    then begin
           Writeln;
           Writeln ('Time of calculations: ', time)
         end
    else begin
           Writeln (results, ' ');
           Writeln (results, 'Time of calculations: ', time)
         end;
  if output='s'
    then begin
           Writeln;
           Writeln ('Press Enter ...');
           Readln;
           Writeln ('Enclosures of the exact solution obtained by the'
                    +' Adams-Bashforth method of fourth order');
           Writeln
         end
    else begin
           Writeln (results, ' ');
           Writeln (results, 'Enclosures of the exact solution obtained by the'
                    +' Adams-Bashforth method of fourth order');
           Writeln (results, ' ')
         end;
  k:=3;
  real_t:=h[1]+h[2]+h[3];
  Str (real_t:26, left);
  T1:=int_read(left);
  w:=Exp(0.5*real_t);
  Str (w:26, left);
  Y1:=int_read(left);
  w:=h[1]+h[2];
  Str (w:26, left);
  T2:=int_read(left);
  w:=Exp(0.5*w);
  Str (w:26, left);
  Y2:=int_read(left);
  w:=h[1];
  Str (w:26, left);
  T3:=int_read(left);
  w:=Exp(0.5*w);
  Str (w:26,left);
  Y3:=int_read(left);
  T4:=i0;
  Y4:=i1;
  iends_to_strings (T4, left, right);
  if output='s'
    then begin
           Writeln ('Initial intervals:');
           Writeln ('T[0] = [', left, ',', right, ']')
         end
    else begin
           Writeln (results, 'Initial intervals:');
           Writeln (results, 'T[0] = [', left, ',', right, ']')
         end;
  iends_to_strings (Y4, left, right);
  if output='s'
    then Writeln ('Y[0] = [', left, ',', right, ']')
    else Writeln (results, 'Y[0] = [', left, ',', right, ']');
  iends_to_strings (T3, left, right);
  if output='s'
    then Writeln ('T[1] = [', left, ',', right, ']')
    else Writeln (results, 'T[1] = [', left, ',', right, ']');
  iends_to_strings (Y3, left, right);
  if output='s'
    then Writeln ('Y[1] = [', left, ',', right, ']')
    else Writeln (results, 'Y[1] = [', left, ',', right, ']');
  iends_to_strings (T2, left, right);
  if output='s'
    then Writeln ('T[2] = [', left, ',', right, ']')
    else Writeln (results, 'T[2] = [', left, ',', right, ']');
  iends_to_strings (Y2, left, right);
  if output='s'
    then Writeln ('Y[2] = [', left, ',', right, ']')
    else Writeln (results, 'Y[2] = [', left, ',', right, ']');
  iends_to_strings (T1, left, right);
  if output='s'
    then Writeln ('T[3] = [', left, ',', right, ']')
    else Writeln (results, 'T[3] = [', left, ',', right, ']');
  iends_to_strings (Y1, left, right);
  if output='s'
    then begin
           Writeln ('Y[3] = [', left, ',', right, ']');
           Writeln
         end
    else begin
           Writeln (results, 'Y[3] = [', left, ',', right, ']');
           Writeln (results, ' ')
         end;
  time1:=Now;
  repeat
    k:=k+1;
    real_t:=real_t+h[k];
    Str (h[k-3]:26, left);
    ih3:=int_read(left);
    Str (h[k-2]:26, left);
    ih2:=int_read(left);
    Str (h[k-1]:26, left);
    ih1:=int_read(left);
    Str (h[k]:26, left);
    ih:=int_read(left);
    Th:=iadd(ih3,iadd(ih2,ih1));
    Th.a:=-Th.b;
    Th.b:=ih.b;
    PSI:=iadd(Y1,imul(Th,F(delta_t,delta_y)));
    PSI:=PSI4(iadd(T1,Th),PSI);
    Y:=isub(F(T3,Y3),F(T4,Y4));
    Y:=imul(idiv(ih2,ih3),Y);
    Y:=isub(isub(F(T2,Y2),F(T3,Y3)),Y);
    g:=imul(idiv(ih1,ih2),idiv(iadd(ih1,ih2),iadd(ih2,ih3)));
    g:=imul(g,Y);
    Y:=isub(F(T2,Y2),F(T3,Y3));
    Y:=imul(idiv(ih1,ih2),Y);
    Y:=isub(isub(isub(F(T1,Y1),F(T2,Y2)),Y),g);
    g:=idiv(iadd(iadd(ih,ih1),ih2),iadd(iadd(ih1,ih2),ih3));
    g:=imul(idiv(iadd(ih,ih1),iadd(ih1,ih2)),g);
    Y:=imul(imul(idiv(ih,ih1),g),Y);
    g:=idiv(ih,iadd(iadd(ih,ih1),ih2));
    g:=imul(idiv(isub(i1,idiv(ih,imul(i2,iadd(ih,ih1)))),i6),g);
    g:=isub(idiv(isub(i1,idiv(ih,imul(i3,iadd(ih,ih1)))),i2),g);
    Y:=imul(g,Y);
    g:=imul(idiv(ih1,ih2),isub(F(T2,Y2),F(T3,Y3)));
    g:=isub(isub(F(T1,Y1),F(T2,Y2)),g);
    g:=imul(imul(idiv(ih,ih1),idiv(iadd(ih,ih1),iadd(ih1,ih2))),g);
    g:=idiv(imul(isub(i1,idiv(ih,imul(i3,iadd(ih,ih1)))),g),i2);
    Y:=iadd(g,Y);
    g:=isub(F(T1,Y1),F(T2,Y2));
    g:=idiv(imul(idiv(ih,ih1),g),i2);
    Y:=iadd(iadd(F(T1,Y1),g),Y);
    Y:=iadd(Y1,imul(ih,Y));
    g:=idiv(iadd(iadd(ih1,ih2),ih3),ih);
    g:=imul(idiv(iadd(ih1,ih2),ih),g);
    g:=idiv(imul(idiv(ih1,ih),g),i2);
    g:=iadd(idiv(imul(idiv(ih2,ih),idiv(iadd(ih2,ih3),ih)),i3),g);
    g:=iadd(idiv(imul(idiv(ih1,ih),idiv(iadd(iadd(imul(i3,ih1),imul(i4,ih2)),
       imul(i2,ih3)),ih)),i3),g);
    g:=iadd(idiv(iadd(iadd(imul(i3,ih1),imul(i2,ih2)),ih3),imul(i4,ih)),g);
    g:=idiv(iadd(idiv(i1,i5),g),i24);
    Y:=iadd(Y,imul(imul(imul(imul(imul(ih,ih),imul(ih,ih)),ih),g),PSI));
    iends_to_strings (Y, left, right);
    w:=int_width(Y);
    if output='s'
      then begin
             Writeln ('t = ', real_t:5:2, '  Y[', k:2, '] = [', left, ', ',
                      right, ']  width = ', w:11);
             Writeln ('           exact =  ', Exp(0.5*real_t):25)
           end
      else begin
             Writeln (results, 't = ', real_t:5:2, '  Y[', k:2, '] = [', left,
                      ', ', right, ']  width = ', w:11);
             Writeln (results, '           exact =  ', Exp(0.5*real_t):25)
           end;
    T4:=T3;
    Y4:=Y3;
    T3:=T2;
    Y3:=Y2;
    T2:=T1;
    Y2:=Y1;
    T1:=iadd(T1,ih);
    Y1:=Y
  until k=m;
  time2:=Now;
  time:=modify_time(time1,time2);
  if output='s'
    then begin
           Writeln;
           Writeln ('Time of calculations: ', time)
         end
    else begin
           Writeln (results, ' ');
           Writeln (results, 'Time of calculations: ', time)
         end;
  Finalize (h);
  if output='f'
    then CloseFile (results);
  Writeln;
  Writeln ('END OF PROGRAM');
  Readln
end.