program ABVSS_Example_2_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 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_2 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}; procedure IntervalAB1 (const T1, Y1, ih, delta_t, delta_y : interval; out Y : interval); var PSI, Th : interval; begin 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)) end {IntervalAB1}; procedure IntervalAB2 (const T1, Y1, T2, Y2, ih1, ih, delta_t, delta_y : interval; out Y : interval); var PSI, Th, g : interval; begin 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)) end {IntervalAB2}; procedure IntervalAB3 (const T1, Y1, T2, Y2, T3, Y3, ih2, ih1, ih, delta_t, delta_y : interval; out Y : interval); var PSI, Th, g : interval; begin Th:=iadd(ih1,ih2); 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))) end {IntervalAB3}; 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 p : Integer; 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 order p used to solve } { the initial value problem } { y' = 0.5y, y(0) = 1, } { with the exact solution y = exp(0.5t). To find this step the Newton } { iteration method is used. } { Data: } { p - order of method equal to the number k of method steps (p must } { be equal to 1, 2, 3 or 4), } { 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: } { 1 - if p < 1 or p > 4, } { 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 (p<1) or (p>4) then status:=1 else 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>p) or (status=2) end; if status<>0 then begin h[p]:=0; it:=0 end else begin Yk_widths:=0; for k:=p-1 downto 0 do Yk_widths:=Yk_widths+(k+1)*int_width(Y[k]); it:=0; h0:=h[p]; status:=3; repeat it:=it+1; case p of 1 : h1:=h0-(h0*h0*PSI_width/2+(1+h0*lambda)*Yk_widths-Y_width) /(h0*PSI_width+Y_width/2); 2 : begin if h0/h[p-1]<=1 then begin q:=h0; q_prime:=1; end else begin q:=h0*h0/h[p-1]; q_prime:=2*h0/h[p-1] end; h1:=h0-(h0*h0*(2*h0+3*h[p-1])*PSI_width/12 +lambda*q*Yk_widths+int_width(Y[p-1])-Y_width) /(h0*(h0+h[p-1])*PSI_width/2+lambda*q_prime*Yk_widths) end; 3 : begin if (h0/h[p-1]<=1) and (h0*(h0+h[p-1])/(h[p-1]*h[p-2])<=1) then begin q:=h0; q_prime:=1 end else if h0*(h0+h[p-1])/(h[p-1]*h[p-2])<=h0/h[p-1] then begin q:=h0*h0/h[p-1]; q_prime:=2*h0/h[p-1] end else begin q:=h0*h0*(h0+h[p-1])/(h[p-1]*h[p-2]); q_prime:=h0*(3*h0+2*h[p-1])/(h[p-1]*h[p-2]) end; h1:=h0-(h0*h0*((3*h0+4*(2*h[p-1]+h[p-2]))*h0 +6*h[p-1]*(h[p-1]+h[p-2]))*PSI_width/72 +lambda*q*Yk_widths+int_width(Y[p-1])-Y_width) /(h0*(2*h0*h0+2*h0*(2*h[p-1]+h[p-2])+h[p-1] *(h[p-1]+h[p-2]))*PSI_width/12+lambda*q_prime*Yk_widths) end; 4 : begin if (h0/h[p-1]<=1) and (h0*(h0+h[p-1])/(h[p-1]*h[p-2])<=1) and (h0*(h0+h[p-1])*(h0+h[p-1]+h[p-2]) /(h[p-1]*h[p-2]*(h[p-2]+h[p-3]))<=1) and (h0*(h0+h[p-1])*(h0+h[p-1]+h[p-2]) /((h[p-1]+h[p-2])*h[p-2]*h[p-3])<=1) then begin q:=h0; q_prime:=1 end else if (h0*(h0+h[p-1])/(h[p-1]*h[p-2])<=h0/h[p-1]) and (h0*(h0+h[p-1])*(h0+h[p-1]+h[p-2]) /(h[p-1]*h[p-2]*(h[p-2]+h[p-3]))<=h0/h[p-1]) and (h0*(h0+h[p-1])*(h0+h[p-1]+h[p-2]) /((h[p-1]+h[p-2])*h[p-2]*h[p-3])<=h0/h[p-1]) then begin q:=h0*h0/h[p-1]; q_prime:=2*h0/h[p-1] end else if (h0*(h0+h[p-1])*(h0+h[p-1]+h[p-2]) /(h[p-1]*h[p-2]*(h[p-2]+h[p-3])) <=h0*(h0+h[p-1])/(h[p-1]*h[p-2])) and (h0*(h0+h[p-1])*(h0+h[p-1]+h[p-2]) /((h[p-1]+h[p-2])*h[p-2]*h[p-3]) <=h0*(h0+h[p-1])/(h[p-1]*h[p-2])) then begin q:=h0*h0*(h0+h[p-1])/(h[p-1]*h[p-2]); q_prime:=h0*(3*h0+2*h[p-1]) /(h[p-1]*h[p-2]) end else if h0*(h0+h[p-1])*(h0+h[p-1]+h[p-2]) /((h[p-1]+h[p-2])*h[p-2]*h[p-3]) <=h0*(h0+h[p-1])*(h0+h[p-1]+h[p-2]) /(h[p-1]*h[p-2]*(h[p-2]+h[p-3])) then begin q:=h0*h0*(h0+h[p-1]) *(h0+h[p-1]+h[p-2]) /(h[p-1]*h[p-2] *(h[p-2]+h[p-3])); q_prime:=h0*(h0 *(4*h0+6*h[p-1] +3*h[p-2]) +2*h[p-1] *(h[p-1]+h[p-2])) /(h[p-1]*h[p-2] *(h[p-2]+h[p-3])) end else begin q:=h0*h0*(h0+h[p-1]) *(h0+h[p-1]+h[p-2]) /((h[p-1]+h[p-2]) *h[p-2]*h[p-3]); q_prime:=h0*(h0 *(4*h0+6*h[p-1] +3*h[p-2]) +2*h[p-1] *(h[p-1]+h[p-2])) /((h[p-1]+h[p-2]) *h[p-2]*h[p-3]) end; h1:=h0-(h0*h0*(((12*h0+15*(3*h[p-1]+2*h[p-2]+h[p-3]))*h0 +20*(h[p-1]*(h[p-1]+h[p-2])+(2*h[p-1]+h[p-2]) *(h[p-1]+h[p-2]+h[p-3])))*h0 +30*h[p-1]*(h[p-1]+h[p-2])*(h[p-1]+h[p-2]+h[p-3])) *PSI_width/1440 +lambda*q*Yk_widths+int_width(Y[p-1])-Y_width) /(h0*(12*h0*h0*h0+9*h0*h0*(3*h[p-1]+2*h[p-2]+h[p-3]) +8*h0*(h[p-1]*(h[p-1]+h[p-2]) +(2*h[p-1]+h[p-2])*(h[p-1]+h[p-2]+h[p-3])) +6*h[p-1]*(h[p-1]+h[p-2])*(h[p-1]+h[p-2]+h[p-3])) *PSI_width/288+lambda*q_prime*Yk_widths) end; end; 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[p]:=h1 else begin h[p]:=0; it:=0 end end end {Find_stepsize}; var it, k, max_it, p, 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'' = 0.5y, y(0) = 1'); Writeln ('by interval k-step Adams-Bashforth methods of order p = k (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 (' order p = '); Readln (p); if (p<1) or (p>4) then Writeln ('p must be equal to 1, 2, 3 or 4') until (p>0) and (p<5); 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, p+1); SetLength (Yk, p); for k:=1 to p 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:=0.5; 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'' = 0.5y, y(0) = 1'); Writeln (results, 'by an interval version of ', p, '-step ' +'Adams-Bashforth method of order p = ', p); 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 p 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 Write ('Enclosures of the exact solution obtained by the interval ' +'Adams-Bashforth method of ') else Write (results, 'Enclosures of the exact solution obtained by the ' +'interval Adams-Bashforth method of '); case p of 1 : begin if output='s' then begin Writeln ('first order'); Writeln end else begin Writeln (results, 'first order '); Writeln (results, ' ') end; k:=1; 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; Yk[0]:=Y1; 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 end; 2 : begin if output='s' then begin Writeln ('second order'); Writeln end else begin Writeln (results, 'second order '); Writeln (results, ' ') end; k:=2; real_t:=h[1]; Str (h[1]:26, left); ih1:=int_read(left); T1:=ih1; w:=Exp(0.5*real_t); Str (w:26, left); Y1:=int_read(left); Yk[1]:=Y1; T2:=i0; Y2:=i1; Yk[0]:=Y2; 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 end; 3 : begin if output='s' then begin Writeln ('third order'); Writeln end else begin Writeln (results, 'third order '); Writeln (results, ' ') end; k:=3; real_t:=h[1]+h[2]; Str (h[1]:26, left); ih1:=int_read(left); Str (h[2]:26, left); ih2:=int_read(left); Str (real_t:26, left); T1:=int_read(left); w:=Exp(0.5*real_t); Str (w:26, left); Y1:=int_read(left); Yk[2]:=Y1; T2:=ih1; w:=Exp(0.5*h[1]); Str (w:26, left); Y2:=int_read(left); Yk[1]:=Y2; T3:=i0; Y3:=i1; Yk[0]:=Y3; 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 end; 4 : begin if output='s' then begin Writeln ('fourth order'); Writeln end else begin Writeln (results, 'fourth order '); Writeln (results, ' ') end; k:=4; 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); Yk[3]:=Y1; 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); Yk[2]:=Y2; w:=h[1]; Str (w:26, left); T3:=int_read(left); w:=Exp(0.5*w); Str (w:26,left); Y3:=int_read(left); Yk[1]:=Y3; T4:=i0; Y4:=i1; 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 end end; Writeln; if output='f' then Writeln ('Calculating. Please, wait ...'); time1:=Now; case p of 1 : PSI_width:=int_width(PSI1(delta_t,delta_y)); 2 : PSI_width:=int_width(PSI2(delta_t,delta_y)); 3 : PSI_width:=int_width(PSI3(delta_t,delta_y)); 4 : PSI_width:=int_width(PSI4(delta_t,delta_y)) end; finish:=False; repeat Find_stepsize (p, Yk, lambda, PSI_width, Y_width, eps, max_it, h, it, status); if status=0 then begin real_h:=h[p]; if real_t+real_h