program ABVSS_Example_1;
{------------------------------------------------------------------------------}
{                                                                              }
{ 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 one can use iadd, isub, imul and     }
{       idiv functions from our IntervalArithmetic32and64 unit (see            }
{       ABVSS_Example_1_ASMD program). In this program we use overloading      }
{       operators also implemented in our unit.                                }
{                                                                              }
{------------------------------------------------------------------------------}
{$APPTYPE CONSOLE}
{$R *.res}
uses System.SysUtils, System.DateUtils, IntervalArithmetic32and64;

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:=Y/2
end {F};
function PSI1(const T, Y : interval) : interval;
begin
  Result:=Y/4
end {PSI1};
function PSI2(const T, Y : interval) : interval;
begin
  Result:=Y/8
end {PSI2};
function PSI3(const T, Y : interval) : interval;
begin
  Result:=Y/16
end {PSI3};
function PSI4(const T, Y : interval) : interval;
begin
  Result:=Y/32
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:=0;
  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:=1;
  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:=Y1+Th*F(delta_t,delta_y);
    PSI:=PSI1(T1+Th,PSI);
    Y:=Y1+ih*F(T1,Y1)+(ih*ih/2)*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:=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:=0;
  Y2:=1;
  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:=Y1+Th*F(delta_t,delta_y);
    PSI:=PSI2(T1+Th,PSI);
    Y:=F(T1,Y1)-F(T2,Y2);
    Y:=F(T1,Y1)+(ih/ih1)*Y/2;
    Y:=Y1+ih*Y;
    g:=(1/3+(ih1/ih)/2)/2;
    Y:=Y+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;

    T2:=T1;
    Y2:=Y1;
    T1:=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:=0;
  Y3:=1;
  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:=ih2+ih1;
    Th.a:=-Th.b;
    Th.b:=ih.b;
    PSI:=Y1+Th*F(delta_t,delta_y);
    PSI:=PSI3(T1+Th,PSI);
    Y:=F(T2,Y2)-F(T3,Y3);
    Y:=(ih1/ih2)*(F(T2,Y2)-F(T3,Y3));
    Y:=F(T1,Y1)-F(T2,Y2)-Y;
    g:=(ih/ih1)*((ih+ih1)/(ih1+ih2));
    g:=(1-ih/(3*(ih+ih1)))*g;
    Y:=g*Y/2;
    g:=(ih/ih1)*(F(T1,Y1)-F(T2,Y2))/2;
    Y:=Y1+ih*(F(T1,Y1)+g+Y);
    g:=(ih1/ih)*((ih1+ih2)/ih)/2;
    g:=((2*ih1+ih2)/ih)/3+g;
    g:=(1/4+g)/6;
    Y:=Y+ih*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;
    T3:=T2;
    Y3:=Y2;
    T2:=T1;
    Y2:=Y1;
    T1:=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:=0;
  Y4:=1;
  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:=ih3+ih2+ih1;
    Th.a:=-Th.b;
    Th.b:=ih.b;
    PSI:=Y1+Th*F(delta_t,delta_y);
    PSI:=PSI4(T1+Th,PSI);
    Y:=(ih2/ih3)*(F(T3,Y3)-F(T4,Y4));
    g:=(ih1/ih2)*((ih1+ih2)/(ih2+ih3))*(F(T2,Y2)-F(T3,Y3)-Y);
    Y:=(ih1/ih2)*(F(T2,Y2)-F(T3,Y3));
    Y:=F(T1,Y1)-F(T2,Y2)-Y-g;
    g:=(ih+ih1+ih2)/(ih1+ih2+ih3);
    g:=((ih+ih1)/(ih1+ih2))*g;
    Y:=(ih/ih1)*g*Y;
    g:=ih/(ih+ih1+ih2);
    g:=((1-ih/(2*(ih+ih1)))/6)*g;
    Y:=((1-ih/(3*(ih+ih1)))/2-g)*Y;
    g:=(ih1/ih2)*(F(T2,Y2)-F(T3,Y3));
    g:=F(T1,Y1)-F(T2,Y2)-g;
    g:=(ih/ih1)*((ih+ih1)/(ih1+ih2))*g;
    Y:=(1-ih/(3*(ih+ih1)))*g/2+Y;
    g:=(ih/ih1)*(F(T1,Y1)-F(T2,Y2))/2;
    Y:=F(T1,Y1)+g+Y;
    Y:=Y1+ih*Y;
    g:=(ih1+ih2+ih3)/ih;
    g:=((ih1+ih2)/ih)*g;
    g:=((ih1/ih)*g)/2;
    g:=((ih2/ih)*((ih2+ih3)/ih))/3+g;
    g:=((ih1/ih)*(3*ih1+4*ih2+2*ih3)/ih)/3+g;
    g:=(3*ih1+2*ih2+ih3)/(4*ih)+g;
    g:=(1/5+g)/24;
    Y:=Y+ih*ih*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.