program IMM_Example_4; {$APPTYPE CONSOLE} uses System.SysUtils, IntervalArithmetic32and64; // The program solves the simple pendulum motion problem given by the initial // value problem y'[1] = -u^2*y[2], y'[2] = y[1], y[1](0) = 0, y[2] = phi0, // where u = Sqrt(g/L), g = 9.80665 , L = 1, phi0 = Pi/6, by four-step // interval method of Milne type. It is assume h = 0.001, and all starting // intervals are defined within the program. // 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 array2interval = array [1..2] of interval; function FTY (const u2, T : interval; const Y : array2interval) : array2interval; begin Result[1]:=-u2*Y[2]; Result[2]:=Y[1] end {FTY}; function PSI (const u2, T : interval; const Y : array2interval) : array2interval; var u4 : interval; begin u4:=u2*u2; Result[1]:=-u4*u2*Y[2]; Result[2]:=u4*Y[1] end {PSI}; var i, n, m : Integer; w : Extended; delta_t, h, h5, hh, T, T0, T1, T2, T3, T4, u2 : Interval; i1, i2, i3, i4, delta_y, Fdelta, Y, Y0, Y1, Y2, Y3, Y4 : array2interval; left, right, s : string; find : Boolean; z : Char; plik : TextFile; procedure Results (const T : interval; const Y : array2interval); var hreal : Extended; hrealstring : string; begin hreal:=(T.a+T.b)/2; Str (hreal:4:2, hrealstring); iends_to_strings (T, left, right); if z='f' then Writeln (plik, ' T('+hrealstring+') = ['+left+','+right+']') else Writeln (' T('+hrealstring+') = ['+left+','+right+']'); if z='f' then begin iends_to_strings (Y[1], left, right); Writeln (plik, 'Y[1]('+hrealstring+') = ['+left+','+right+']'); w:=int_width(Y[1]); Writeln (plik, 'width = ', w:11); iends_to_strings (Y[2], left, right); Writeln (plik, 'Y[2]('+hrealstring+') = ['+left+','+right+']'); w:=int_width(Y[2]); Writeln (plik, 'width = ', w:11) end else begin iends_to_strings (Y[1], left, right); Writeln ('Y[1]('+hrealstring+') = ['+left+','+right+']'); w:=int_width(Y[1]); Writeln ('width = ', w:11); iends_to_strings (Y[2], left, right); Writeln ('Y[2]('+hrealstring+') = ['+left+','+right+']'); w:=int_width(Y[2]); Writeln ('width = ', w:11) end end {Results}; begin repeat Write ('Output (s - screen, f - file) : '); Readln (z) until (z='s') or (z='f'); if z='f' then begin AssignFile (plik, 'IMM4RESULTS.TXT'); Rewrite (plik); Writeln ('Output to IMM4RESULTS.TXT in current directory') end else Writeln; m:=2000; s:='9.80665'; u2.a:=left_read(s); u2.b:=right_read(s); s:='0.001'; h.a:=left_read(s); h.b:=right_read(s); delta_t.a:=0; delta_t.b:=2; s:='-1.8'; delta_y[1].a:=left_read(s); s:='1.8'; delta_y[1].b:=right_read(s); s:='-0.6'; delta_y[2].a:=left_read(s); s:='0.6'; delta_y[2].b:=right_read(s); Fdelta:=FTY(u2,delta_t,delta_y); w:=Sqrt(9.80665); if z='s' then begin Writeln ('Initial value problem: y''[1] = -u^2*y[2], y''[2] = ' +'y[1],'); Writeln (' y[1](0) = 0, y[2](0) = ' +'phi0,'); Writeln ('where u = Sqrt(g/L), g = 9.80665, L = 1, phi0 = Pi/6'); Writeln; Writeln ('Exact solution:'); Writeln ('y[1](0.5) = ', -Pi*w*sin(0.5*w)/6:26, ' y[2](0.5) = ', Pi*cos(0.5*w)/6:26); Writeln ('y[1](1.0) = ', -Pi*w*sin(w)/6:26, ' y[2](1.0) = ', Pi*cos(w)/6:26); Writeln ('y[1](1.5) = ', -Pi*w*sin(1.5*w)/6:26, ' y[2](1.5) = ', Pi*cos(1.5*w)/6:26); Writeln ('y[1](2.0) = ', -Pi*w*sin(2*w)/6:26, ' y[2](2.0) = ', Pi*cos(2*w)/6:26); Writeln; Readln end else begin Writeln (plik, 'Initial value problem: y''[1] = -u^2*y[2], y''[2] = ' +'y[1],'); Writeln (plik, ' y[1](0) = 0, y[2](0) = ' +'phi0,'); Writeln (plik, 'where u = Sqrt(g/L), g = 9.80665, L = 1, ' +'phi0 = Pi/6'); Writeln (plik, ' '); Writeln (plik, 'Exact solution:'); Writeln (plik, '==============='); Writeln (plik, 'y[1](0.5) = ', -Pi*sin(0.5*w)/6:26, ' y[2](0.5) = ', Pi*cos(0.5*w)/6:26); Writeln (plik, 'y[1](1.0) = ', -Pi*sin(w)/6:26, ' y[2](1.0) = ', Pi*cos(w)/6:26); Writeln (plik, 'y[1](1.5) = ', -Pi*sin(1.5*w)/6:26, ' y[2](1.5) = ', Pi*cos(1.5*w)/6:26); Writeln (plik, 'y[1](2.0) = ', -Pi*sin(2*w)/6:26, ' y[2](2.0) = ', Pi*cos(2*w)/6:26) end; T0:=0; Y0[1]:=0; s:='0.52359877559829887'; Y0[2].a:=left_read(s); s:='0.52359877559829888'; Y0[2].b:=right_read(s); T1:=T0+h; s:='-0.00513474154027594'; Y1[1].a:=left_read(s); s:='-0.00513474154027593'; Y1[1].b:=right_read(s); s:='0.52359620822543062'; Y1[2].a:=left_read(s); s:='0.52359620822543063'; Y1[2].b:=right_read(s); T2:=T1+h; s:='-0.01026943272597990'; Y2[1].a:=left_read(s); s:='-0.01026943272597989'; Y2[1].b:=right_read(s); s:='0.52358850613200318'; Y2[2].a:=left_read(s); s:='0.52358850613200319'; Y2[2].b:=right_read(s); s:='-0.01540402320303372'; T3:=T2+h; Y3[1].a:=left_read(s); s:='-0.01540402320303371'; Y3[1].b:=right_read(s); s:='0.52357566939354822'; Y3[2].a:=left_read(s); s:='0.52357566939354823'; Y3[2].b:=right_read(s); s:='-0.003'; hh.a:=left_read(s); s:='0.001'; hh.b:=right_read(s); h5:=h*h; h5:=h5*h5; h5:=h5*h; for n:=4 to m do begin T:=T1+hh; for i:=1 to 2 do Y[i]:=Y1[i]+hh*Fdelta[i]; i1:=PSI(u2,T,Y); for i:=1 to 2 do begin i2[i]:=27*i1[i]-251*i1[i]; i1[i]:=h5*i2[i]/720 end; i2:=FTY(u2,T3,Y3); i3:=FTY(u2,T2,Y2); i4:=FTY(u2,T1,Y1); T4:=T3+h; for i:=1 to 2 do begin i2[i]:=2*i2[i]-i3[i]+2*i4[i]; i2[i]:=4*h*i2[i]/3; Y4[i]:=Y0[i]+(i2[i]-i1[i]) end; if (n=500) or (n=1000) or (n=1500) or (n=2000) then Results (T4, Y4); if n