program IPCM_Example_7_WrappingEffect; {$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 four-step explicit interval // method of Nystrom. It is assumed h = 0.0005 and 0 <= t <= 5. At t = 5 the // wrapping effect is self-evident. // 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 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 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, i5, 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: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+']'); 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'); if z='f' then begin AssignFile (plik, 'IPCM7WERESULTS.TXT'); Rewrite (plik); Writeln ('Output to IPCM7WERESULTS.TXT in current directory') end else Writeln; if z='s' then begin Writeln ('Exact solution:'); Writeln ('y[1](1.0) = ', Cos(1.0):26, ' y[2](1.0) = ', Sin(1.0):26); Writeln ('y[3](1.0) = ', -Sin(1.0):26, ' y[4](1.0) = ', Cos(1.0):26); Writeln ('y[1](3.0) = ', Cos(3.0):26, ' y[2](3.0) = ', Sin(3.0):26); Writeln ('y[3](3.0) = ', -Sin(3.0):26, ' y[4](3.0) = ', Cos(3.0):26); Writeln ('y[1](5.0) = ', Cos(5.0):26, ' y[2](5.0) = ', Sin(5.0):26); Writeln ('y[3](5.0) = ', -Sin(5.0):26, ' y[4](5.0) = ', Cos(5.0):26); Writeln; Readln end else begin Writeln (plik, 'Exact solution:'); Writeln (plik, '==============='); Writeln (plik, 'y[1](1.0) = ', Cos(1.0):26, ' y[2](1.0) = ', Sin(1.0):26); Writeln (plik, 'y[3](1.0) = ', -Sin(1.0):26, ' y[4](1.0) = ', Cos(1.0):26); Writeln (plik, 'y[1](3.0) = ', Cos(3.0):26, ' y[2](3.0) = ', Sin(3.0):26); Writeln (plik, 'y[3](3.0) = ', -Sin(3.0):26, ' y[4](3.0) = ', Cos(3.0):26); Writeln (plik, 'y[1](5.0) = ', Cos(5.0):26, ' y[2](5.0) = ', Sin(5.0):26); Writeln (plik, 'y[3](5.0) = ', -Sin(5.0):26, ' y[4](5.0) = ', Cos(5.0):26) end; Writeln ('Be patient ... Calculations may take a lot of minutes.'); Writeln; if z='f' then begin Writeln (plik, ' '); Writeln (plik, 'k = 4, Nystrom (explicit)'); Writeln (plik, '=========================') end else begin Writeln ('k = 4, Nystrom (explicit)'); Writeln end; m:=2000; s:='0.0005'; h.a:=left_read(s); h.b:=right_read(s); 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 end; case k of 0 : s:='1.001'; 1 : s:='2.001'; 2 : s:='3.001'; 3 : s:='4.001'; 4 : s:='5.001' end; delta_t.b:=right_read(s); if k=3 then delta_y[1].a:=-1 else begin case k of 0 : s:='0.4'; 1 : s:='-0.5'; 2 : s:='-1.1'; 4 : s:='-0.6' end; delta_y[1].a:=left_read(s) end; case k of 0 : s:='1.1'; 1 : s:='0.6'; 2 : s:='-0.3'; 3 : s:='-0.6'; 4 : s:='0.4' end; delta_y[1].b:=right_read(s); delta_y[4]:=delta_y[1]; if k=4 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' end; delta_y[2].a:=left_read(s) end; if (k=0) or (k=2) then delta_y[2].b:=1 else begin case k of 1 : s:='1.1'; 3 : s:='0.2'; 4 : s:='-0.7' end; delta_y[2].b:=right_read(s) end; if k<=2 then delta_y[3].a:=-1 else begin case k of 3 : s:='-0.2'; 4 : s:='0.7' end; delta_y[3].a:=left_read(s) end; if k=4 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' 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 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.54030230583219599'; 2 : s:='-0.41614684036671793'; 3 : s:='-0.98999635060230819'; 4 : s:='-0.65503034304622548' end; Y0[1].a:=left_read(s); case k of 1 : s:='0.54030230590115955'; 2 : s:='-0.41614683273255191'; 3 : s:='-0.98998864274042803'; 4 : s:='-0.65225665909786173' end; Y0[1].b:=right_read(s); case k of 1 : s:='0.84147098477955367'; 2 : s:='0.90929741250467681'; 3 : s:='0.14111737625646685'; 4 : s:='-0.75607167693982105' end; Y0[2].a:=left_read(s); case k of 1 : s:='0.84147098483777535'; 2 : s:='0.90929744115396114'; 3 : s:='0.14112263994600178'; 4 : s:='-0.75607167693982105' end; Y0[2].b:=right_read(s); case k of 1 : s:='-0.84147098487205889'; 2 : s:='-0.90929743529475427'; 3 : s:='-0.14112806615015036'; 4 : s:='0.75431354812173348' end; Y0[3].a:=left_read(s); case k of 1 : s:='-0.84147098474731947'; 2 : s:='-0.90929741835590930'; 3 : s:='-0.14111194976199424'; 4 : s:='0.75929412518838869' end; Y0[3].b:=right_read(s); case k of 1 : s:='0.54030230580557908'; 2 : s:='-0.41614686480058189'; 3 : s:='-0.98999623632066398'; 4 : s:='-0.65538913295807498' end; Y0[4].a:=left_read(s); case k of 1 : s:='0.54030230593241650'; 2 : s:='-0.41614680828925584'; 3 : s:='-0.98998875691905698'; 4 : s:='-0.65189605434989212' end; Y0[4].b:=right_read(s) end; T1:=T0+h; case k of 0 : s:='0.99999987500000260'; 1 : s:='0.53988150281940306'; 2 : s:='-0.41660143705524524'; 3 : s:='-0.99006679879541348'; 4 : s:='-0.65465558402265320' end; Y1[1].a:=left_read(s); case k of 0 : s:='0.99999987500000261'; 1 : s:='0.53988150288855708'; 2 : s:='-0.41660142939569300'; 3 : s:='-0.99005906678717803'; 4 : s:='-0.65187443661413741' end; Y1[1].b:=right_read(s); case k of 0 : s:='0.00049999997916666'; 1 : s:='0.84174103073723529'; 2 : s:='0.90908922538794385'; 3 : s:='0.14062235671701562'; 4 : s:='-0.75786249474417313' end; Y1[2].a:=left_read(s); case k of 0 : s:='0.00049999997916667'; 1 : s:='0.84174103079565032'; 2 : s:='0.90908925412192514'; 3 : s:='0.14062763161811609'; 4 : s:='-0.75639578308343477' end; Y1[2].b:=right_read(s); case k of 0 : s:='-0.00049999997916667'; 1 : s:='-0.84174103082997911'; 2 : s:='-0.90908924825725691'; 3 : s:='-0.14063307688722787'; 4 : s:='0.75463390682909712' end; Y1[3].a:=left_read(s); case k of 0 : s:='-0.00049999997916666'; 1 : s:='-0.84174103070493290'; 2 : s:='-0.90908923125235061'; 3 : s:='-0.14061691168917846'; 4 : s:='0.75962721878187472' end; Y1[3].b:=right_read(s); if k=0 then Y1[4]:=Y1[1] else begin case k of 1 : s:='0.53988150279276498'; 2 : s:='-0.41660146155445341'; 3 : s:='-0.99006667913196906'; 4 : s:='-0.65501671067704284' end; Y1[4].a:=left_read(s); case k of 1 : s:='0.53988150292001072'; 2 : s:='-0.41660140488078198'; 3 : s:='-0.99005918652260017'; 4 : s:='-0.65151151420019268' end; Y1[4].b:=right_read(s) end; T2:=T1+h; case k of 0 : s:='0.99999950000004166'; 1 : s:='0.53946056483631381'; 2 : s:='-0.41705592959401476'; 3 : s:='-0.99013699977373680'; 4 : s:='-0.65428068552896519' end; Y2[1].a:=left_read(s); case k of 0 : s:='0.99999950000004167'; 1 : s:='0.53946056490565187'; 2 : s:='-0.41705592190900989'; 3 : s:='-0.99012924354625621'; 4 : s:='-0.65149205556896229' end; Y2[1].b:=right_read(s); case k of 0 : s:='0.00099999983333334'; 1 : s:='0.84201086625972458'; 2 : s:='0.90888081100402671'; 3 : s:='0.14012730214362568'; 4 : s:='-0.75819166177775069' end; Y2[2].a:=left_read(s); case k of 0 : s:='0.00099999983333335'; 1 : s:='0.84201086631832717'; 2 : s:='0.90888083982289961'; 3 : s:='0.14013258827578291'; 4 : s:='-0.75671969937217923' end; Y2[2].b:=right_read(s); case k of 0 : s:='-0.00099999983333335'; 1 : s:='-0.84201086635273221'; 2 : s:='-0.90888083394510055'; 3 : s:='-0.14013805213967895'; 4 : s:='0.75495407743972587' end; Y2[3].a:=left_read(s); case k of 0 : s:='-0.00099999983333334'; 1 : s:='-0.84201086622736447'; 2 : s:='-0.90888081687382780'; 3 : s:='-0.14012183798814049'; 4 : s:='0.75996015478892618' end; Y2[3].b:=right_read(s); if k=0 then Y2[4]:=Y2[1] else begin case k of 1 : s:='0.53946056480946804'; 2 : s:='-0.41705595416507771'; 3 : s:='-0.99013687452353425'; 4 : s:='-0.65464412892234860' end; Y2[4].a:=left_read(s); case k of 1 : s:='0.53946056493713680'; 2 : s:='-0.41705589732851296'; 3 : s:='-0.99012936869293396'; 4 : s:='-0.65112677432064900' end; Y2[4].b:=right_read(s) end; T3:=T2+h; case k of 0 : s:='0.99999887500021093'; 1 : s:='0.53903949198800898'; 2 : s:='-0.41751031786830209'; 3 : s:='-0.99020695299047841'; 4 : s:='-0.65390561904318064' end; Y3[1].a:=left_read(s); case k of 0 : s:='0.99999887500021094'; 1 : s:='0.53903949205753844'; 2 : s:='-0.41751031015771265'; 3 : s:='-0.99019917247016361'; 4 : s:='-0.65110948740809495' end; Y3[1].b:=right_read(s); case k of 0 : s:='0.00149999943750006'; 1 : s:='0.84228049127943993'; 2 : s:='0.90867216939453299'; 3 : s:='0.13963221239549290'; 4 : s:='-0.75852064000500906' end; Y3[2].a:=left_read(s); case k of 0 : s:='0.00149999943750007'; 1 : s:='0.84228049133823712'; 2 : s:='0.90867219829859024'; 3 : s:='0.13963750977884418'; 4 : s:='-0.75704340865518897' end; Y3[2].b:=right_read(s); case k of 0 : s:='-0.00149999943750007'; 1 : s:='-0.84228049137268754'; 2 : s:='-0.90867219241536222'; 3 : s:='-0.13964299283213669'; 4 : s:='0.75527402702264953' end; Y3[3].a:=left_read(s); case k of 0 : s:='-0.00149999943750006'; 1 : s:='-0.84228049124701158'; 2 : s:='-0.90867217527751340'; 3 : s:='-0.13962672958418007'; 4 : s:='0.76029290031515090' end; Y3[3].b:=right_read(s); if k=0 then Y3[4]:=Y3[1] else begin case k of 1 : s:='0.53903949196114155'; 2 : s:='-0.41751034250504894'; 3 : s:='-0.99020682229883774'; 4 : s:='-0.65427142038106588' end; Y3[4].a:=left_read(s); case k of 1 : s:='0.53903949208922128'; 2 : s:='-0.41751028550524352'; 3 : s:='-0.99019930323451345'; 4 : s:='-0.65074186727943624' end; Y3[4].b:=right_read(s) end; Fdelta:=FTY(delta_t,delta_y); s:='-0.0015'; hh.a:=left_read(s); s:='0.0005'; hh.b:=right_read(s); for n:=4 to m do begin T4:=T3+h; T:=T3+hh; for i:=1 to 4 do Y[i]:=Y3[i]+hh*Fdelta[i]; i1:=Psi4Explicit(T,Y); for i:=1 to 4 do begin i1[i]:=251*i1[i]-19*i1[i]; i2[i]:=(h*h)*(h*h); i2[i]:=i2[i]*h*i1[i]; i1[i]:=i2[i]/720 end; i2:=FTY(T3,Y3); i3:=FTY(T2,Y2); i4:=FTY(T1,Y1); i5:=FTY(T0,Y0); for i:=1 to 4 do begin i2[i]:=8*i2[i]-5*i3[i]; i2[i]:=i2[i]+4*i4[i]; i2[i]:=i2[i]-i5[i]; i2[i]:=h*i2[i]/3; Y4[i]:=Y2[i]+(i2[i]+i1[i]) end; if n=m then Results (T4, Y4, 'e'); T0:=T1; Y0:=Y1; T1:=T2; Y1:=Y2; T2:=T3; Y2:=Y3; T3:=T4; Y3:=Y4 end; k:=k+1 until k=5; if z='f' then CloseFile (plik); Writeln; Writeln ('The end. Press Enter ...'); Readln end.