program IPCM_Example_7_LongTime; {$APPTYPE CONSOLE} uses System.SysUtils, IntervalArithmetic32and64; // The program solves a special case of Hill equations given by the initial // value problem y'[l] = y[l] (l = 1, 2), y'[3] = -(kappa/r^3)y[1], // y'[4] = -(kappa/r^3)y[2], y[1](0) = y[4](0) = 1, y[2](0) = y[3](0) = 0, // where r = Sqrt(y[1]^2 + y[2]^2), kappa = 1, by three-step interval // predictor-corrector method of Nystrom-Milne-Simpson type. It is assumed // h = 0.0005 and 0 <= t <= 16. // All starting intervals are defined within the program. // You can direct output to the screen (s) or to the file (f). // Directing output to the screen, press Enter after each result displayed. // 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 array4interval = array [1..4] of interval; function Psi3Explicit (const T : interval; const Y : array4interval) : array4interval; var psi, r2, r4, r5, r6, r7, r8, sqr_y1, sqr_y2, sqr_y3, sqr_y4, sqr2_y1, sqr2_y2, w1, w2, w3, w4, w5, y1y2, y3y4 : interval; begin sqr_y1:=Y[1]*Y[1]; sqr_y2:=Y[2]*Y[2]; r2:=sqr_y1+sqr_y2; w1.a:=Sqrt(r2.a); w1.b:=Sqrt(r2.b); r4:=r2*r2; r6:=r4*r2; r8:=r4*r4; r5:=w1*r4; r7:=w1*r6; sqr_y3:=Y[3]*Y[3]; sqr_y4:=Y[4]*Y[4]; y1y2:=Y[1]*Y[2]; y3y4:=Y[3]*Y[4]; w1:=3-5*sqr_y1/r2; w2:=1-5*sqr_y1/r2; w3:=1-5*sqr_y2/r2; w4:=1-3*sqr_y1/r2; psi:=(3*Y[1]/r5)*(w1*sqr_y3); w1:=(6*Y[2]/r5)*(w2*y3y4); psi:=psi+w1; w1:=(3*Y[1]/r5)*(w3*sqr_y4); psi:=psi+w1; w1:=(1/r6)*(w4*Y[1]); psi:=psi+w1; w1:=(3*y1y2/r8)*Y[2]; Result[1]:=psi-w1; w1:=3-5*sqr_y2/r2; w4:=1-3*sqr_y2/r2; psi:=(3*Y[2]/r5)*(w1*sqr_y4); w1:=(3*Y[2]/r5)*(w2*sqr_y3); psi:=psi+w1; w1:=(6*Y[1]/r5)*(w3*y3y4); psi:=psi+w1; w1:=(3*y1y2/r8)*Y[1]; psi:=psi-w1; w1:=(1/r6)*(w4*Y[2]); Result[2]:=psi+w1; w1:=35*(sqr_y1*sqr_y2)/r4; w1:=w1-5*(sqr_y1+sqr_y2)/r2; w1:=1+w1; w2:=9*(sqr_y1+sqr_y2)/r2; w2:=w2-4; w3:=15*y1y2/r7; w4:=6*y1y2/r8; sqr2_y1:=sqr_y1*sqr_y1; sqr2_y2:=sqr_y2*sqr_y2; w5:=3-30*sqr_y1/r2; w5:=w5+35*sqr2_y1/r4; psi:=(3/r5)*(w5*(sqr_y3*Y[3])); w5:=21*sqr_y1/r2-9; w5:=w3*(w5*(sqr_y3*Y[4])); psi:=psi+w5; w5:=(9/r5)*(w1*(Y[3]*sqr_y4)); psi:=psi+w5; w5:=7*sqr_y2/r2-3; w5:=w3*(w5*(sqr_y4*Y[4])); psi:=psi+w5; w5:=54*(sqr_y1*sqr_y2)/r4; w5:=54*sqr2_y1/r4+w5; w5:=w5-9*sqr_y2/r2; w5:=w5-33*sqr_y1/r2; w5:=1+w5; w5:=(1/r6)*(w5*Y[3]); psi:=psi+w5; w5:=w4*(w2*Y[4]); Result[3]:=psi+w4; w5:=7*sqr_y1/r2-3; psi:=w3*w5; w5:=(9/r5)*w1; w5:=w5*(sqr_y3*Y[4]); psi:=psi+w5; w5:=21*sqr_y2/r2-9; psi:=psi+w3*w5; w5:=35*sqr2_y2/r4; w5:=w5-30*sqr_y2/r2; w5:=(3/r5)*(3+w5); psi:=psi+(w5*(sqr_y4*Y[4])); w5:=w4*(w2*Y[3]); psi:=psi+w5; w5:=54*sqr2_y2/r4; w5:=w5+54*(sqr_y1*sqr_y2)/r4; w5:=w5-33*sqr_y2/r2; w5:=w5-9*sqr_y1/r2; w5:=((1/r6)*(1+w5))*Y[4]; Result[4]:=psi+w5 end {Psi3Explicit}; function Psi4Explicit (const T : interval; const Y : array4interval) : array4interval; var psi, r2, r4, r5, r6, r7, r8, r9, sqr_y1, sqr_y2, sqr_y3, sqr_y4, sqr2_y1, sqr2_y2, sqr2_y3, sqr2_y4, w1, w2, w3, w4, w5, w6, y1y2, y3y4 : interval; begin sqr_y1:=Y[1]*Y[1]; sqr_y2:=Y[2]*Y[2]; r2:=sqr_y1+sqr_y2; w1.a:=Sqrt(r2.a); w1.b:=Sqrt(r2.b); r4:=r2*r2; r6:=r4*r2; r8:=r4*r4; r5:=w1*r4; r7:=w1*r6; r9:=w1*r8; sqr_y3:=Y[3]*Y[3]; sqr_y4:=Y[4]*Y[4]; y1y2:=Y[1]*Y[2]; y3y4:=Y[3]*Y[4]; w1:=35*(sqr_y1*sqr_y2)/r4; w1:=w1-5*(sqr_y1+sqr_y2)/r2; w1:=1+w1; w2:=9*(sqr_y1+sqr_y2)/r2; w2:=w2-4; w3:=15*y1y2/r7; w4:=6*y1y2/r8; sqr2_y1:=sqr_y1*sqr_y1; sqr2_y2:=sqr_y2*sqr_y2; sqr2_y3:=sqr_y3*sqr_y3; sqr2_y4:=sqr_y4*sqr_y4; w5:=3-30*sqr_y1/r2; w5:=w5+35*sqr2_y1/r4; psi:=(3/r5)*(w5*(sqr_y3*Y[3])); w5:=21*sqr_y1/r2-9; w5:=w3*(w5*(sqr_y3*Y[4])); psi:=psi+w5; w5:=(9/r5)*(w1*(Y[3]*sqr_y4)); psi:=psi+w5; w5:=7*sqr_y2/r2-3; w5:=w3*(w5*(sqr_y4*Y[4])); psi:=psi+w5; w5:=54*(sqr_y1*sqr_y2)/r4; w5:=54*sqr2_y1/r4+w5; w5:=w5-9*sqr_y2/r2; w5:=w5-33*sqr_y1/r2; w5:=1+w5; w5:=(1/r6)*(w5*Y[3]); psi:=psi+w5; w5:=w4*(w2*Y[4]); Result[1]:=psi+w4; w5:=7*sqr_y1/r2-3; psi:=w3*(w5*(sqr_y3*Y[3])); w5:=(9/r5)*w1; w5:=w5*(sqr_y3*Y[4]); psi:=psi+w5; w5:=21*sqr_y2/r2-9; psi:=psi+w3*(w5*(Y[3]*sqr_y4)); w5:=35*sqr2_y2/r4; w5:=w5-30*sqr_y2/r2; w5:=(3/r5)*(3+w5); psi:=psi+w5*(sqr_y4*Y[4]); w5:=w4*(w2*Y[3]); psi:=psi+w5; w5:=54*sqr2_y2/r4; w5:=w5+54*(sqr_y1*sqr_y2)/r4; w5:=w5-33*sqr_y2/r2; w5:=w5-9*sqr_y1/r2; w5:=(1/r6)*(1+w5)*Y[4]; Result[2]:=psi+w5; w1:=14*sqr_y1/r2; w1:=w1-21*sqr2_y1/r4; w1:=w1-1; w2:=14*sqr_y2/r2; w2:=w1-21*sqr2_y2/r4; w2:=w1-1; w3:=7*sqr_y1/r2; w3:=w3+21*sqr_y2/r2; w4:=63*(sqr_y1*sqr_y2)/r4; w3:=w3-w4-3; w4:=21*sqr_y1/r2-w4; w4:=w4+7*sqr_y2/r2; w4:=w4-3; w5:=33*(sqr_y1+sqr_y2)/r2; w5:=w5-54*(sqr2_y1+sqr2_y2)/r4; w5:=w5-108*(sqr_y1*sqr_y2)/r4; w5:=w5-1; w6:=70*sqr_y1/r2; w6:=w6-63*sqr2_y1/r4; w6:=(w6-15)*sqr2_y3; psi:=(15*Y[1]/r7)*w6; w6:=w1*((sqr_y3*Y[3])*Y[4]); w6:=(180*Y[2]/r7)*w6; psi:=psi+w6; w6:=w3*(sqr_y3*sqr_y4); w6:=(90*Y[1]/r7)*w6; psi:=psi+w6; w6:=w4*(Y[3]*(sqr_y4*Y[4])); w6:=(60*Y[2]/r7)*w6; psi:=psi+w6; w6:=w2*sqr2_y4; w6:=(45*Y[1]/r7)*w6; psi:=psi+w6; w6:=(250*sqr_y1+105*sqr_y2)/r2; w6:=w6-285*sqr2_y1/r4; w6:=w6-285*(sqr_y1*sqr_y2)/r4; w6:=(w6-33)*sqr_y3; w6:=3*Y[1]*w6/r8; psi:=psi+w6; w6:=(181*sqr_y1+36*sqr_y2)/r2; w6:=w6-285*sqr2_y1/r4; w6:=w6-285*(sqr_y1*sqr_y2)/r4; w6:=w6*(Y[3]*Y[4]); w6:=6*Y[2]*w6/r8; psi:=psi+w6; w6:=(33*sqr_y1+178*sqr_y2)/r2; w6:=w6-285*(sqr_y1*sqr_y2)/r4; w6:=w6-285*sqr2_y2/r4; w6:=(w6-11)*sqr2_y4; w6:=3*Y[1]*w6/r8; psi:=psi+w6; w6:=Y[1]*w5/r9; Result[3]:=psi+w6; w6:=w1*sqr2_y3; psi:=45*(Y[2]*w6)/r7; w6:=w3*((sqr_y3*Y[3])*Y[4]); w6:=60*(Y[1]*w6)/r7; psi:=psi+w6; w6:=w4*(sqr_y3*sqr_y3); w6:=90*(Y[2]*w6)/r7; psi:=psi+w6; w6:=w2*(Y[3]*(sqr_y4*Y[4])); w6:=180*(Y[1]*w6)/r7; psi:=psi+w6; w6:=70*sqr_y2/r2; w6:=w6-63*sqr2_y2/r4; w6:=(w6-15)*sqr2_y4; w6:=15*(Y[2]*w6)/r7; psi:=psi+w6; w6:=(178*sqr_y1+33*sqr_y2)/r2; w6:=w6-285*sqr2_y1/r4; w6:=w6-285*(sqr_y1*sqr_y2)/r4; w6:=(w6-11)*sqr2_y3; w6:=3*(Y[2]*w6)/r8; psi:=psi+w6; w6:=(36*sqr_y1+181*sqr_y2)/r2; w6:=w6-285*(sqr_y1*sqr_y2)/r4; w6:=w6-285*sqr2_y2/r4; w6:=(w6-11)*(Y[3]*Y[4]); w6:=6*(Y[1]*w6)/r8; psi:=psi+w6; w6:=(105*sqr_y1+250*sqr_y2)/r2; w6:=w6-285*(sqr_y1*sqr_y2)/r4; w6:=w6-285*sqr2_y2/r4; w6:=(w6-33)*sqr2_y4; w6:=3*(Y[2]*w6)/r8; psi:=psi+w6; w6:=Y[2]*w5/r9; Result[4]:=psi+w6 end {Psi4Explicit}; function Psi3Implicit (const T : interval; const Y : array4interval) : array4interval; begin Result:=Psi4Explicit(T,Y) end {Psi3Implicit}; function FTY (const T : interval; const Y : array4interval) : array4interval; var r, r2 : interval; begin Result[1]:=Y[3]; Result[2]:=Y[4]; r2:=Y[1]*Y[1]+Y[2]*Y[2]; r.a:=Sqrt(r2.a); r.b:=Sqrt(r2.b); r:=r2*r; Result[3]:=-Y[1]/r; Result[4]:=-Y[2]/r end {FTY}; var i, k, l, n, m, max_l : Integer; w : Extended; delta_t, h, hh, hh1, T, T0, T1, T2, T3, T4 : interval; delta_y, i1, i2, i3, i4, Fdelta, Y, Y0, Y1, Y2, Y3, Y4 : array4interval; left, right, s : string; find : Boolean; kind, z : Char; plik : TextFile; procedure Results (const T : interval; const Y : array4interval; const kind : Char); var i : Integer; hreal : Extended; hrealstring : string; begin hreal:=(T.a+T.b)/2; Str (hreal:4:1, hrealstring); iends_to_strings (T, left, right); if z='f' then Writeln (plik, 'T('+hrealstring+') = ['+left+','+right+']') else Writeln ('T('+hrealstring+') = ['+left+','+right+']'); for i:=1 to 4 do begin iends_to_strings (Y[i], left, right); if z='f' then Writeln (plik, 'Y[', i, ']('+hrealstring+') = ['+left+','+right+']') else Writeln ('Y[', i, ']('+hrealstring+') = ['+left+','+right+']'); w:=int_width(Y[i]); if z='f' then if (kind='i') and (i=4) then Writeln (plik, 'width = ', w:11, ', max number of iterations : ', max_l) else Writeln (plik, 'width = ', w:11) else if (kind='i') and (i=4) then Writeln ('width = ', w:11, ', max number of iterations : ', max_l) else Writeln ('width = ', w:11) end; if z='s' then Readln else Writeln (plik, ' ') end {Results}; begin repeat Write ('Output (s - screen, f - file) : '); Readln (z) until (z='s') or (z='f'); Writeln; Writeln ('Be patient ... Calculations may take a lot of minutes.'); if z='f' then begin AssignFile (plik, 'IPCM7LTRESULTS.TXT'); Rewrite (plik); Writeln ('Output to IPCM7LTRESULTS.TXT in current directory'); Writeln ('Please wait ...'); Writeln (plik, ' '); Writeln (plik, 'k = 3, Nystrom-Milne-Simpson predictor-corrector'); Writeln (plik, '================================================'); Writeln (plik, ' ') end else begin Writeln; Writeln ('k = 3, Nystrom-Milne-Simpson predictor-corrector'); Writeln end; m:=2000; s:='0.0005'; h.a:=left_read(s); h.b:=right_read(s); s:='-0.001'; hh.a:=left_read(s); s:='0.0005'; hh.b:=right_read(s); s:='-0.0015'; hh1.a:=left_read(s); hh1.b:=0; k:=0; repeat case k of 0 : delta_t.a:=0; 1 : delta_t.a:=1; 2 : delta_t.a:=2; 3 : delta_t.a:=3; 4 : delta_t.a:=4; 5 : delta_t.a:=5; 6 : delta_t.a:=6; 7 : delta_t.a:=7; 8 : delta_t.a:=8; 9 : delta_t.a:=9; 10 : delta_t.a:=10; 11 : delta_t.a:=11; 12 : delta_t.a:=12; 13 : delta_t.a:=13; 14 : delta_t.a:=14; 15 : delta_t.a:=15; end; case k of 0 : s:='1.001'; 1 : s:='2.001'; 2 : s:='3.001'; 3 : s:='4.001'; 4 : s:='5.001'; 5 : s:='6.001'; 6 : s:='7.001'; 7 : s:='8.001'; 8 : s:='9.001'; 9 : s:='10.001'; 10 : s:='11.001'; 11 : s:='12.001'; 12 : s:='13.001'; 13 : s:='14.001'; 14 : s:='15.001'; 15 : s:='16.001' end; delta_t.b:=right_read(s); if (k=3) or (k=8) or (k=9) or (k=15) then delta_y[1].a:=-1 else if k=11 then delta_y[1].a:=0 else begin case k of 0 : s:='0.4'; 1 : s:='-0.5'; 2 : s:='-1.1'; 4 : s:='-0.6'; 5 : s:='0.2'; 6 : s:='0.7'; 7 : s:='-0.2'; 10 : s:='-0.9'; 12 : s:='0.8'; 13 : s:='0.1'; 14 : s:='-0.8' end; delta_y[1].a:=left_read(s) end; if (k=5) or (k=6) or (k=12) or (k=13) then delta_y[1].b:=1 else begin case k of 0 : s:='1.1'; 1 : s:='0.6'; 2 : s:='-0.3'; 3 : s:='-0.6'; 4 : s:='0.4'; 7 : s:='0.8'; 8 : s:='-0.1'; 9 : s:='-0.8'; 10 : s:='0.1'; 11 : s:='0.9'; 14 : s:='0.2'; 15 : s:='-0.7' end; delta_y[1].b:=right_read(s) end; delta_y[4]:=delta_y[1]; if (k=4) or (k=5) or (k=10) or (k=11) then delta_y[2].a:=-1 else if k=2 then delta_y[2].a:=0 else begin case k of 0 : s:='-0.1'; 1 : s:='0.8'; 3 : s:='-0.9'; 6, 15 : s:='-0.3'; 7, 14 : s:='0.6'; 8, 13 : s:='0.4'; 9, 12 : s:='-0.6' end; delta_y[2].a:=left_read(s) end; if (k=0) or (k=2) or (k=7) or (k=8) or (k=13) or (k=14) then delta_y[2].b:=1 else if (k=9) or (k=12) then delta_y[2].b:=0.5 else if (k=10) or (k=11) then delta_y[2].b:=-0.5 else begin case k of 1 : s:='1.1'; 3 : s:='0.2'; 4 : s:='-0.7'; 5 : s:='-0.2'; 6, 15 : s:='0.7' end; delta_y[2].b:=right_read(s) end; if (k<=2) or (k=7) or (k=8) or (k=13) or (k=14) then delta_y[3].a:=-1 else if (k=10) or (k=11) then delta_y[3].a:=0.5 else if k=12 then delta_y[3].a:=-0.5 else begin case k of 3 : s:='-0.2'; 4 : s:='0.7'; 5 : s:='0.2'; 6, 15 : s:='-0.7'; 9 : s:='-0.4'; end; delta_y[3].a:=left_read(s) end; if (k=4) or (k=5) or (k=10) or (k=11) then delta_y[3].b:=1 else if k=2 then delta_y[3].b:=0 else begin case k of 0 : s:='0.1'; 1 : s:='-0.7'; 3 : s:='0.9'; 6, 15 : s:='0.3'; 7, 14 : s:='-0.6'; 8 : s:='-0.4'; 9, 12 : s:='0.6'; 13 : s:='-0.4' end; delta_y[3].b:=right_read(s) end; Fdelta:=FTY(delta_t,delta_y); case k of 0 : T0:=0; 1 : T0:=1; 2 : T0:=2; 3 : T0:=3; 4 : T0:=4; 5 : T0:=5; 6 : T0:=6; 7 : T0:=7; 8 : T0:=8; 9 : T0:=9; 10 : T0:=10; 11 : T0:=11; 12 : T0:=12; 13 : T0:=13; 14 : T0:=14; 15 : T0:=15 end; if k=0 then begin Y0[1]:=1; Y0[2]:=0; Y0[3]:=0; Y0[4]:=1 end else begin case k of 1 : s:='0.54030230586806193'; 2 : s:='-0.41614683654775671'; 3 : s:='-0.98999249660561484'; 4 : s:='-0.65364362090705469'; 5 : s:='0.28366218524571665'; 6 : s:='0.96017028528834229'; 7 : s:='0.75390224306024896'; 8 : s:='-0.14550009494259557'; 9 : s:='-0.91113059596678088'; 10 : s:='-0.83907433930908552'; 11 : s:='0.00440885388991498'; 12 : s:='0.84377178449471551'; 13 : s:='0.90675155978951351'; 14 : s:='0.13216376315652315'; 15 : s:='-0.77989719797627575' end; Y0[1].a:=left_read(s); case k of 1 : s:='0.54030230586832392'; 2 : s:='-0.41614683654635373'; 3 : s:='-0.98999249659487845'; 4 : s:='-0.65364362081898057'; 5 : s:='0.28366218568220511'; 6 : s:='0.96017028801279997'; 7 : s:='0.75390226562544140'; 8 : s:='-0.14549997267595006'; 9 : s:='-0.91112992780268283'; 10 : s:='-0.83906871883065234'; 11 : s:='0.00444254242689546'; 12 : s:='0.84393612041229826'; 13 : s:='0.90814130320436176'; 14 : s:='0.14128459834577134'; 15 : s:='-0.73884604744310250' end; Y0[1].b:=right_read(s); case k of 1 : s:='0.84147098480775095'; 2 : s:='0.90929742682451305'; 3 : s:='0.14112000805250249'; 4 : s:='-0.75680249533945081'; 5 : s:='-0.95892427491974987'; 6 : s:='-0.27941550000768980'; 7 : s:='0.65698659093572109'; 8 : s:='0.98935818378424983'; 9 : s:='0.41211801563208774'; 10 : s:='-0.54402314343935587'; 11 : s:='-1.0000057070508080'; 12 : s:='-0.53669349945768798'; 13 : s:='0.41962116724074424'; 14 : s:='0.98679166313019064'; 15 : s:='0.61903860843909704' end; Y0[2].a:=left_read(s); case k of 1 : s:='0.84147098480799766'; 2 : s:='0.90929742682671995'; 3 : s:='0.14112000806669103'; 4 : s:='-0.75680249527699240'; 5 : s:='-0.95892427440600963'; 6 : s:='-0.27941549638871765'; 7 : s:='0.65698660650295722'; 8 : s:='0.98935830946220317'; 9 : s:='0.41211895484934122'; 10 : s:='-0.54401907833988212'; 11 : s:='-0.99997470562124056'; 12 : s:='-0.53645231613732828'; 13 : s:='0.42071296838809284'; 14 : s:='0.99439485439173265'; 15 : s:='0.68020167857772875' end; Y0[2].b:=right_read(s); case k of 1 : s:='-0.84147098480813335'; 2 : s:='-0.90929742682673177'; 3 : s:='-0.14112000807149823'; 4 : s:='0.75680249522115091'; 5 : s:='0.95892427437251778'; 6 : s:='0.27941549519945945'; 7 : s:='-0.65698662170412729'; 8 : s:='-0.98935832664528820'; 9 : s:='-0.41211921023087121'; 10 : s:='0.54401526049331531'; 11 : s:='0.99996685508387995'; 12 : s:='0.53640123583206923'; 13 : s:='-0.42162827395353277'; 14 : s:='-0.99755851328019976'; 15 : s:='-0.68857745643195872' end; Y0[3].a:=left_read(s); case k of 1 : s:='-0.84147098480752950'; 2 : s:='-0.90929742682461389'; 3 : s:='-0.14112000804768938'; 4 : s:='0.75680249539553139'; 5 : s:='0.95892427495332990'; 6 : s:='0.27941550119692548'; 7 : s:='-0.65698657573443962'; 8 : s:='-0.98935816660123022'; 9 : s:='-0.41211776024987189'; 10 : s:='0.54402696133480519'; 11 : s:='1.0000135588090740'; 12 : s:='0.53674453493846494'; 13 : s:='-0.41870865920767607'; 14 : s:='-0.98372616583052569'; 15 : s:='-0.60818530301414816' end; Y0[3].b:=right_read(s); case k of 1 : s:='0.54030230586781509'; 2 : s:='-0.41614683654955448'; 3 : s:='-0.98999249661136999'; 4 : s:='-0.65364362092317612'; 5 : s:='0.28366218492638227'; 6 : s:='0.96017028367999888'; 7 : s:='0.75390224057570812'; 8 : s:='-0.14550016504118860'; 9 : s:='-0.91113108954202910'; 10 : s:='-0.83907475245953394'; 11 : s:='0.00439289421692623'; 12 : s:='0.84362853489013391'; 13 : s:='0.90666758855457258'; 14 : s:='0.12842327253346220'; 15 : s:='-0.82233747023865439' end; Y0[4].a:=left_read(s); case k of 1 : s:='0.54030230586841365'; 2 : s:='-0.41614683654496195'; 3 : s:='-0.98999249659005160'; 4 : s:='-0.65364362080341119'; 5 : s:='0.28366218600140824'; 6 : s:='0.96017028962111593'; 7 : s:='0.75390226810983664'; 8 : s:='-0.14549990257768570'; 9 : s:='-0.91112943422946400'; 10 : s:='-0.83906830567000993'; 11 : s:='0.00445850326738771'; 12 : s:='0.84407945562812538'; 13 : s:='0.90822503300187967'; 14 : s:='0.14494281245391239'; 15 : s:='-0.70217556215033967' end; Y0[4].b:=right_read(s) end; T1:=T0+h; case k of 0 : s:='0.99999987500000260'; 1 : s:='0.53988150285540168'; 2 : s:='-0.41660143322387290'; 3 : s:='-0.99006293285765106'; 4 : s:='-0.65326513796975984'; 5 : s:='0.28414161190515307'; 6 : s:='0.96030987300883725'; 7 : s:='0.75357365552530015'; 8 : s:='-0.14599475589781258'; 9 : s:='-0.91133654167224251'; 10 : s:='-0.83880222680768283'; 11 : s:='0.0049088367412610'; 12 : s:='0.84403987956798485'; 13 : s:='0.90654063185008491'; 14 : s:='0.13166496602099398'; 15 : s:='-0.78024140337658676' end; Y1[1].a:=left_read(s); case k of 0 : s:='0.99999987500000261'; 1 : s:='0.53988150285566393'; 2 : s:='-0.41660143322246896'; 3 : s:='-0.99006293284690298'; 4 : s:='-0.65326513788159893'; 5 : s:='0.28414161234193124'; 6 : s:='0.96030987573629385'; 7 : s:='0.75357367811348522'; 8 : s:='-0.14599463355112602'; 9 : s:='-0.91133587278273002'; 10 : s:='-0.83879660047593158'; 11 : s:='0.00494254863440416'; 12 : s:='0.84420438723870621'; 13 : s:='0.90793183582400716'; 14 : s:='0.14079271895492912'; 15 : s:='-0.73915003144691220' end; Y1[1].b:=right_read(s); case k of 0 : s:='0.00049999997916666'; 1 : s:='0.84174103076555766'; 2 : s:='0.90908923975273205'; 3 : s:='0.14062499418482109'; 4 : s:='-0.75712922253598473'; 5 : s:='-0.95878232396766428'; 6 : s:='-0.27893537995891652'; 7 : s:='0.65736345991697643'; 8 : s:='0.98928531003494974'; 9 : s:='0.41166239859135305'; 10 : s:='-0.54444261279643568'; 11 : s:='-1.0000033856128163'; 12 : s:='-0.53627161817722694'; 13 : s:='0.42007444828580482'; 14 : s:='0.98685574878802483'; 15 : s:='0.61862734483609187' end; Y1[2].a:=left_read(s); case k of 0 : s:='0.00049999997916667'; 1 : s:='0.84174103076580463'; 2 : s:='0.90908923975494112'; 3 : s:='0.14062499419902004'; 4 : s:='-0.75712922247346684'; 5 : s:='-0.95878232345338707'; 6 : s:='-0.27893537633697460'; 7 : s:='0.65736347549798555'; 8 : s:='0.98928543584419948'; 9 : s:='0.41166333863653429'; 10 : s:='-0.54443854447163166'; 11 : s:='-0.99997235136312896'; 12 : s:='-0.53603020931444261'; 13 : s:='0.42116702856286073'; 14 : s:='0.99446720403923067'; 15 : s:='0.67985052112590267' end; Y1[2].b:=right_read(s); case k of 0 : s:='-0.00049999997916667'; 1 : s:='-0.84174103076594050'; 2 : s:='-0.90908923975495312'; 3 : s:='-0.14062499420383427'; 4 : s:='0.75712922241758130'; 5 : s:='0.95878232341992252'; 6 : s:='0.27893537514590531'; 7 : s:='-0.65736349071332147'; 8 : s:='-0.98928545300775326'; 9 : s:='-0.41166359445205973'; 10 : s:='0.54443472247103245'; 11 : s:='0.99996450869843526'; 12 : s:='0.53597903085612526'; 13 : s:='-0.42208346654579323'; 14 : s:='-0.99762997136747911'; 15 : s:='-0.68824405810930523' end; Y1[3].a:=left_read(s); case k of 0 : s:='-0.00049999997916666'; 1 : s:='-0.84174103076533594'; 2 : s:='-0.90908923975283270'; 3 : s:='-0.14062499418000074'; 4 : s:='0.75712922259210957'; 5 : s:='0.95878232400121698'; 6 : s:='0.27893538114996332'; 7 : s:='-0.65736344470152891'; 8 : s:='-0.98928529287146138'; 9 : s:='-0.41166214277514017'; 10 : s:='0.54444643484602580'; 11 : s:='1.0000112294976891'; 12 : s:='0.53632275170616154'; 13 : s:='-0.41916081451069583'; 14 : s:='-0.98379126325706513'; 15 : s:='-0.60774976147931451' end; Y1[3].b:=right_read(s); if k=0 then Y1[4]:=Y1[1] else begin case k of 1 : s:='0.53988150285515458'; 2 : s:='-0.41660143322567263'; 3 : s:='-0.99006293286340553'; 4 : s:='-0.65326513798591313'; 5 : s:='0.28414161158539092'; 6 : s:='0.96030987140032718'; 7 : s:='0.75357365303492832'; 8 : s:='-0.14599482609334472'; 9 : s:='-0.91133703542743327'; 10 : s:='-0.83880264087710436'; 11 : s:='0.00489285760885398'; 12 : s:='0.84389654769810143'; 13 : s:='0.90645657422145496'; 14 : s:='0.12791943802105755'; 15 : s:='-0.82271715363215553' end; Y1[4].a:=left_read(s); case k of 1 : s:='0.53988150285575403'; 2 : s:='-0.41660143322107518'; 3 : s:='-0.99006293284207678'; 4 : s:='-0.65326513786599710'; 5 : s:='0.28414161266156274'; 6 : s:='0.96030987734477706'; 7 : s:='0.75357368060371201'; 8 : s:='-0.14599456335592174'; 9 : s:='-0.91133537902956982'; 10 : s:='-0.83879618639626863'; 11 : s:='0.00495852893611134'; 12 : s:='0.84434780486321909'; 13 : s:='0.90801564986010116'; 14 : s:='0.14445583191028884'; 15 : s:='-0.70245474862055725' end; Y1[4].b:=right_read(s) end; T2:=T1+h; case k of 0 : s:='0.99999950000004166'; 1 : s:='0.53946056487236853'; 2 : s:='-0.41705592574963295'; 3 : s:='-0.99013312159395930'; 4 : s:='-0.65288649171618417'; 5 : s:='0.28462096752918731'; 6 : s:='0.96044922065186588'; 7 : s:='0.75324487959692578'; 8 : s:='-0.14648938035438312'; 9 : s:='-0.91154225954408045'; 10 : s:='-0.83852990460930026'; 11 : s:='0.00540881835694527'; 12 : s:='0.84430776350319501'; 13 : s:='0.90632947634026246'; 14 : s:='0.13116613322659680'; 15 : s:='-0.78058544201617655' end; Y2[1].a:=left_read(s); case k of 0 : s:='0.99999950000004167'; 1 : s:='0.53946056487263114'; 2 : s:='-0.41705592574822784'; 3 : s:='-0.99013312158319907'; 4 : s:='-0.65288649162793551'; 5 : s:='0.28462096796625707'; 6 : s:='0.96044922338232763'; 7 : s:='0.75324490220813002'; 8 : s:='-0.14648925792760128'; 9 : s:='-0.91154158992830521'; 10 : s:='-0.83852427241849060'; 11 : s:='0.0054425536147291544'; 12 : s:='0.84447244314170193'; 13 : s:='0.90772214240758773'; 14 : s:='0.14030080712360624'; 15 : s:='-0.73945379716579882' end; Y2[1].b:=right_read(s); case k of 0 : s:='0.00099999983333334'; 1 : s:='0.84201086628811103'; 2 : s:='0.90888082540864577'; 3 : s:='0.14012994516089170'; 4 : s:='-0.75745576045021735'; 5 : s:='-0.95864013332000371'; 6 : s:='-0.27845519017630239'; 7 : s:='0.65774016455735711'; 8 : s:='0.98921218896423957'; 9 : s:='0.41120667863463263'; 10 : s:='-0.54486194604528397'; 11 : s:='-1.0000008141934110'; 12 : s:='-0.53584960294516966'; 13 : s:='0.42052762377691990'; 14 : s:='0.98691958256144451'; 15 : s:='0.61821589131215455' end; Y2[2].a:=left_read(s); case k of 0 : s:='0.00099999983333335'; 1 : s:='0.84201086628835837'; 2 : s:='0.90888082541085727'; 3 : s:='0.14012994517510158'; 4 : s:='-0.75745576038763901'; 5 : s:='-0.95864013280518730'; 6 : s:='-0.27845518655138578'; 7 : s:='0.65774018015216203'; 8 : s:='0.98921231490493039'; 9 : s:='0.41120761950828401'; 10 : s:='-0.54485787449132743'; 11 : s:='-0.99996974709250307'; 12 : s:='-0.53560796836762878'; 13 : s:='0.42162098400036999'; 14 : s:='0.99453931021899248'; 15 : s:='0.67949922386409007' end; Y2[2].b:=right_read(s); case k of 0 : s:='-0.00099999983333335'; 1 : s:='-0.84201086628849430'; 2 : s:='-0.90888082541086921'; 3 : s:='-0.14012994517992221'; 4 : s:='0.75745576033171060'; 5 : s:='0.95864013277175199'; 6 : s:='0.27845518535850788'; 7 : s:='-0.65774019538166029'; 8 : s:='-0.98921233204898786'; 9 : s:='-0.41120787575836684'; 10 : s:='0.54485404833418288'; 11 : s:='0.99996191230411854'; 12 : s:='0.53555669162047298'; 13 : s:='-0.42253855536959966'; 14 : s:='-0.99770118080798525'; 15 : s:='-0.68791055053526634' end; Y2[3].a:=left_read(s); case k of 0 : s:='-0.00099999983333334'; 1 : s:='-0.84201086628788946'; 2 : s:='-0.90888082540874666'; 3 : s:='-0.14012994515606490'; 4 : s:='0.75745576050638495'; 5 : s:='0.95864013335352702'; 6 : s:='0.27845519136915788'; 7 : s:='-0.65774014932774752'; 8 : s:='-0.98921217182024776'; 9 : s:='-0.41120642238386052'; 10 : s:='0.54486577225152785'; 11 : s:='1.0000086502022759'; 12 : s:='0.53590083465889381'; 13 : s:='-0.41961286339374584'; 14 : s:='-0.98385611413184558'; 15 : s:='-0.60731398724691633' end; Y2[3].b:=right_read(s); if k=0 then Y2[4]:=Y2[1] else begin case k of 1 : s:='0.53946056487212131'; 2 : s:='-0.41705592575143449'; 3 : s:='-0.99013312159971258'; 4 : s:='-0.65288649173236853'; 5 : s:='0.28462096720899887'; 6 : s:='0.96044921904319494'; 7 : s:='0.75324487710071645'; 8 : s:='-0.14648945064701035'; 9 : s:='-0.91154275347888794'; 10 : s:='-0.83853031960067200'; 11 : s:='0.00539281972969396'; 12 : s:='0.84416434942739940'; 13 : s:='0.90624533161355768'; 14 : s:='0.12741556290556430'; 15 : s:='-0.82309667688728127' end; Y2[4].a:=left_read(s); case k of 1 : s:='0.53946056487272125'; 2 : s:='-0.41705592574683262'; 3 : s:='-0.99013312157837445'; 4 : s:='-0.65288649161230267'; 5 : s:='0.28462096828631441'; 6 : s:='0.96044922499097135'; 7 : s:='0.75324490470419362'; 8 : s:='-0.14648918763530325'; 9 : s:='-0.91154109599553146'; 10 : s:='-0.83852385741683526'; 11 : s:='0.00545855341307715'; 12 : s:='0.84461594311595907'; 13 : s:='0.90780604122061913'; 14 : s:='0.14396882346389341'; 15 : s:='-0.70273372519990820' end; Y2[4].b:=right_read(s) end; n:=2; max_l:=0; repeat n:=n+1; T3:=T2+h; T:=T2+hh; for i:=1 to 4 do Y[i]:=Y2[i]+hh*Fdelta[i]; i1:=Psi3Explicit(T,Y); for i:=1 to 4 do begin i1[i]:=9*i1[i]-i1[i]; i2[i]:=(h*h)*(h*h); i1[i]:=i2[i]*i1[i]/24 end; i2:=FTY(T2,Y2); i3:=FTY(T1,Y1); i4:=FTY(T0,Y0); for i:=1 to 4 do begin i2[i]:=7*i2[i]-2*i3[i]; i2[i]:=i2[i]+i4[i]; i2[i]:=h*i2[i]/3; Y3[i]:=Y1[i]+(i2[i]+i1[i]) end; T:=T3+hh1; // iterations l:=0; // number of iterations repeat l:=l+1; find:=True; for i:=1 to 4 do Y[i]:=Y3[i]+hh1*Fdelta[i]; i1:=Psi3Implicit(T,Y); for i:=1 to 4 do begin i1[i]:=11*i1[i]-19*i1[i]; i2[i]:=(h*h)*(h*h); i2[i]:=i2[i]*h; i1[i]:=i2[i]*i1[i]/720 end; i2:=FTY(T3,Y3); i3:=FTY(T2,Y2); for i:=1 to 4 do i2[i]:=i2[i]+4*i3[i]; i3:=FTY(T1,Y1); for i:=1 to 4 do begin i2[i]:=i2[i]+i3[i]; i2[i]:=h*i2[i]/3; Y4[i]:=Y1[i]+(i2[i]+i1[i]) end; for i:=1 to 4 do if (Y3[i].a=0) or (Y3[i].b=0) then if (Abs(Y4[i].a-Y3[i].a)>=1e-18) or (Abs(Y4[i].b-Y3[i].b)>=1e-18) then find:=False // 1e-18 - accuracy else else if (Abs((Y4[i].a-Y3[i].a)/Y3[i].a)>=1e-18) or (Abs((Y4[i].b-Y3[i].b)/Y3[i].b)>=1e-18) then find:=False; Y3:=Y4; until (l=20) or find; // 20 - max number of iterations if not find then Writeln ('For n = ', n, ' accuracy not achieved after 20 ' +'iterations') else if n=2000 then Results (T3, Y4, 'i'); T0:=T1; Y0:=Y1; T1:=T2; Y1:=Y2; T2:=T3; Y2:=Y3; if l>max_l then max_l:=l until (n=m) or not find; k:=k+1 until k=16; if z='f' then CloseFile (plik); Writeln; Writeln ('The end. Press Enter ...'); Readln end.