program ABVSS_Example_4_ASMD;
{------------------------------------------------------------------------------}
{                                                                              }
{ The program solves the initial value problem (the motion of a simple         }
{ pendulum                                                                     }
{                    y'[1] = -u^2*y[2],  y'[2] = y[1],                         }
{                     y[1](0) = 0,  y[2](0) = phi_0,                           }
{ where u = sqrt(g/L), g = 9.80665 (the gravitational acceleration at the      }
{ Earth surface), L = 1 (the pendulum length), phi_0 = pi/6 (an initial        }
{ angle), by an interval third order Adams-Bashforth method with finding a     }
{ step size for each integration step. The given problem has the solution      }
{               y[1] = -u*ph_0*sin(ut), y[2] = phi_0*cos(ut).                  }
{ 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_Example_4 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);
      i6  : interval = (a:6; b:6);
      im1 : interval = (a:-1; b:-1);

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

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 u, T : interval;
            const Y    : interval_array) : interval_array;
var u2 : interval;
begin
  u2:=imul(u,u);
  Result[1]:=imul(imul(im1,u2),Y[2]);
  Result[2]:=Y[1]
end {F};

function PSI3(const u, T : interval;
              const Y    : interval_array) : interval_array;
var u4 : interval;
begin
  u4:=imul(u,u);
  u4:=imul(u4,u4);
  Result[1]:=imul(u4,Y[1]);
  Result[2]:=imul(u4,Y[2])
end {PSI3};

procedure IntervalAB3 (const u, T1, T2, T3,
                             ih2, ih1, ih, delta_t : interval;
                       const Y1, Y2, Y3, delta_y   : interval_array;
                       out Y                       : interval_array);
var p     : Integer;
    PSI   : interval_array;
    Th, g : interval;
begin
  Th:=iadd(ih2,ih1);
  Th.a:=-Th.b;
  Th.b:=ih.b;
  for p:=1 to 2 do
    Y[p]:=iadd(Y1[p],imul(Th,F(u,delta_t,delta_y)[p]));
  PSI:=PSI3(u,iadd(T1,Th),Y);
  for p:=1 to 2 do
    begin
      Y[p]:=isub(F(u,T2,Y2)[p],F(u,T3,Y3)[p]);
      Y[p]:=imul(idiv(ih1,ih2),Y[p]);
      Y[p]:=isub(isub(F(u,T1,Y1)[p],F(u,T2,Y2)[p]),Y[p]);
      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[p]:=idiv(imul(g,Y[p]),i2);
      g:=isub(F(u,T1,Y1)[p],F(u,T2,Y2)[p]);
      g:=idiv(imul(idiv(ih,ih1),g),i2);
      Y[p]:=iadd(iadd(F(u,T1,Y1)[p],g),Y[p]);
      Y[p]:=iadd(Y1[p],imul(ih,Y[p]));
      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[p]:=iadd(Y[p],imul(imul(imul(ih,ih),imul(ih,ih)),imul(g,PSI[p])))
    end
end {IntervalAB3};

procedure Find_stepsize (const Y                  : array of interval_array;
                         const lambda, PSI_width,
                               Y_width, eps       : Extended;
                         const max_it             : Integer;
                         var h                    : array of Extended;
                         out it, status           : Integer);
{------------------------------------------------------------------------------}
{                                                                              }
{ The function Find_stepsize calculates the step size for the next integration }
{ step for an interval method of Adams_Bashforth type of fourth order used to  }
{ solve the initial value problem                                              }
{                    y'[1] = -u^2*y[2],  y'[2] = y[1],                         }
{                     y[1](0) = 0,  y[2](0) = phi_0,                           }
{ where u = sqrt(g/L), g = 9.80665 (the gravitational acceleration at the      }
{ Earth surface), L = 1 (the pendulum length), phi_0 = pi/6 (an initial        }
{ angle). The given problem has the solution                                   }
{               y[1] = -u*ph_0*sin(ut), y[2] = phi_0*cos(ut).                  }
{ To find this step the Newton iteration method is used.                       }
{ Data:                                                                        }
{   Y         - two-dimensional array containing interval enclosures of        }
{               solution at previous k steps (this array must be indexed from  }
{               zero to k-1 in the first dimension and from 1 to 2 in the      }
{               second one,                                                    }
{   lambda    - constant occuring in the estimation                            }
{                 width(F(T,Y)) <= lambda*(width(T)+width(Y)),                 }
{               where width(F(T,Y) = max(p=1,2) width(F(t,Y)[p]) and           }
{               width(Y) = max(p=1,2) width(Y[p]),                             }
{   PSI_width - max(p=1,2) width(PSI3(u, delta_t, delta_y)[p]),                }
{   Y_width   - required max(p=1,2) width(Y[p]),                               }
{   eps       - accuracy for finding h for the next inegration step,           }
{   max_it    - maximal number of iterations,                                  }
{   h         - array (indexed from zero) containing step sizes at previous    }
{               k-1 steps and an initial approximation of h for the next step  }
{               included in h[k].                                              }
{ Results:                                                                     }
{   h[k]   - step size for the next integration step,                          }
{   it     - number of iterations performed to find h[k],                      }
{   status - a variable which within the procedure is assigned the value of:   }
{              2 - if PSI_width <= 0 or any of h[l] for l = 1, 2, ... , k-1 is }
{                  less than zero or Y_width < 0 or Y_width > 0.1 or eps <= 0  }
{                  or eps > 0.1,                                               }
{              3 - if the required accuracy eps is not achieved in max_it      }
{                  iterations,                                                 }
{              4 - if the calulated h[k] was less or equal to zero,            }
{              0 - otherwise.                                                  }
{ Note: If status = 1, 2, 3 or 4, then h[k] = 0 and it = 0.                    }
{                                                                              }
{------------------------------------------------------------------------------}
var k, p                          : Integer;
    h0, h1, q, q_prime, Yk_widths : Extended;
begin
  status:=0;
  if PSI_width<=0
    then status:=2
    else if (Y_width<0) or (Y_width>0.1)
           then status:=2
           else if (eps<=0) or (eps>0.1)
                  then status:=2
                  else begin
                         k:=1;
                         repeat
                           if h[k]<0
                             then status:=2
                             else k:=k+1
                         until (k>3) or (status=2)
                       end;
 if status<>0
   then begin
          h[3]:=0;
          it:=0
        end
   else begin
          Yk_widths:=0;
          for k:=2 downto 0 do
            begin
              q:=int_width(Y[k,1]);
              if int_width(Y[k,2])>q
                then q:=int_width(Y[k,2]);
              Yk_widths:=Yk_widths+(k+1)*q
            end;
          it:=0;
          h0:=h[2];
          status:=3;
          repeat
            it:=it+1;
            if (h0/h[2]<=1) and (h0*(h0+h[2])/(h[2]*h[2])<=1)
              then begin
                     q:=h0;
                     q_prime:=1
                   end
              else if h0*(h0+h[2])/(h[2]*h[1])<=h0/h[2]
                     then begin
                            q:=h0*h0/h[2];
                            q_prime:=2*h0/h[2]
                          end
                     else begin
                            q:=h0*h0*(h0+h[2])/(h[2]*h[1]);
                            q_prime:=h0*(3*h0+2*h[2])/(h[2]*h[1])
                          end;
            h1:=int_width(Y[2,1]);
            if int_width(Y[2,2])>h1
              then h1:=int_width(Y[2,2]);
            h1:=h0-(h0*h0*((3*h0+4*(2*h[2]+h[1]))*h0
                +6*h[2]*(h[2]+h[1]))*PSI_width/72+lambda*q*Yk_widths+h1-Y_width)
                /(h0*(2*h0*h0+2*h0*(2*h[2]+h[1])+h[2]
                *(h[2]+h[1]))*PSI_width/12+lambda*q_prime*Yk_widths);
            if h1<=0
              then status:=4
              else if Abs(h1-h0)>eps
                     then h0:=h1
                     else status:=0
          until (it=max_it) or (status=0) or (status=4);
          if status=0
            then h[3]:=h1
            else begin
                   h[3]:=0;
                   it:=0
                 end
        end
end {Find_stepsize};

var it, k, max_it, p, status             : Integer;
    a, eps, lambda, PSI_width, real_h,
    real_t, real_u, w, Y_width           : Extended;
    file_name, left, right, time         : string;
    agree, output                        : Char;
    finish                               : Boolean;
    b_left, b_right, h                   : array of Extended;
    delta_t, ih, ih1, ih2, T1, T2, T3, u : interval;
    delta_y, Y1, Y2, Y3, Yh              : interval_array;
    Yk                                   : array of interval_array;
    results                              : TextFile;
    time1, time2                         : TDateTime;

begin
  Writeln ('Solving the initial value problem');
  Writeln ('y'' = -u^2*y[2], y''[2] = y[1],');
  Writeln ('y[1](0) = 0, y[2](0) = phi_0,');
  Writeln ('where u = sqrt(g/L), g = 9.80665, L = 1, phi_0 = pi/6');
  Writeln ('for t in [0, a] and y[l] in {b_left[l], b_right[l]} (l = 1,2)');
  Writeln ('by an interval third order Adams-Bashforth method with finding');
  Writeln ('a step size for each integration step.');
  Writeln;
  Writeln ('Enter data:');
  repeat
    Write ('                  a = ');
    Readln (a);
    if a<=0
      then Writeln ('a must be greater then 0')
  until a>0;
  SetLength (b_left, 3);
  SetLength (b_right, 3);
  for p:=1 to 2 do
    begin
      Write ('          b_left[', p, '] = ');
      Readln (b_left[p]);
      repeat
        Write ('         b_right[', p, '] = ');
        Readln (b_right[p]);
        if b_right[p]<=b_left[p]
          then Writeln ('b_right[', p, '] must be greater then b_left[', p, ']')
      until b_right[p]>b_left[p]
    end;
  Str (a:26, right);
  delta_t.a:=0;
  delta_t.b:=right_read(right);
  for p:=1 to 2 do
    begin
      Str (b_left[p]:26, left);
      delta_y[p].a:=left_read(left);
      Str (b_right[p]:26, right);
      delta_y[p].b:=right_read(right)
    end;
  SetLength (h, 4);
  SetLength (Yk, 3);
  h[1]:=0.0001;
  h[2]:=h[1];
  h[3]:=h[1];
  real_u:=Sqrt(9.80665);
  Str (real_u:26, left);
  u:=int_read(left);
  repeat
    Write ('  required Y width = ');
    Readln (Y_width);
    if (Y_width<=0) or (Y_width>0.1)
      then Writeln ('required Y width must be in the interval (0, 0.1]')
  until (Y_width>0) and (Y_width<=0.1);
  repeat
    Write ('accuracy for h, eps = ');
    Readln (eps);
    if (eps<=0) or (eps>0.1)
      then Writeln ('accuracy eps for finding h must be in the interval '
                    +'(0, 0.1]')
  until (eps>0) and (eps<=0.1);
  repeat
    Write ('             max_it = ');
    Readln (max_it);
    if max_it<=0
      then Writeln ('max number of iterations must be greater than 0')
  until max_it>0;
  lambda:=Sqr(real_u)*Pi/6;
  Writeln;
  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 (agree)
           until agree='y';
           Assign (results, file_name);
           Rewrite (results);
           Writeln (results, 'Solving the initial value problem');
           Writeln (results, 'y'' = -u^2*y[2], y''[2] = y[1],');
           Writeln (results, 'y[1](0) = 0, y[2](0) = phi_0,');
           Writeln (results, 'where u = sqrt(g/L), g = 9.80665, L = 1, phi_0'
                    +' = pi/6');
           Writeln (results, 'for t in [0, a] and y[l] in {b_left[l], '
                    +'b_right[l]}, (l = 1,2)');
           Writeln (results, 'by an interval third order Adams-Bashforth method'
                    +' with finding');
           Writeln (results, 'a step size for each integration step.');
           Writeln (results, ' ');
           Writeln (results, 'a          = ', a:4:2);
           for p:=1 to 2 do
             begin
               Writeln (results, 'b_left[', p, ']  = ', b_left[p]:4:2);
               Writeln (results, 'b_right[', p, '] = ', b_right[p]:4:2)
             end;
           for k:=1 to 3 do
             Writeln (results, 'h[', k:2, '] = ', h[k]:4:2);
           Writeln (results, '  required Y width = ', Y_width);
           Writeln (results, 'accuracy for h, eps = ', eps);
           Writeln (results, 'max_it              = ', max_it)
         end;
  Writeln;
  Writeln ('Press Enter ...');
  Readln;
  if output='s'
    then begin
           Write ('Enclosures of the exact solution obtained by the interval '
                  +'Adams-Bashforth method of third order');
           Writeln
         end
    else begin
           Write (results, 'Enclosures of the exact solution obtained by the '
                  +'interval Adams-Bashforth method of third order');
           Writeln (results, ' ')
         end;
  k:=3;
  real_t:=h[1]+h[2];
  Str (real_t:26, left);
  T1:=int_read(left);
  left:='-1.0269499194046190E-0003';
  right:='-1.0269499194046189E-0003';
  Y1[1].a:=left_read(left);
  Y1[1].b:=right_read(right);
  left:='5.2359867290330357E-0001';
  right:='5.2359867290330358E-0001';
  Y1[2].a:=left_read(left);
  Y1[2].b:=right_read(right);
  for p:=1 to 2 do
    Yk[2,p]:=Y1[p];
  w:=h[1];
  Str (w:26, left);
  T2:=int_read(left);
  left:='-5.1347498487965657E-0004';
  right:='-5.1347498487965656E-0004';
  Y2[1].a:=left_read(left);
  Y2[1].b:=right_read(right);
  left:='5.2359874992454941E-0001';
  right:='5.2359874992454942E-0001';
  Y2[2].a:=left_read(left);
  Y2[2].b:=right_read(right);
  for p:=1 to 2 do
    Yk[1,p]:=Y2[p];
  T3:=i0;
  Y3[1]:=i0;
  left:='5.2359877559829887E-0001';
  right:='5.2359877559829888E-0001';
  Y3[2].a:=left_read(left);
  Y3[2].b:=right_read(right);
  for p:=1 to 2 do
    Yk[0,p]:=Y3[p];
  iends_to_strings (T3, left, right);
  if output='s'
    then begin
           Writeln;
           Writeln ('Initial intervals:');
           Writeln ('  T[0] = [', left, ',', right, ']')
         end
    else begin
           Writeln (results, ' ');
           Writeln (results, 'Initial intervals:');
           Writeln (results, '  T[0] = [', left, ',', right, ']')
         end;
  for p:=1 to 2 do
    begin
      iends_to_strings (Y3[p], left, right);
      if output='s'
        then Writeln ('Y[', p, ',0] = [', left, ',', right, ']')
        else Writeln (results, 'Y[', p, ',0] = [', left, ',', right, ']')
    end;
  iends_to_strings (T2, left, right);
  if output='s'
    then Writeln ('  T[1] = [', left, ',', right, ']')
    else Writeln (results, '  T[1] = [', left, ',', right, ']');
  for p:=1 to 2 do
    begin
      iends_to_strings (Y2[p], left, right);
      if output='s'
        then Writeln ('Y[', p, ',1] = [', left, ',', right, ']')
        else Writeln (results, 'Y[', p, ',1] = [', left, ',', right, ']')
    end;
  iends_to_strings (T1, left, right);
  if output='s'
    then Writeln ('  T[2] = [', left, ',', right, ']')
    else Writeln (results, '  T[2] = [', left, ',', right, ']');
  for p:=1 to 2 do
    begin
      iends_to_strings (Y1[p], left, right);
      if output='s'
        then begin
               Writeln ('Y[', p, ',2] = [', left, ',', right, ']');
               if p=2
                 then Writeln
             end
        else begin
               Writeln (results, 'Y[', p, ',2] = [', left, ',', right, ']');
               if p=2
                 then Writeln (results, ' ')
             end
    end;
  Writeln;
  if output='f'
    then Writeln ('Calculating. Please, wait ...');
  time1:=Now;
  PSI_width:=int_width(PSI3(u,delta_t,delta_y)[1]);
  if int_width(PSI3(u,delta_t,delta_y)[2])>PSI_width
    then PSI_width:=int_width(PSI3(u,delta_t,delta_y)[2]);
  finish:=False;
  repeat
    Find_stepsize (Yk, lambda, PSI_width, Y_width, eps, max_it, h, it, status);
    if status=0
      then begin
             real_h:=h[3];
             if real_t+real_h<a
               then begin
                      real_t:=real_t+real_h;
                      if output='s'
                        then Writeln ('Step size for the next integration step '
                                      +'h =', real_h:26, ' found in ', it,
                                      ' iterations')
                        else Writeln (results, 'Step size for the next '
                                      +'integration step h =', real_h:26,
                                      'found in ', it, ' iterations');
                      Str (real_h:26, left);
                      ih:=int_read(left)
                    end
               else begin
                      real_h:=a-real_t;
                      real_t:=a;
                      if output='s'
                        then Writeln ('Step size for the last integration step '
                                      +'h =', real_h:26)
                        else Writeln (results, 'Step size for the last '
                                      +'integration step h =', real_h:26);
                      Str (real_h:26, left);
                      ih:=int_read(left);
                      finish:=True
                    end;
             Str (h[2]:26, left);
             ih1:=int_read(left);
             Str (h[1]:26, left);
             ih2:=int_read(left);
             IntervalAB3 (u, T1, T2, T3, ih2, ih1, ih, delta_t, Y1, Y2, Y3,
                         delta_y, Yh);
             if not (real_t=a)
               then begin
                      for p:=1 to 2 do
                        begin
                          Y3[p]:=Y2[p];
                          Yk[0,p]:=Y3[p]
                        end;
                      h[1]:=h[2];
                      T3:=T2;
                      for p:=1 to 2 do
                        begin
                          Y2[p]:=Y1[p];
                          Yk[1,p]:=Y2[p]
                        end;
                      h[2]:=h[3];
                      T2:=T1;
                      for p:=1 to 2 do
                        begin
                          Y1[p]:=Yh[p];
                          Yk[2,p]:=Y1[p]
                        end;
                      Str (real_t:26, left);
                      T1:=int_read(left)
                    end;
             for p:=1 to 2 do
               begin
                 iends_to_strings (Yh[p], left, right);
                 w:=int_width(Yh[p]);
                 if output='s'
                   then if p=1
                          then begin
                                 Writeln ('t = ', real_t:10:8, '  Y[', p, ',',
                                          k:3, '] = [', left, ', ', right,
                                          ']  width = ', w:11);
                                 Writeln ('                   exact =  ',
                                          -real_u*(Pi/6)*Sin(real_u*real_t):25)
                               end
                          else begin
                                 Writeln ('                Y[', p, ',', k:3,
                                          '] = [', left, ', ', right,
                                          ']  width = ', w:11);
                                 Writeln ('                   exact =  ',
                                          (Pi/6)*Cos(real_u*real_t):25)
                               end
                   else if p=1
                          then begin
                                 Writeln (results, 't = ', real_t:10:8, '  Y[',
                                          p, ',', k:3, '] = [', left, ', ',
                                          right, ']  width = ', w:11);
                                 Writeln (results, '                   exact ='
                                          +'  ',
                                          -real_u*(Pi/6)*Sin(real_u*real_t):25)
                               end
                          else begin
                                 Writeln (results, '                Y[', p, ',',
                                          k:3, '] = [', left, ', ', right,
                                          ']  width = ', w:11);
                                 Writeln (results, '                   exact ='
                                          +'  ', (Pi/6)*Cos(real_u*real_t):25)
                               end
               end;
             k:=k+1
           end
      else begin
             if output='s'
               then Writeln ('status = ', status, ' for k = ', k)
               else Writeln (results, 'status = ', status, ' for k = ', k);
             finish:=True
           end;
  until finish;
  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 (b_left);
  Finalize (b_right);
  Finalize (h);
  Finalize (Yk);
  if output='f'
    then CloseFile (results);
  Writeln;
  Writeln ('END OF PROGRAM');
  Readln
end.