program IMM_Example_2; {$APPTYPE CONSOLE} uses System.SysUtils, Winapi.Windows, IntervalArithmetic32and64; // The program solves the initial value problem y' = 0.5y, y(0) = 1 // by the four-step explicit interval methods of Adams-Bashforth type // or by the four-step explicit interval method of Nystrom 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 interval_function = function (const T, Y : interval) : interval; function FTY (const T, Y : interval) : interval; begin Result:=Y/2 end {FTY}; function PSI (const T, Y : interval) : interval; begin Result:=Y/32 end {PSI}; procedure AB4 (const T0, Y0, Y1, Y2, Y3, Fdelta : interval; const h : Extended; const FTY, PSI : interval_function; out T1, T2, T3, T4, Y4 : interval); // explicit interval four-step method of Adams-Bashforth type var i1, i2, i3, i4, i5, ih, hh, T, Y : interval; s : string; begin Str (h:26, s); ih.a:=left_read(s); ih.b:=right_read(s); Str (-3*h:26, s); hh.a:=left_read(s); hh.b:=ih.b; T1:=T0+ih; T2:=T1+ih; T3:=T2+ih; T4:=T3+ih; T:=T3+hh; Y:=Y3+hh*Fdelta; i1:=PSI(T,Y); i2:=ih*ih; i2:=i2*i2; i2:=i2*ih; i2:=271*i2*i1; i1:=i2/720; i2:=FTY(T3,Y3); i3:=FTY(T2,Y2); i4:=FTY(T1,Y1); i5:=FTY(T0,Y0); i2:=55*i2-59*i3; i2:=i2+37*i4; i2:=i2-9*i5; i2:=ih*i2/24; Y4:=Y3+(i2+i1) end {AB4}; procedure N4 (const T0, Y0, Y1, Y2, Y3, Fdelta : interval; const h : Extended; const FTY, PSI : interval_function; out T1, T2, T3, T4, Y4 : interval); // explicit interval four-step method of Nystrom type var i1, i2, i3, i4, i5, ih, hh, T, Y : interval; s : string; begin Str (h:26, s); ih.a:=left_read(s); ih.b:=right_read(s); Str (-3*h:26, s); hh.a:=left_read(s); hh.b:=ih.b; T1:=T0+ih; T2:=T1+ih; T3:=T2+ih; T4:=T3+ih; T:=T3+hh; Y:=Y3+hh*Fdelta; i1:=PSI(T,Y); i1:=251*i1-19*i1; i2:=ih*ih; i2:=i2*i2; i2:=i2*ih*i1; i1:=i2/720; i2:=FTY(T3,Y3); i3:=FTY(T2,Y2); i4:=FTY(T1,Y1); i5:=FTY(T0,Y0); i2:=8*i2-5*i3; i2:=i2+4*i4; i2:=i2-i5; i2:=ih*i2/3; Y4:=Y2+(i2+i1) end {N4}; var n, m, method, output_step : Integer; amh, real_h, w : Extended; delta_t, delta_y, Fdelta, T0, Y0, T1, Y1, T2, Y2, T3, Y3, T4, Y4 : Interval; left, right, strg : string; s : AnsiString; find : Boolean; z : Char; plik : TextFile; procedure Results (const T, Y : interval); var hreal : Extended; hrealstring : string; begin hreal:=(T.a+T.b)/2; Str (hreal:3:1, hrealstring); iends_to_strings (T, left, right); if z='f' then Writeln (plik, 'T('+hrealstring+') = ['+left+','+right+']') else Writeln ('T('+hrealstring+') = ['+left+','+right+']'); iends_to_strings (Y, left, right); if z='f' then Writeln (plik, 'Y('+hrealstring+') = ['+left+','+right+']') else Writeln ('Y('+hrealstring+') = ['+left+','+right+']'); w:=int_width(Y); if z='f' then Writeln (plik, 'width = ', w:11) else Writeln ('width = ', w:11) end {Results}; begin Writeln; Writeln ('Solving the initial value problem y'' = 0.5y, y(0) = 1'); Writeln; Writeln ('for delta_t = {t in R : 0 <= t <= 1}'); Writeln ('and delta_y = {y in R : 1 <= y <= 1.65}'); Writeln; Writeln ('by an interval multistep method.'); Writeln; Writeln ('Choose an interval method:'); Writeln; Writeln (' 1 - explicit four-step fourth order of Adams-Bashforth type'); Writeln (' 2 - explicit four-step fourth order of Nystrom type'); Writeln; repeat Write ('Your choice of method (1 or 2): '); Readln (method) until (method=1) or (method=2); Writeln; Writeln ('You have choosen'); Write ('the '); case method of 1 : Write ('explicit four-step fourth order method of Adams-Bashforth'); 2 : Write ('explicit four-step fourth order method of Nystrom') end; Writeln (' type.'); Writeln; Writeln ('Enter data for your problem:'); Writeln; w:=1; amh:=w; Str (w:26, s); delta_t.a:=0; delta_t.b:=left_read(s); delta_y.a:=left_read(s); w:=1.65; Str (w:26, s); delta_y.b:=right_read(s); Fdelta:=FTY(delta_t,delta_y); Write ('step size h = '); Readln (real_h); Write ('Y0 = [left_y, right_y] left_y = '); Readln (w); Str (w:26, s); Y0.a:=left_read(s); Write (' right_y = '); Readln (w); Str (w:26, s); Y0.b:=right_read(s); Write ('Y1 = [left_y, right_y] left_y = '); Readln (w); Str (w:26, s); Y1.a:=left_read(s); Write (' right_y = '); Readln (w); Str (w:26, s); Y1.b:=right_read(s); Write ('Y2 = [left_y, right_y] left_y = '); Readln (w); Str (w:26, s); Y2.a:=left_read(s); Write (' right_y = '); Readln (w); Str (w:26, s); Y2.b:=right_read(s); Write ('Y3 = [left_y, right_y] left_y = '); Readln (w); Str (w:26, s); Y3.a:=left_read(s); Write (' right_y = '); Readln (w); Str (w:26, s); Y3.b:=right_read(s); find:=False; repeat Write ('maximum nunmber of steps (m*h <= 1) m = '); Readln (m); if (m*real_h<=amh) and (m>0) then find:=True until find; repeat Write ('output every k steps (k <= m) output_step = '); Readln (output_step) until output_step<=m; Writeln; repeat Write ('Output (s - screen, f - file) : '); Readln (z) until (z='s') or (z='f'); if z='f' then begin case method of 1 : strg:='AB4'; 2 : strg:='N4' end; strg:='IMM2'+strg+'RESULTS.TXT'; AssignFile (plik, strg); Rewrite (plik); Writeln ('Output to '+strg+' in current directory') end else Writeln; case method of 1 : strg:='k = 4, Adams-Bashforth (explicit)'; 2 : strg:='k = 4, Nystrom (explicit)' end; if z='f' then begin Writeln (plik, ' '); Writeln (plik, strg) end else Writeln (strg); for n:=4 to m do begin if method=1 then AB4 (T0, Y0, Y1, Y2, Y3, Fdelta, real_h, FTY, PSI, T1, T2, T3, T4, Y4) else N4 (T0, Y0, Y1, Y2, Y3, Fdelta, real_h, FTY, PSI, T1, T2, T3, T4, Y4); if (n mod output_step)=0 then Results (T4, Y4); if n