program ABVSS_Example_3_ASMD; {------------------------------------------------------------------------------} { } { The program solves the initial value problem } { y' = (y-t)/(y+t), y(0) = 4, } { with unknownc solution, by an interval fourth order method of } { Adams-Bashforth type with finding a step size for each integration step. } { 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_3 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); i10 : interval = (a:10; b:10); i13 : interval = (a:13; b:13); i16 : interval = (a:16; b:16); i24 : interval = (a:24; b:24); i40 : interval = (a:40; b:40); 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(isub(Y,T),iadd(Y,T)) end {F}; function PSI4(const T, Y : interval) : interval; var t2, x, y2 : interval; begin y2:=imul(Y,Y); t2:=imul(T,T); x:=isub(imul(imul(i16,y2),Y),imul(imul(i13,y2),T)); x:=iadd(x,imul(imul(i10,Y),T)); x:=isub(x,imul(imul(i3,t2),T)); x:=imul(i40,imul(iadd(y2,t2),x)); y2:=iadd(Y,T); t2:=imul(y2,y2); t2:=imul(t2,t2); t2:=imul(t2,t2); Result:=idiv(x,imul(t2,y2)) end {PSI4}; procedure IntervalAB4 (const T1, Y1, T2, Y2, T3, Y3, T4, Y4, ih3, ih2, ih1, ih, delta_t, delta_y : interval; out Y : interval); var PSI, Th, g : interval; begin 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)) end {IntervalAB4}; procedure Find_stepsize (const Y : array of interval; 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' = (y-t)/(y+t), y(0) = 4, } { with unknownc solution. To find this step the Newton iteration method is } { used. } { Data: } { Y - array containing interval enclosures of solution at previous } { k steps (this array must be indexed from zero to k-1), } { lambda - constant occuring in the estimation } { width(F(T,Y)) <= lambda*(width(T)+width(Y)), } { PSI_width - width of PSI(delta_t, delta_y), } { Y_width - required width of Y at t_k, } { 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 : 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>4) or (status=2) end; if status<>0 then begin h[4]:=0; it:=0 end else begin Yk_widths:=0; for k:=3 downto 0 do Yk_widths:=Yk_widths+(k+1)*int_width(Y[k]); it:=0; h0:=h[4]; status:=3; repeat it:=it+1; if (h0/h[3]<=1) and (h0*(h0+h[3])/(h[3]*h[2])<=1) and (h0*(h0+h[3])*(h0+h[3]+h[2])/(h[3]*h[2]*(h[2]+h[1]))<=1) and (h0*(h0+h[3])*(h0+h[3]+h[2])/((h[3]+h[2])*h[2]*h[1])<=1) then begin q:=h0; q_prime:=1 end else if (h0*(h0+h[3])/(h[3]*h[2])<=h0/h[3]) and (h0*(h0+h[3])*(h0+h[3]+h[2]) /(h[3]*h[2]*(h[2]+h[1]))<=h0/h[3]) and (h0*(h0+h[3])*(h0+h[3]+h[2]) /((h[3]+h[2])*h[2]*h[1])<=h0/h[3]) then begin q:=h0*h0/h[3]; q_prime:=2*h0/h[3] end else if (h0*(h0+h[3])*(h0+h[3]+h[2]) /(h[3]*h[2]*(h[2]+h[1]))<=h0*(h0+h[3])/(h[3]*h[2])) and (h0*(h0+h[3])*(h0+h[3]+h[2]) /((h[3]+h[2])*h[2]*h[1]) <=h0*(h0+h[3])/(h[3]*h[2])) then begin q:=h0*h0*(h0+h[3])/(h[3]*h[2]); q_prime:=h0*(3*h0+2*h[3])/(h[3]*h[2]) end else if h0*(h0+h[3])*(h0+h[3]+h[2]) /((h[3]+h[2])*h[2]*h[1]) <=h0*(h0+h[3])*(h0+h[3]+h[2]) /(h[3]*h[2]*(h[2]+h[1])) then begin q:=h0*h0*(h0+h[3])*(h0+h[3]+h[2]) /(h[3]*h[2]*(h[2]+h[1])); q_prime:=h0*(h0*(4*h0+6*h[3] +3*h[2])+2*h[3] *(h[3]+h[2]))/(h[3]*h[2] *(h[2]+h[1])) end else begin q:=h0*h0*(h0+h[3])*(h0+h[3]+h[2]) /((h[3]+h[2])*h[2]*h[1]); q_prime:=h0*(h0*(4*h0+6*h[3] +3*h[2])+2*h[3] *(h[3]+h[2]))/((h[3]+h[2]) *h[2]*h[1]) end; h1:=h0-(h0*h0*(((12*h0+15*(3*h[3]+2*h[2]+h[1]))*h0 +20*(h[3]*(h[3]+h[2])+(2*h[3]+h[2]) *(h[3]+h[2]+h[1])))*h0 +30*h[3]*(h[3]+h[2])*(h[3]+h[2]+h[1]))*PSI_width/1440 +lambda*q*Yk_widths+int_width(Y[3])-Y_width) /(h0*(12*h0*h0*h0+9*h0*h0*(3*h[3]+2*h[2]+h[1]) +8*h0*(h[3]*(h[3]+h[2])+(2*h[3]+h[2])*(h[3]+h[2]+h[1])) +6*h[3]*(h[3]+h[2])*(h[3]+h[2]+h[1])) *PSI_width/288+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[4]:=h1 else begin h[4]:=0; it:=0 end end end {Find_stepsize}; var it, k, max_it, status : Integer; a, b_left, b_right, eps, lambda, PSI_width, real_h, real_t, w, Y_width : Extended; file_name, left, right, time : string; agree, output : Char; finish : Boolean; h : array of Extended; delta_t, delta_y, ih, ih1, ih2, ih3, T1, T2, T3, T4, Yh, Y1, Y2, Y3, Y4 : interval; Yk : array of interval; results : TextFile; time1, time2 : TDateTime; begin Writeln ('Solving the initial value problem'); Writeln ('y'' = (y-t)/(y+t), y(0) = 4'); Writeln ('by an interval four-step Adams-Bashforth methods of order four'); 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); SetLength (h, 5); SetLength (Yk, 4); h[1]:=0.081746227283888863; h[2]:=h[1]; h[3]:=h[1]; h[4]:=h[1]; // Using the below instruction you can insert your own values for h[k] // for k:=1 to 4 do // repeat // Write (' h[', k, '] = '); // Readln (h[k]); // if h[k]<=0 // then Writeln ('h[', k, '] must be greater then zero') // until h[k]>0; 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:=1; 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'' = (y-t)/(y+t), y(0) = 4'); Writeln (results, 'by an interval version of four-step ' +'Adams-Bashforth method of order four'); 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); for k:=1 to 4 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 fourth order'); Writeln end else begin Write (results, 'Enclosures of the exact solution obtained by the ' +'interval Adams-Bashforth method of fourth order'); Writeln (results, ' ') end; k:=4; real_t:=h[1]+h[2]+h[3]; Str (real_t:26, left); T1:=int_read(left); // You should change the below value of Y1, Y2 and Y3 for your own h[k] left:='4.2313071906453238'; right:='4.2313071906453243'; Y1.a:=left_read(left); Y1.b:=right_read(right); Yk[3]:=Y1; w:=h[1]+h[2]; Str (w:26, left); T2:=int_read(left); left:='4.1571485011046702'; right:='4.1571485011046705'; Y2.a:=left_read(left); Y2.b:=right_read(right); Yk[2]:=Y2; w:=h[1]; Str (w:26, left); T3:=int_read(left); left:='4.0801194662264153'; right:='4.0801194662264155'; Y3.a:=left_read(left); Y3.b:=right_read(right); Yk[1]:=Y3; T4:=i0; Y4:=i4; Yk[0]:=Y4; 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; Writeln; if output='f' then Writeln ('Calculating. Please, wait ...'); time1:=Now; PSI_width:=int_width(PSI4(delta_t,delta_y)); 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[4]; if real_t+real_h