!***************************************************************************************** !> ! Dense output coefficients for the variable-step RK formulas. ! ! Each `*_dense_coefficients` routine returns the matrix `b` (`number_of_stages` x `degree`) ! defining the interpolant over the step \( [t, t+h] \), along with its form: either ! the polynomial form (`dense_form_polynomial`), or the form used in `DOPRI5` and `DOP853` ! (`dense_form_dopri`). See [[dense_coefficients]]. So, the coefficients can be entered ! as they are published. ! ! If the dense output needs function evaluations in addition to the ones from the step, ! these are computed in the `*_dense_extra_stages` routine, and stored in `me%funcs` ! after the ones from the step. submodule(rklib_module) rklib_variable_dense implicit none contains !***************************************************************************************** !***************************************************************************************** !> ! Dense output for the Dormand-Prince 5(4) method [[rkdp54]]. ! This is the 4th order continuous extension due to Shampine (as used in `DOPRI5`). ! The coefficients are given here in the polynomial form \( b_i(\theta) \) ! (row `i` contains the coefficients of \( \theta^1 ... \theta^4 \)). ! !### References ! * L. F. Shampine, "Some Practical Runge-Kutta Formulas", ! Mathematics of Computation, 46(173), 1986, p. 135-150. ! * E. Hairer, S. P. Norsett, G. Wanner, "Solving Ordinary Differential Equations I: ! Nonstiff Problems", 2nd ed., Springer, 1993, p. 191-192. module procedure rkdp54_dense_coefficients form = dense_form_polynomial allocate(b(7,4)) b(1,:) = [ 1.0_wp, & -8048581381.0_wp / 2820520608.0_wp, & 8663915743.0_wp / 2820520608.0_wp, & -12715105075.0_wp / 11282082432.0_wp ] b(2,:) = 0.0_wp b(3,:) = [ 0.0_wp, & 131558114200.0_wp / 32700410799.0_wp, & -68118460800.0_wp / 10900136933.0_wp, & 87487479700.0_wp / 32700410799.0_wp ] b(4,:) = [ 0.0_wp, & -1754552775.0_wp / 470086768.0_wp, & 14199869525.0_wp / 1410260304.0_wp, & -10690763975.0_wp / 1880347072.0_wp ] b(5,:) = [ 0.0_wp, & 127303824393.0_wp / 49829197408.0_wp, & -318862633887.0_wp / 49829197408.0_wp, & 701980252875.0_wp / 199316789632.0_wp ] b(6,:) = [ 0.0_wp, & -282668133.0_wp / 205662961.0_wp, & 2019193451.0_wp / 616988883.0_wp, & -1453857185.0_wp / 822651844.0_wp ] b(7,:) = [ 0.0_wp, & 40617522.0_wp / 29380423.0_wp, & -110615467.0_wp / 29380423.0_wp, & 69997945.0_wp / 29380423.0_wp ] end procedure rkdp54_dense_coefficients !***************************************************************************************** !***************************************************************************************** !> ! Dense output for the Dormand-Prince 8(5) method [[rkdp85]]. ! This is the 7th order continuous extension from `DOP853`. It requires 4 additional ! function evaluations (see [[rkdp85_dense_extra_stages]]). ! ! The coefficients are given here in the form used in `DOP853`: ! ! \( x(t + \theta h) = r_0 + \theta(r_1 + \theta_1(r_2 + \theta(r_3 + \theta_1(r_4 + ! \theta(r_5 + \theta_1(r_6 + \theta r_7)))))) \), where \( \theta_1 = 1 - \theta \) ! ! where \( r_0 = x(t) \) and \( r_j = h \sum_i b_{ij} f_i \). ! !### References ! * E. Hairer, S. P. Norsett, G. Wanner, "Solving Ordinary Differential Equations I: ! Nonstiff Problems", 2nd ed., Springer, 1993. ! * [DOP853](https://github.com/jacobwilliams/dop853) module procedure rkdp85_dense_coefficients ! weights of the 8th order solution (see [[rkdp85]]): real(wp),dimension(12),parameter :: w = [ 5.42937341165687622380535766363e-2_wp, & 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, & 4.45031289275240888144113950566_wp, & 1.89151789931450038304281599044_wp, & -5.8012039600105847814672114227_wp, & 3.1116436695781989440891606237e-1_wp, & -1.52160949662516078556178806805e-1_wp, & 2.01365400804030348374776537501e-1_wp, & 4.47106157277725905176885569043e-2_wp ] form = dense_form_dopri allocate(b(16,7)) b = 0.0_wp ! r1 = x(t+h) - x(t) b(1:12,1) = w ! r2 = h*f1 - r1 b(1:12,2) = -w b(1,2) = b(1,2) + 1.0_wp ! r3 = r1 - h*f13 - r2 b(1:12,3) = 2.0_wp*w b(1,3) = b(1,3) - 1.0_wp b(13,3) = -1.0_wp ! r4 - r7: b([1,6,7,8,9,10,11,12,13,14,15,16],4) = [ & -0.84289382761090128651353491142e+01_wp, 0.56671495351937776962531783590e+00_wp, & -0.30689499459498916912797304727e+01_wp, 0.23846676565120698287728149680e+01_wp, & 0.21170345824450282767155149946e+01_wp, -0.87139158377797299206789907490e+00_wp, & 0.22404374302607882758541771650e+01_wp, 0.63157877876946881815570249290e+00_wp, & -0.88990336451333310820698117400e-01_wp, 0.18148505520854727256656404962e+02_wp, & -0.91946323924783554000451984436e+01_wp, -0.44360363875948939664310572000e+01_wp ] b([1,6,7,8,9,10,11,12,13,14,15,16],5) = [ & 0.10427508642579134603413151009e+02_wp, 0.24228349177525818288430175319e+03_wp, & 0.16520045171727028198505394887e+03_wp, -0.37454675472269020279518312152e+03_wp, & -0.22113666853125306036270938578e+02_wp, 0.77334326684722638389603898808e+01_wp, & -0.30674084731089398182061213626e+02_wp, -0.93321305264302278729567221706e+01_wp, & 0.15697238121770843886131091075e+02_wp, -0.31139403219565177677282850411e+02_wp, & -0.93529243588444783865713862664e+01_wp, 0.35816841486394083752465898540e+02_wp ] b([1,6,7,8,9,10,11,12,13,14,15,16],6) = [ & 0.19985053242002433820987653617e+02_wp, -0.38703730874935176555105901742e+03_wp, & -0.18917813819516756882830838328e+03_wp, 0.52780815920542364900561016686e+03_wp, & -0.11573902539959630126141871134e+02_wp, 0.68812326946963000169666922661e+01_wp, & -0.10006050966910838403183860980e+01_wp, 0.77771377980534432092869265740e+00_wp, & -0.27782057523535084065932004339e+01_wp, -0.60196695231264120758267380846e+02_wp, & 0.84320405506677161018159903784e+02_wp, 0.11992291136182789328035130030e+02_wp ] b([1,6,7,8,9,10,11,12,13,14,15,16],7) = [ & -0.25693933462703749003312586129e+02_wp, -0.15418974869023643374053993627e+03_wp, & -0.23152937917604549567536039109e+03_wp, 0.35763911791061412378285349910e+03_wp, & 0.93405324183624310003907691704e+02_wp, -0.37458323136451633156875139351e+02_wp, & 0.10409964950896230045147246184e+03_wp, 0.29840293426660503123344363579e+02_wp, & -0.43533456590011143754432175058e+02_wp, 0.96324553959188282948394950600e+02_wp, & -0.39177261675615439165231486172e+02_wp, -0.14972683625798562581422125276e+03_wp ] end procedure rkdp85_dense_coefficients !***************************************************************************************** !***************************************************************************************** !> ! The additional function evaluations for the dense output of [[rkdp85]]. ! These are stored in `me%funcs(:,13:16)`. Stage 13 is `f` at the end of the step, ! which is cached for the first stage of the next step (as in `DOP853`). module procedure rkdp85_dense_extra_stages real(wp),parameter :: c14 = 0.1_wp real(wp),parameter :: c15 = 0.2_wp real(wp),parameter :: c16 = 7.0_wp / 9.0_wp real(wp),parameter :: a141 = 5.61675022830479523392909219681e-2_wp real(wp),parameter :: a147 = 2.53500210216624811088794765333e-1_wp real(wp),parameter :: a148 = -2.46239037470802489917441475441e-1_wp real(wp),parameter :: a149 = -1.24191423263816360469010140626e-1_wp real(wp),parameter :: a1410 = 1.5329179827876569731206322685e-1_wp real(wp),parameter :: a1411 = 8.20105229563468988491666602057e-3_wp real(wp),parameter :: a1412 = 7.56789766054569976138603589584e-3_wp real(wp),parameter :: a1413 = -8.298e-3_wp real(wp),parameter :: a151 = 3.18346481635021405060768473261e-2_wp real(wp),parameter :: a156 = 2.83009096723667755288322961402e-2_wp real(wp),parameter :: a157 = 5.35419883074385676223797384372e-2_wp real(wp),parameter :: a158 = -5.49237485713909884646569340306e-2_wp real(wp),parameter :: a1511 = -1.08347328697249322858509316994e-4_wp real(wp),parameter :: a1512 = 3.82571090835658412954920192323e-4_wp real(wp),parameter :: a1513 = -3.40465008687404560802977114492e-4_wp real(wp),parameter :: a1514 = 1.41312443674632500278074618366e-1_wp real(wp),parameter :: a161 = -4.28896301583791923408573538692e-1_wp real(wp),parameter :: a166 = -4.69762141536116384314449447206e0_wp real(wp),parameter :: a167 = 7.68342119606259904184240953878e0_wp real(wp),parameter :: a168 = 4.06898981839711007970213554331e0_wp real(wp),parameter :: a169 = 3.56727187455281109270669543021e-1_wp real(wp),parameter :: a1613 = -1.39902416515901462129418009734e-3_wp real(wp),parameter :: a1614 = 2.9475147891527723389556272149e0_wp real(wp),parameter :: a1615 = -9.15095847217987001081870187138e0_wp associate (t => me%dense_t, & h => me%dense_h, & x => me%dense_x, & f1 => me%funcs(:,1), & f6 => me%funcs(:,6), & f7 => me%funcs(:,7), & f8 => me%funcs(:,8), & f9 => me%funcs(:,9), & f10 => me%funcs(:,10), & f11 => me%funcs(:,11), & f12 => me%funcs(:,12), & f13 => me%funcs(:,13), & f14 => me%funcs(:,14), & f15 => me%funcs(:,15), & f16 => me%funcs(:,16)) ! f at the end of the step (this is also the first stage of the next step, ! so it is saved in the function evaluation cache): call me%set_fsal_cache(t+h, me%dense_xf, f13) if (me%stopped) return call me%f(t+c14*h, x+h*(a141*f1+a147*f7+a148*f8+a149*f9+a1410*f10+& a1411*f11+a1412*f12+a1413*f13), f14) if (me%stopped) return call me%f(t+c15*h, x+h*(a151*f1+a156*f6+a157*f7+a158*f8+a1511*f11+& a1512*f12+a1513*f13+a1514*f14), f15) if (me%stopped) return call me%f(t+c16*h, x+h*(a161*f1+a166*f6+a167*f7+a168*f8+a169*f9+& a1613*f13+a1614*f14+a1615*f15), f16) end associate end procedure rkdp85_dense_extra_stages !***************************************************************************************** ! BEGIN GENERATED DENSE OUTPUT (generated by scripts/generate_dense.py -- do not edit by hand) !***************************************************************************************** !> ! Dense output for Verner's "most efficient" 6(5) method [[rkv65e]]. ! This is the interpolant of order 6 (`bi6`) from the coefficient file, which requires 3 additional ! function evaluations (see [[rkv65e_dense_extra_stages]]). ! !@note The interpolant is published in the polynomial form. Here it has been ! converted (exactly, using high precision arithmetic) ! to the form used in `DOPRI5` and `DOP853` (see [[dense_coefficients]]), which has ! much smaller roundoff error for these coefficients. This code is generated by ! `scripts/generate_dense.py`. ! !### References ! * J.H. Verner, "Numerically optimal Runge--Kutta pairs with interpolants", ! Numerical Algorithms, 53, (2010) pp. 383-396. [doi](https://doi.org/10.1007/s11075-009-9290-3) ! * Coefficients: ! * [RKV65.IIIXb.Efficient.00000144617.081204.CoeffsOnlyFLOAT](https://web.archive.org/web/https://www.sfu.ca/~jverner/RKV65.IIIXb.Efficient.00000144617.081204.CoeffsOnlyFLOAT) module procedure rkv65e_dense_coefficients form = dense_form_dopri allocate(b(12,6)) b = 0.0_wp b(1,:) = [ 3.438957868357036009278820124728322385900e-2_wp, & 9.656104213164296399072117987527167761410e-1_wp, & -9.312208426328592798144235975054335522500e-1_wp, & -5.881762592546286216504782819308229424370e0_wp, & 3.174650111594311357773083170065202627230e0_wp, & 1.118304243294735597460358736881373385894e1_wp ] b(4,:) = [ 2.582624555633503404659558098586120858100e-1_wp, & -2.582624555633503404659558098586120858100e-1_wp, & 5.165249111267006809319116197172241719100e-1_wp, & 1.585731477158971090460968672531878207274e1_wp, & -6.198298933520408171182939436606690061060e0_wp, & -4.782638065987969267888070552937260849568e1_wp ] b(5,:) = [ 4.209371189673537150642551514069801966100e-1_wp, & -4.209371189673537150642551514069801966100e-1_wp, & 8.418742379347074301285103028139603936200e-1_wp, & 2.584553910459551810494526629638858407746e1_wp, & -1.010249085521648916154212363376752472094e1_wp, & -7.795131832728772501189910211240374013022e1_wp ] b(6,:) = [ 4.405396469669310170148836816197095665100e0_wp, & -4.405396469669310170148836816197095665100e0_wp, & 8.810792939338620340297673632394191329600e0_wp, & 2.704913432376956444471385805145016738243e2_wp, & -1.057295152720634440835720835887302959577e2_wp, & -8.158141610498722537312660770735362342391e2_wp ] b(7,:) = [ -1.764831190242986576151740942499002125500e2_wp, & 1.764831190242986576151740942499002125500e2_wp, & -3.529662380485973152303481884998004248800e2_wp, & -1.083606350809193757757168938694387304775e4_wp, & 4.235594856583167782764178261997605100050e3_wp, & 3.268205907857382548429149893516670601905e4_wp ] b(8,:) = [ 1.723641334014150730294022582711902413500e2_wp, & -1.723641334014150730294022582711902413500e2_wp, & 3.447282668028301460588045165423804826100e2_wp, & 1.058315779084688548400529865785108081779e4_wp, & -4.136739201633961752705654198508565791940e3_wp, & -3.191928396322501352396338116133152617251e4_wp ] b(9,:) = [ -6.000000000000000000000000000000000000000e-39_wp, & 6.000000000000000000000000000000000000000e-39_wp, & -9.999999999999999999999999999999999999880e-1_wp, & -2.002282515073212747631352282515073212754e0_wp, & 3.999999999999999999999999999999999999996e0_wp, & 5.941557405812358439404089705553960506588e0_wp ] b(10,:) = [ -3.000000000000000000000000000000000000000e-38_wp, & 3.000000000000000000000000000000000000000e-38_wp, & 3.000000000000000000000000000000000000000e-38_wp, & -1.961034737620103473762010347376201034738e1_wp, & 1.600000000000000000000000000000000000002e1_wp, & 3.791300539268018942815690783170457967206e1_wp ] b(11,:) = [ 4.000000000000000000000000000000000000000e-38_wp, & -4.000000000000000000000000000000000000000e-38_wp, & -1.200000000000000000000000000000000000000e-37_wp, & -1.823029074639408932621452105403485845596e1_wp, & 4.000000000000000000000000000000000000000e-38_wp, & 6.028535299733495147557711988768141023822e1_wp ] b(12,:) = [ -5.000000000000000000000000000000000000000e-38_wp, & 5.000000000000000000000000000000000000000e-38_wp, & 1.400000000000000000000000000000000000000e-37_wp, & -1.356379663861415686233204481379663861425e1_wp, & -5.000000000000000000000000000000000000000e-38_wp, & 6.349378645945285577618640608637905204545e1_wp ] end procedure rkv65e_dense_coefficients !***************************************************************************************** !***************************************************************************************** !> ! The additional function evaluations for the dense output of [[rkv65e]]. ! These are stored in `me%funcs(:,10:12)`. module procedure rkv65e_dense_extra_stages real(wp),parameter :: c10 = 0.5_wp real(wp),parameter :: a1001 = .1652415901357280684383619367363197852645e-1_wp real(wp),parameter :: a1004 = .3053128187514178931377105638345583032476_wp real(wp),parameter :: a1005 = .2071200938201978848991082158995582390341_wp real(wp),parameter :: a1006 = -1.293879140655123187129665774327355723229_wp real(wp),parameter :: a1007 = 57.11988411588149149650780257779402737914_wp real(wp),parameter :: a1008 = -55.87979207510932290773937033203265749155_wp real(wp),parameter :: a1009 = .2483002829776601348057855515823731483430e-1_wp real(wp),parameter :: c11 = 0.828_wp real(wp),parameter :: a1101 = .3815008181862774566390122945290392930051e-1_wp real(wp),parameter :: a1104 = .2502358252513705267411717435114864276070_wp real(wp),parameter :: a1105 = .3249441447817607965203943927614998837766_wp real(wp),parameter :: a1106 = 1.822460665832796218309507694895413940004_wp real(wp),parameter :: a1107 = -67.71372332692620404266416384783466979895_wp real(wp),parameter :: a1108 = 66.03587911808127345206414399075146898997_wp real(wp),parameter :: a1109 = -.3638810874951269663495520353810337170236e-1_wp real(wp),parameter :: a1110 = .106441599909888_wp real(wp),parameter :: c12 = 0.28_wp real(wp),parameter :: a1201 = .1117816803966601180374383193782622585818_wp real(wp),parameter :: a1204 = .2575750510934521211743051533509279738017e-1_wp real(wp),parameter :: a1205 = 3.785140856363645858969575735188195099293_wp real(wp),parameter :: a1206 = 92.34088993695727177430761424880372341560_wp real(wp),parameter :: a1207 = -3819.461508432344125782385301764384773807_wp real(wp),parameter :: a1208 = 3732.492711530704227040059052662831877264_wp real(wp),parameter :: a1209 = -1.075694020996303327721922083840048543815_wp real(wp),parameter :: a1210 = -3.231539970732086111815287724591251178493_wp real(wp),parameter :: a1211 = -4.707539085458634781568599908721077305888_wp associate (t => me%dense_t, h => me%dense_h, x => me%dense_x, k => me%funcs) call me%f(t+c10*h, x+h*(a1001*k(:,1) + a1004*k(:,4) + a1005*k(:,5) + a1006*k(:,6) + & a1007*k(:,7) + a1008*k(:,8) + a1009*k(:,9)), k(:,10)) if (me%stopped) return call me%f(t+c11*h, x+h*(a1101*k(:,1) + a1104*k(:,4) + a1105*k(:,5) + a1106*k(:,6) + & a1107*k(:,7) + a1108*k(:,8) + a1109*k(:,9) + a1110*k(:,10)), k(:,11)) if (me%stopped) return call me%f(t+c12*h, x+h*(a1201*k(:,1) + a1204*k(:,4) + a1205*k(:,5) + a1206*k(:,6) + & a1207*k(:,7) + a1208*k(:,8) + a1209*k(:,9) + a1210*k(:,10) + & a1211*k(:,11)), k(:,12)) end associate end procedure rkv65e_dense_extra_stages !***************************************************************************************** !***************************************************************************************** !> ! Dense output for Verner's "most robust" 6(5) method [[rkv65r]]. ! This is the interpolant of order 6 (`bi6`) from the coefficient file, which requires 3 additional ! function evaluations (see [[rkv65r_dense_extra_stages]]). ! !@note The interpolant is published in the polynomial form. Here it has been ! converted (exactly, using high precision arithmetic) ! to the form used in `DOPRI5` and `DOP853` (see [[dense_coefficients]]), which has ! much smaller roundoff error for these coefficients. This code is generated by ! `scripts/generate_dense.py`. ! !### References ! * J.H. Verner, "Numerically optimal Runge--Kutta pairs with interpolants", ! Numerical Algorithms, 53, (2010) pp. 383-396. [doi](https://doi.org/10.1007/s11075-009-9290-3) ! * Coefficients: ! * [RKV65.IIIXb.Robust.00010102836.081204.CoeffsOnlyFLOAT](https://web.archive.org/web/https://www.sfu.ca/~jverner/RKV65.IIIXb.Robust.00010102836.081204.CoeffsOnlyFLOAT) module procedure rkv65r_dense_coefficients form = dense_form_dopri allocate(b(12,6)) b = 0.0_wp b(1,:) = [ 7.638888888888888888888888888888888888900e-2_wp, & 9.236111111111111111111111111111111111110e-1_wp, & -8.472222222222222222222222222222222222230e-1_wp, & -2.749606387764282501124606387764282501124e0_wp, & 2.166666666666666666666666666666666666667e0_wp, & -1.105862947968211126105862947968211126106e0_wp ] b(4,:) = [ 3.694083694083694083694083694083694083300e-1_wp, & -3.694083694083694083694083694083694083300e-1_wp, & 7.388167388167388167388167388167388168600e-1_wp, & 1.346919746919746919746919746919746919739e1_wp, & -8.865800865800865800865800865800865800910e0_wp, & -5.683205683205683205683205683205683205683e1_wp ] b(6,:) = [ 2.480158730158730158730158730158730159320e-1_wp, & -2.480158730158730158730158730158730159320e-1_wp, & 4.960317460317460317460317460317460316000e-1_wp, & 9.043040293040293040293040293040293040380e0_wp, & -5.952380952380952380952380952380952380920e0_wp, & -3.815628815628815628815628815628815628816e1_wp ] b(7,:) = [ 2.367424242424242424242424242424242424300e-1_wp, & -2.367424242424242424242424242424242424300e-1_wp, & 4.734848484848484848484848484848484847900e-1_wp, & 8.631993006993006993006993006993006993060e0_wp, & -5.681818181818181818181818181818181818140e0_wp, & -3.642191142191142191142191142191142191142e1_wp ] b(8,:) = [ 6.944444444444444444444444444444444444500e-2_wp, & -6.944444444444444444444444444444444444500e-2_wp, & 1.388888888888888888888888888888888889000e-1_wp, & 2.532051282051282051282051282051282051270e0_wp, & -1.666666666666666666666666666666666666680e0_wp, & -1.068376068376068376068376068376068376068e1_wp ] b(9,:) = [ -6.300000000000000000000000000000000000000e-38_wp, & 6.300000000000000000000000000000000000000e-38_wp, & -9.999999999999999999999999999999999998400e-1_wp, & 8.378986866791744840525328330206378986770e0_wp, & 3.999999999999999999999999999999999999970e0_wp, & -3.689806128830519074421513445903689806129e1_wp ] b(10,:) = [ 3.000000000000000000000000000000000000000e-38_wp, & -3.000000000000000000000000000000000000000e-38_wp, & -1.000000000000000000000000000000000000000e-37_wp, & -1.954807692307692307692307692307692307685e1_wp, & 1.600000000000000000000000000000000000005e1_wp, & 6.695156695156695156695156695156695156695e1_wp ] b(11,:) = [ 6.000000000000000000000000000000000000000e-38_wp, & -6.000000000000000000000000000000000000000e-38_wp, & -1.400000000000000000000000000000000000000e-37_wp, & -1.724591198275408801724591198275408801717e1_wp, & 3.000000000000000000000000000000000000000e-38_wp, & 7.789481473692000007789481473692000007789e1_wp ] b(12,:) = [ 2.600000000000000000000000000000000000000e-38_wp, & -2.600000000000000000000000000000000000000e-38_wp, & -3.000000000000000000000000000000000000000e-38_wp, & -2.511673624478502527283015087893136673620e0_wp, & -2.000000000000000000000000000000000000000e-38_wp, & 3.525155964180354424256863281253525155964e1_wp ] end procedure rkv65r_dense_coefficients !***************************************************************************************** !***************************************************************************************** !> ! The additional function evaluations for the dense output of [[rkv65r]]. ! These are stored in `me%funcs(:,10:12)`. module procedure rkv65r_dense_extra_stages real(wp),parameter :: c10 = 0.5_wp real(wp),parameter :: a1001 = .8384427865256290015733788189661258860630e-1_wp real(wp),parameter :: a1004 = .3246035582961368344664582851907258080863_wp real(wp),parameter :: a1005 = .2549935289615053977690984763272451791636_wp real(wp),parameter :: a1006 = -.2083300359719391581689534767215651661231_wp real(wp),parameter :: a1007 = .5594054201101300154461755512581943785651e-1_wp real(wp),parameter :: a1008 = .8068756448575586934400376221808702561795e-2_wp real(wp),parameter :: a1009 = -.1912062839785456270295909804064655015142e-1_wp real(wp),parameter :: c11 = 0.95_wp real(wp),parameter :: a1101 = .7702148739298602319228763650413263409805e-1_wp real(wp),parameter :: a1104 = .3660758457775556965544924825400616256967_wp real(wp),parameter :: a1105 = -.1607021082606043255126845533097619914719e-1_wp real(wp),parameter :: a1106 = .2565093765357899021858046924635752307515_wp real(wp),parameter :: a1107 = .2199918837388982340830864930547493396401_wp real(wp),parameter :: a1108 = .6247528293046205754409862849312624917397e-1_wp real(wp),parameter :: a1109 = -.3224866554963148100850147772466888021308e-1_wp real(wp),parameter :: a1110 = .16245e-1_wp real(wp),parameter :: c12 = 0.18_wp real(wp),parameter :: a1201 = .7800698111046033567745282845296820331593e-1_wp real(wp),parameter :: a1204 = .1811172942641407315027221158488386380810_wp real(wp),parameter :: a1205 = .5645435705965205370619057911637952373440e-1_wp real(wp),parameter :: a1206 = .4731750428236662786851828049771877530203e-1_wp real(wp),parameter :: a1207 = .1125370152418598004794505204249672699349_wp real(wp),parameter :: a1208 = .3437883070149390182751763620421893535207e-1_wp real(wp),parameter :: a1209 = .7109159773950015472868039358391738634605e-1_wp real(wp),parameter :: a1210 = -.2177060563732161788892773024487246345251_wp real(wp),parameter :: a1211 = -.1831975240262574269012550516802840975412_wp associate (t => me%dense_t, h => me%dense_h, x => me%dense_x, k => me%funcs) call me%f(t+c10*h, x+h*(a1001*k(:,1) + a1004*k(:,4) + a1005*k(:,5) + a1006*k(:,6) + & a1007*k(:,7) + a1008*k(:,8) + a1009*k(:,9)), k(:,10)) if (me%stopped) return call me%f(t+c11*h, x+h*(a1101*k(:,1) + a1104*k(:,4) + a1105*k(:,5) + a1106*k(:,6) + & a1107*k(:,7) + a1108*k(:,8) + a1109*k(:,9) + a1110*k(:,10)), k(:,11)) if (me%stopped) return call me%f(t+c12*h, x+h*(a1201*k(:,1) + a1204*k(:,4) + a1205*k(:,5) + a1206*k(:,6) + & a1207*k(:,7) + a1208*k(:,8) + a1209*k(:,9) + a1210*k(:,10) + & a1211*k(:,11)), k(:,12)) end associate end procedure rkv65r_dense_extra_stages !***************************************************************************************** !***************************************************************************************** !> ! Dense output for Verner's "most efficient" 7(6) method [[rkv76e]]. ! This is the interpolant of order 7 (`bi7`) from the coefficient file, which requires 6 additional ! function evaluations (see [[rkv76e_dense_extra_stages]]). ! !@note The interpolant is published in the polynomial form. Here it has been ! converted (exactly, using high precision arithmetic) ! to the form used in `DOPRI5` and `DOP853` (see [[dense_coefficients]]), which has ! much smaller roundoff error for these coefficients. This code is generated by ! `scripts/generate_dense.py`. ! !### References ! * J.H. Verner, "Numerically optimal Runge--Kutta pairs with interpolants", ! Numerical Algorithms, 53, (2010) pp. 383-396. [doi](https://doi.org/10.1007/s11075-009-9290-3) ! * Coefficients: ! * [RKV76.IIa.Efficient.00001675585.240712.CoeffsOnlyFLOAT](https://web.archive.org/web/https://www.sfu.ca/~jverner/RKV76.IIa.Efficient.00001675585.240712.CoeffsOnlyFLOAT) module procedure rkv76e_dense_coefficients form = dense_form_dopri allocate(b(16,7)) b = 0.0_wp b(1,:) = [ 4.715561848627222170431765108838175679000e-2_wp, & 9.528443815137277782956823489116182432100e-1_wp, & -9.056887630274555565913646978232364864100e-1_wp, & -6.507603564276723372369540449322724521640e0_wp, & 5.892770895323481590261325543930968013490e0_wp, & 1.364197995800900074206670061519242599876e1_wp, & -1.128508943677596076849260971330592556729e1_wp ] b(4,:) = [ 2.575056429843415189596436101037687581000e-1_wp, & -2.575056429843415189596436101037687581000e-1_wp, & 5.150112859686830379192872202075375162000e-1_wp, & 1.305036462316410850131105537042233350770e1_wp, & -4.165923804371634166612301904744348080780e1_wp, & -4.376652209331779289377480632786606178470e1_wp, & 1.511334163133759631666637566286118558727e2_wp ] b(5,:) = [ 2.621665397741262047713863095764527711300e-1_wp, & -2.621665397741262047713863095764527711300e-1_wp, & 5.243330795482524095427726191529055422000e-1_wp, & 1.328657848579128197486903476473833551600e1_wp, & -4.241327747606633744339811782365370765810e1_wp, & -4.455870373237099915437215158718941244490e1_wp, & 1.538689573553524332103726607127558372301e2_wp ] b(6,:) = [ 1.521609265673855740323133199165117534960e-1_wp, & -1.521609265673855740323133199165117534960e-1_wp, & 3.043218531347711480646266398330235073100e-1_wp, & 7.711503134801713250990254599635954060160e0_wp, & -2.461657999940846333853583565454619037286e1_wp, & -2.586178103582818667910666866851085994882e1_wp, & 8.930522995543090743868723335345884982106e1_wp ] b(7,:) = [ 4.939969170032484246907175893227876844100e-1_wp, & -4.939969170032484246907175893227876844100e-1_wp, & 9.879938340064968493814351786455753688000e-1_wp, & 2.503572277056216162882446737658922055630e1_wp, & -7.991877350645755558358970267411571783970e1_wp, & -8.396137160911947777093082122447357412020e1_wp, & 2.899332257332942839838911097057931625569e2_wp ] b(8,:) = [ -2.943031171403250441557244744092703429700e-1_wp, & 2.943031171403250441557244744092703429700e-1_wp, & -5.886062342806500883114489488185406856000e-1_wp, & -1.491525756058313178084518157903887432990e1_wp, & 4.761232985757159814364646622851366088860e1_wp, & 5.002074412496497926839983016688341959190e1_wp, & -1.727303332447657961368012115808350057416e2_wp ] b(9,:) = [ 8.131747232495109999734599440136761885300e-2_wp, & -8.131747232495109999734599440136761885300e-2_wp, & 1.626349446499021999946919888027352380800e-1_wp, & 4.121162751137067054203752922400085135340e0_wp, & -1.315553281643792667454779228069697606564e1_wp, & -1.382099012602666039388673116245715251628e1_wp, & 4.772628380491678816433752531852051270954e1_wp ] b(11,:) = [ 5.800000000000000000000000000000000000000e-38_wp, & -5.800000000000000000000000000000000000000e-38_wp, & -1.000000000000000000000000000000000000210e0_wp, & -2.274137494621028294117921462492887425110e0_wp, & 1.350004216864734643978305010464529379608e1_wp, & 8.107773996135029032330387735860697715060e0_wp, & -4.025624390076203635639388769027397381848e1_wp ] b(12,:) = [ 5.000000000000000000000000000000000000000e-39_wp, & -5.000000000000000000000000000000000000000e-39_wp, & -1.200000000000000000000000000000000000000e-38_wp, & 5.921404526401594125717004582969923833258e0_wp, & 3.177881772328029994825234392478111343014e0_wp, & -1.888807125659099312686665549661614362658e1_wp, & -1.895782173277762779516963812746246780056e1_wp ] b(13,:) = [ -3.000000000000000000000000000000000000000e-39_wp, & 3.000000000000000000000000000000000000000e-39_wp, & -5.000000000000000000000000000000000000000e-39_wp, & 8.487477767791544489695951710080762138470e0_wp, & 4.555034328602948446656001189774133752145e0_wp, & -2.707332089067715696671560778020912445733e1_wp, & -2.717329811960778939941155595778967136139e1_wp ] b(14,:) = [ -3.000000000000000000000000000000000000000e-38_wp, & 3.000000000000000000000000000000000000000e-38_wp, & 1.000000000000000000000000000000000000000e-37_wp, & -1.918681975136084566194177555055440130000e1_wp, & 4.445558244677403600994761077346504577710e1_wp, & 8.432060941152660306911484033272216318690e1_wp, & -1.767759230365517641222212235669933633475e2_wp ] b(15,:) = [ 9.000000000000000000000000000000000000000e-38_wp, & -9.000000000000000000000000000000000000000e-38_wp, & -2.000000000000000000000000000000000000000e-37_wp, & -1.588375234934266311207464684269517278940e1_wp, & 4.483823277676953244042532492832832592390e1_wp, & 3.303461482677318813218485106899798191110e1_wp, & -1.234864347071902008205124473040986284057e2_wp ] b(16,:) = [ 4.000000000000000000000000000000000000000e-38_wp, & -4.000000000000000000000000000000000000000e-38_wp, & -2.000000000000000000000000000000000000000e-37_wp, & -1.884664333946507880426245544273255438120e1_wp, & 3.773152759606965164064945431932053325010e1_wp, & 6.880503842652246674155683232766564049480e1_wp, & -1.613019689839392005649497117783811821479e2_wp ] end procedure rkv76e_dense_coefficients !***************************************************************************************** !***************************************************************************************** !> ! The additional function evaluations for the dense output of [[rkv76e]]. ! These are stored in `me%funcs(:,11:16)`. Stage 11 is `f` at the end of the step, ! which is cached for the first stage of the next step (see [[set_fsal_cache]]). module procedure rkv76e_dense_extra_stages real(wp),parameter :: c12 = .6716626365038747067708664674609860468131_wp real(wp),parameter :: a1201 = .4608828480030097784573810160682004587659e-1_wp real(wp),parameter :: a1204 = .2609521155912221411807714016159194601215_wp real(wp),parameter :: a1205 = .2593539886572148954104325364777172658208_wp real(wp),parameter :: a1206 = .1071284853231358590149526776164009620565_wp real(wp),parameter :: a1207 = -.7335737007803501940010147823741528427117e-1_wp real(wp),parameter :: a1208 = .9516817597231148848177657505563031988827e-1_wp real(wp),parameter :: a1209 = -.5557470923733997764040425874100454845262e-1_wp real(wp),parameter :: a1211 = .3190366547506434187770091206691782577319e-1_wp real(wp),parameter :: c13 = .1250000000000000000000000000000000000000_wp real(wp),parameter :: a1301 = .5837333975598387183322580009499934069289e-1_wp real(wp),parameter :: a1304 = .8577524236694000499244056599099302286222e-1_wp real(wp),parameter :: a1305 = -.3429204581002122987958198322824658635148e-1_wp real(wp),parameter :: a1306 = -.2818259731103036151485354567470271399566e-1_wp real(wp),parameter :: a1307 = -.2078116148802372182126410976576048400947e-1_wp real(wp),parameter :: a1308 = -.3670740144441019619026795845974788863064e-2_wp real(wp),parameter :: a1309 = .1140645258436468819102382339952341616568e-1_wp real(wp),parameter :: a1311 = -.2750153913374408098701890663991826389698e-2_wp real(wp),parameter :: a1312 = .5912166395960217591673813569316061988858e-1_wp real(wp),parameter :: c14 = .2500000000000000000000000000000000000000_wp real(wp),parameter :: a1401 = .3781620167894484389097734217537475940628e-1_wp real(wp),parameter :: a1404 = .3134678532631074103775642636872622557979e-1_wp real(wp),parameter :: a1405 = .2225110077526238187226161613568406825821e-1_wp real(wp),parameter :: a1406 = .9958131955421000529301204622879718202274e-2_wp real(wp),parameter :: a1407 = .1522309604389992558287426004143093193551e-1_wp real(wp),parameter :: a1408 = -.7547047454118317125873194132609456732396e-2_wp real(wp),parameter :: a1409 = .1315888206546629705431596397483148558901e-2_wp real(wp),parameter :: a1411 = -.2347519033359896919332911522715968811661e-2_wp real(wp),parameter :: a1412 = -.2426983932792809284254784793342448202006e-1_wp real(wp),parameter :: a1413 = .1662532018290207842691515078471710556231_wp real(wp),parameter :: c15 = .5300000000000000000000000000000000000000_wp real(wp),parameter :: a1501 = .5345075085525607281977266793633803424059e-1_wp real(wp),parameter :: a1504 = .3430173382138450672760919883610065904944_wp real(wp),parameter :: a1505 = .2076742726203398574300787732117594341897_wp real(wp),parameter :: a1506 = .7722675619395836334886386873715328439377e-1_wp real(wp),parameter :: a1507 = .1327268798561135321764135364638696721663e-3_wp real(wp),parameter :: a1508 = .2221995820246240129102290486291650170691e-1_wp real(wp),parameter :: a1509 = -.1741022886840232443411414202292540569931e-1_wp real(wp),parameter :: a1511 = .7197593769915030785550602945672682421140e-2_wp real(wp),parameter :: a1512 = -.8645504868462665761280182814730150657031e-1_wp real(wp),parameter :: a1513 = -.7705411918260392443664124942108348484897e-1_wp real(wp),parameter :: c16 = .7900000000000000000000000000000000000000_wp real(wp),parameter :: a1601 = .5681087911635884235009187585387593516296e-1_wp real(wp),parameter :: a1604 = .3685221901618432443784221615000295848158_wp real(wp),parameter :: a1605 = .2294060530103075331008808070568391539642_wp real(wp),parameter :: a1606 = .8854425014313802531960211711477635403509e-1_wp real(wp),parameter :: a1607 = .2937879283098644165855344854841522051205e-1_wp real(wp),parameter :: a1608 = .5463456745004647616602723887000652830690e-2_wp real(wp),parameter :: a1609 = -.1311750049635594995549114647314453663034e-1_wp real(wp),parameter :: a1611 = .5109093689745296784307624438438730642345e-2_wp real(wp),parameter :: a1612 = .1241964089587753256730552835816994961283_wp real(wp),parameter :: a1613 = -.1043136241598034069260248955079305914610_wp associate (t => me%dense_t, h => me%dense_h, x => me%dense_x, k => me%funcs) ! f at the end of the step (this is also the first stage of the next step, ! so it is saved in the function evaluation cache): call me%set_fsal_cache(t+h, me%dense_xf, k(:,11)) if (me%stopped) return call me%f(t+c12*h, x+h*(a1201*k(:,1) + a1204*k(:,4) + a1205*k(:,5) + a1206*k(:,6) + & a1207*k(:,7) + a1208*k(:,8) + a1209*k(:,9) + a1211*k(:,11)), k(:,12)) if (me%stopped) return call me%f(t+c13*h, x+h*(a1301*k(:,1) + a1304*k(:,4) + a1305*k(:,5) + a1306*k(:,6) + & a1307*k(:,7) + a1308*k(:,8) + a1309*k(:,9) + a1311*k(:,11) + & a1312*k(:,12)), k(:,13)) if (me%stopped) return call me%f(t+c14*h, x+h*(a1401*k(:,1) + a1404*k(:,4) + a1405*k(:,5) + a1406*k(:,6) + & a1407*k(:,7) + a1408*k(:,8) + a1409*k(:,9) + a1411*k(:,11) + & a1412*k(:,12) + a1413*k(:,13)), k(:,14)) if (me%stopped) return call me%f(t+c15*h, x+h*(a1501*k(:,1) + a1504*k(:,4) + a1505*k(:,5) + a1506*k(:,6) + & a1507*k(:,7) + a1508*k(:,8) + a1509*k(:,9) + a1511*k(:,11) + & a1512*k(:,12) + a1513*k(:,13)), k(:,15)) if (me%stopped) return call me%f(t+c16*h, x+h*(a1601*k(:,1) + a1604*k(:,4) + a1605*k(:,5) + a1606*k(:,6) + & a1607*k(:,7) + a1608*k(:,8) + a1609*k(:,9) + a1611*k(:,11) + & a1612*k(:,12) + a1613*k(:,13)), k(:,16)) end associate end procedure rkv76e_dense_extra_stages !***************************************************************************************** !***************************************************************************************** !> ! Dense output for Verner's "most robust" 7(6) method [[rkv76r]]. ! This is the interpolant of order 7 (`bi7`) from the coefficient file, which requires 6 additional ! function evaluations (see [[rkv76r_dense_extra_stages]]). ! !@note The interpolant is published in the polynomial form. Here it has been ! converted (exactly, using high precision arithmetic) ! to the form used in `DOPRI5` and `DOP853` (see [[dense_coefficients]]), which has ! much smaller roundoff error for these coefficients. This code is generated by ! `scripts/generate_dense.py`. ! !### References ! * J.H. Verner, "Numerically optimal Runge--Kutta pairs with interpolants", ! Numerical Algorithms, 53, (2010) pp. 383-396. [doi](https://doi.org/10.1007/s11075-009-9290-3) ! * Coefficients: ! * [RKV76.IIa.Robust.000027015646.240702.CoeffsOnlyFLOAT](https://web.archive.org/web/https://www.sfu.ca/~jverner/RKV76.IIa.Robust.000027015646.240702.CoeffsOnlyFLOAT) module procedure rkv76r_dense_coefficients form = dense_form_dopri allocate(b(16,7)) b = 0.0_wp b(1,:) = [ 4.742583783370675608356917271757453471000e-2_wp, & 9.525741621662932439164308272824254652900e-1_wp, & -9.051483243325864878328616545648509305800e-1_wp, & -9.247769696609901771200432943507423179840e0_wp, & 1.008484494868409955362207214003549261074e1_wp, & 2.423072825123594614621952144410433879113e1_wp, & -2.643772965023979233919406894246137029559e1_wp ] b(4,:) = [ 2.562236165937056265996172745827462345060e-1_wp, & -2.562236165937056265996172745827462345060e-1_wp, & 5.124472331874112531992345491654924687500e-1_wp, & 3.875938560505704911019757003525593265100e0_wp, & -3.320036334624272353976777224485217875430e1_wp, & -1.423079379182164741551616784175434589540e1_wp, & 1.344969468554404353215667450076794927270e2_wp ] b(5,:) = [ 2.695137683307420661947381725807595289110e-1_wp, & -2.695137683307420661947381725807595289110e-1_wp, & 5.390275366614841323894763451615190575800e-1_wp, & 4.076980963533813785615300863577411294600e0_wp, & -3.492244452073484915911976414160485084100e1_wp, & -1.496893577633546677681299933233158046240e1_wp, & 1.414732156929510235031641753915382786420e2_wp ] b(6,:) = [ 1.268662240909278284598913836473917324630e-1_wp, & -1.268662240909278284598913836473917324630e-1_wp, & 2.537324481818556569197827672947834649900e-1_wp, & 1.919127114498253053428817670578791418680e0_wp, & -1.643878418461184000017634575192408471168e1_wp, & -7.046216497083747679445264376713159504310e0_wp, & 6.659464114256471991438515546495341015153e1_wp ] b(7,:) = [ 2.488722594206007162204644942764749276710e-1_wp, & -2.488722594206007162204644942764749276710e-1_wp, & 4.977445188412014324409289885529498554200e-1_wp, & 3.764733320652777545101390681461832815200e0_wp, & -3.224780584010889378401084003978592072780e1_wp, & -1.382249556618862871458536273693198262300e1_wp, & 1.306380711273911686533709699607900281508e2_wp ] b(8,:) = [ 3.074483740820063133530438847909918477650e-3_wp, & -3.074483740820063133530438847909918477650e-3_wp, & 6.148967481640126267060877695819836951200e-3_wp, & 4.650824246067975644877073222847602720000e-2_wp, & -3.983784892834471862756096069951491819120e-1_wp, & -1.707584363751173695218801401601925022450e-1_wp, & 1.613858557592266711452377540896206369143e0_wp ] b(9,:) = [ 4.802380998949694330818906334714312331160e-2_wp, & -4.802380998949694330818906334714312331160e-2_wp, & 9.604761997899388661637812669428624669000e-2_wp, & 7.264644041608727817394616659296532502400e-1_wp, & -6.222720458475440928086060638779301997080e0_wp, & -2.667267545996165184478111760572809322090e0_wp, & 2.520866696763272772148596229997893321296e1_wp ] b(11,:) = [ -3.900000000000000000000000000000000000000e-38_wp, & 3.900000000000000000000000000000000000000e-38_wp, & -9.999999999999999999999999999999999998700e-1_wp, & -1.211426977068005270830123084832586657300e-1_wp, & 1.123550754503722076227627104740100843398e1_wp, & 1.059363194376045721496642931223661793660e0_wp, & -3.490298534551717356011931265592436984878e1_wp ] b(12,:) = [ -4.000000000000000000000000000000000000000e-39_wp, & 4.000000000000000000000000000000000000000e-39_wp, & 5.000000000000000000000000000000000000000e-39_wp, & -2.689843926762734958111868061662373328226e0_wp, & 2.859134803272277717713314303305459761476e0_wp, & 9.953676535514782624714663281792581954677e0_wp, & -9.405048694974597755635902313504801846941e0_wp ] b(13,:) = [ -1.000000000000000000000000000000000000000e-38_wp, & 1.000000000000000000000000000000000000000e-38_wp, & 2.000000000000000000000000000000000000000e-38_wp, & 1.428744953003485588757315864042649968174e1_wp, & -1.518665964031676989448335743597781784354e1_wp, & -5.287022407676068466555335510880435488293e1_wp, & 4.995611723788411149501104419729545343269e1_wp ] b(14,:) = [ 1.800000000000000000000000000000000000000e-38_wp, & -1.800000000000000000000000000000000000000e-38_wp, & -1.000000000000000000000000000000000000000e-37_wp, & -7.739397697899142305765936845963446864600e0_wp, & 4.479929485570897843490005674735774971440e1_wp, & 4.914838447704658193123231430817031966380e1_wp, & -1.987641278148321659042765024584136503760e2_wp ] b(15,:) = [ 7.500000000000000000000000000000000000000e-38_wp, & -7.500000000000000000000000000000000000000e-38_wp, & -3.100000000000000000000000000000000000000e-37_wp, & -3.698276994746987304944919677286868631400e0_wp, & 4.319682040705424529766877626541831495990e1_wp, & -3.536486384630642938900017250523273227800e0_wp, & -1.420948039705731753212788692941391939461e2_wp ] b(16,:) = [ -8.200000000000000000000000000000000000000e-38_wp, & 8.200000000000000000000000000000000000000e-38_wp, & 2.800000000000000000000000000000000000000e-37_wp, & -5.200771122121390853820487420824887082800e0_wp, & 2.644155392001714272573925935640127857710e1_wp, & 2.492102561701874432115001658250079621680e1_wp, & -1.383768221053195484399317741986884163727e2_wp ] end procedure rkv76r_dense_coefficients !***************************************************************************************** !***************************************************************************************** !> ! The additional function evaluations for the dense output of [[rkv76r]]. ! These are stored in `me%funcs(:,11:16)`. Stage 11 is `f` at the end of the step, ! which is cached for the first stage of the next step (see [[set_fsal_cache]]). module procedure rkv76r_dense_extra_stages real(wp),parameter :: c12 = .3067207177086018949109476427952022863650_wp real(wp),parameter :: a1201 = .5084001823323651000647096877290063127003e-1_wp real(wp),parameter :: a1204 = .2340953776302340927060690222130822879757_wp real(wp),parameter :: a1205 = .2957282565868639393723036073846174594572e-1_wp real(wp),parameter :: a1206 = -.7891822081774038060777822301275761921904e-2_wp real(wp),parameter :: a1207 = -.3124515117284615723495332436165903196571e-3_wp real(wp),parameter :: a1208 = .3075900699836012263713563630912227790408e-3_wp real(wp),parameter :: a1209 = .1346518474743767274498197636413097501530e-1_wp real(wp),parameter :: a1211 = -.1335600503747387607704868611157222437932e-1_wp real(wp),parameter :: c13 = .7875e-1_wp real(wp),parameter :: a1301 = .4571638671066279127171353905223268681038e-1_wp real(wp),parameter :: a1304 = .6135139403136968332208211375799160523016e-1_wp real(wp),parameter :: a1305 = .1491507057398255131431717710757995218460e-1_wp real(wp),parameter :: a1306 = .7635076575109071067007610291932885992123e-3_wp real(wp),parameter :: a1307 = -.2047836064685163003552627495384108018023e-2_wp real(wp),parameter :: a1308 = .5715454045443600137343865669766089985382e-5_wp real(wp),parameter :: a1309 = .1261351063159691297340795743480869575660e-2_wp real(wp),parameter :: a1311 = -.6786312294583558454434945884425128592581e-3_wp real(wp),parameter :: a1312 = -.4253695819658754906329560847232154761272e-1_wp real(wp),parameter :: c14 = .2_wp real(wp),parameter :: a1401 = .4895720837021160646063704247238281140243e-2_wp real(wp),parameter :: a1404 = -.2715965016132070652354391268389482837445e-1_wp real(wp),parameter :: a1405 = -.1471118582892732021356233028007905752798e-1_wp real(wp),parameter :: a1406 = -.3548741829515423281693629852942855954300e-2_wp real(wp),parameter :: a1407 = .3086771782083923184307182280403431753151e-2_wp real(wp),parameter :: a1408 = .8149105769962305218267679359657679200786e-4_wp real(wp),parameter :: a1409 = .1992488082891027877223378212655417925868e-2_wp real(wp),parameter :: a1411 = -.2859268825507094536557198895337022154448e-2_wp real(wp),parameter :: a1412 = .5586163713000076148900856473826299382047e-1_wp real(wp),parameter :: a1413 = .1823607377555740483065715654400970625794_wp real(wp),parameter :: c15 = .55_wp real(wp),parameter :: a1501 = .6529904148291268153192118584675719879592e-1_wp real(wp),parameter :: a1504 = .2572420423912244394795082675403296906296_wp real(wp),parameter :: a1505 = .1483811462482641701867308875408674270889_wp real(wp),parameter :: a1506 = .4007284294116317186672184994560870753991e-1_wp real(wp),parameter :: a1507 = -.1000336192875383976289218926700770766316e-1_wp real(wp),parameter :: a1508 = -.5059442698140261792783567717843895596137e-3_wp real(wp),parameter :: a1509 = -.1424880415164983774939393220447878176376e-1_wp real(wp),parameter :: a1511 = .1681180958703984153237987900953773537091e-1_wp real(wp),parameter :: a1512 = .8770029800174049628643329724269360224271e-1_wp real(wp),parameter :: a1513 = -.4074907030212709719213088888252348268140e-1_wp real(wp),parameter :: c16 = .8_wp real(wp),parameter :: a1601 = .1045152970430589502927487092331130990894_wp real(wp),parameter :: a1604 = .5969208192685441448792898447593605742233_wp real(wp),parameter :: a1605 = .3959219441466996682329003000423645782513_wp real(wp),parameter :: a1606 = .1298549335245192828418583437596723989380_wp real(wp),parameter :: a1607 = .8653312743253928887016075996436582821518e-1_wp real(wp),parameter :: a1608 = .3432143531736777732699173398591878436542e-3_wp real(wp),parameter :: a1609 = -.6684380434300498891936757211031674244336e-2_wp real(wp),parameter :: a1611 = -.2656538838207228881702404465288483040129e-2_wp real(wp),parameter :: a1612 = -.2623495692645475071531194269140313972733_wp real(wp),parameter :: a1613 = -.2423988472314797779634692865083841120031_wp associate (t => me%dense_t, h => me%dense_h, x => me%dense_x, k => me%funcs) ! f at the end of the step (this is also the first stage of the next step, ! so it is saved in the function evaluation cache): call me%set_fsal_cache(t+h, me%dense_xf, k(:,11)) if (me%stopped) return call me%f(t+c12*h, x+h*(a1201*k(:,1) + a1204*k(:,4) + a1205*k(:,5) + a1206*k(:,6) + & a1207*k(:,7) + a1208*k(:,8) + a1209*k(:,9) + a1211*k(:,11)), k(:,12)) if (me%stopped) return call me%f(t+c13*h, x+h*(a1301*k(:,1) + a1304*k(:,4) + a1305*k(:,5) + a1306*k(:,6) + & a1307*k(:,7) + a1308*k(:,8) + a1309*k(:,9) + a1311*k(:,11) + & a1312*k(:,12)), k(:,13)) if (me%stopped) return call me%f(t+c14*h, x+h*(a1401*k(:,1) + a1404*k(:,4) + a1405*k(:,5) + a1406*k(:,6) + & a1407*k(:,7) + a1408*k(:,8) + a1409*k(:,9) + a1411*k(:,11) + & a1412*k(:,12) + a1413*k(:,13)), k(:,14)) if (me%stopped) return call me%f(t+c15*h, x+h*(a1501*k(:,1) + a1504*k(:,4) + a1505*k(:,5) + a1506*k(:,6) + & a1507*k(:,7) + a1508*k(:,8) + a1509*k(:,9) + a1511*k(:,11) + & a1512*k(:,12) + a1513*k(:,13)), k(:,15)) if (me%stopped) return call me%f(t+c16*h, x+h*(a1601*k(:,1) + a1604*k(:,4) + a1605*k(:,5) + a1606*k(:,6) + & a1607*k(:,7) + a1608*k(:,8) + a1609*k(:,9) + a1611*k(:,11) + & a1612*k(:,12) + a1613*k(:,13)), k(:,16)) end associate end procedure rkv76r_dense_extra_stages !***************************************************************************************** !***************************************************************************************** !> ! Dense output for Verner's "most efficient" 8(7) method [[rkv87e]]. ! This is the interpolant of order 8 (`bi8`) from the coefficient file, which requires 8 additional ! function evaluations (see [[rkv87e_dense_extra_stages]]). ! !@note The interpolant is published in the polynomial form. Here it has been ! converted (exactly, using high precision arithmetic) ! to the form used in `DOPRI5` and `DOP853` (see [[dense_coefficients]]), which has ! much smaller roundoff error for these coefficients. This code is generated by ! `scripts/generate_dense.py`. ! !### References ! * J.H. Verner, "Numerically optimal Runge--Kutta pairs with interpolants", ! Numerical Algorithms, 53, (2010) pp. 383-396. [doi](https://doi.org/10.1007/s11075-009-9290-3) ! * Coefficients: ! * [RKV87.IIa.Efficient.000000282866.081208.CoeffsOnlyFLOAT](https://web.archive.org/web/https://www.sfu.ca/~jverner/RKV87.IIa.Efficient.000000282866.081208.CoeffsOnlyFLOAT) module procedure rkv87e_dense_coefficients form = dense_form_dopri allocate(b(21,8)) b = 0.0_wp b(1,:) = [ 4.427989419007951074716746668098518858000e-2_wp, & 9.557201058099204892528325333190148114200e-1_wp, & -9.114402116198409785056650666380296225600e-1_wp, & -8.171994333124757515048952495580225724400e0_wp, & 6.218117224597786014516545699875570773300e0_wp, & 3.031855885615617192625492971517856421740e1_wp, & -1.312333874063877708979425960087686120150e1_wp, & -4.037066163211959429657758663355894180800e1_wp ] b(6,:) = [ 3.541049391724448744815552028733568455000e-1_wp, & -3.541049391724448744815552028733568455000e-1_wp, & 7.082098783448897489631104057467136360000e-1_wp, & 1.571353590946602791833084348213441851230e2_wp, & -8.684453516225225897281302987797444548800e1_wp, & -1.142137953989572306294936176630967321549e3_wp, & 3.133840664884543279410228200360555258310e2_wp, & 2.170776389952444021264933974457050280938e3_wp ] b(7,:) = [ 2.479692154956437828667629415370663006000e-1_wp, & -2.479692154956437828667629415370663006000e-1_wp, & 4.959384309912875657335258830741326110000e-1_wp, & 1.100372443614908888134908003296047727080e2_wp, & -6.081465936226426911889323723839788873900e1_wp, & -7.998054280193804142458536875932032445150e2_wp, & 2.194535927614752733203637065660331899130e2_wp, & 1.520130500554341433782477059435641543286e3_wp ] b(8,:) = [ -1.569420203883808405099207034271191196500e1_wp, & 1.569420203883808405099207034271191196500e1_wp, & -3.138840407767616810198414068542382483000e1_wp, & -6.964359432094799837481091812226804369250e3_wp, & 3.849016294408765508661118711572209842100e3_wp, & 5.062042864476473755036991016932779475806e4_wp, & -1.388942178190660596574923609338849580196e4_wp, & -9.621047174510666745715793578288559674281e4_wp ] b(9,:) = [ 2.508406496555856261343930031237186303000e1_wp, & -2.508406496555856261343930031237186303000e1_wp, & 5.016812993111712522687860062474372450000e1_wp, & 1.113114537495147165247342692512124943360e4_wp, & -6.151888101320182430040782986447821661900e3_wp, & -8.090670156835237541683799415555856792670e4_wp, & 2.219948216858710770927295911580358782450e4_wp, & 1.537733309185793956820086499888593205888e5_wp ] b(10,:) = [ -3.173836778626027646833156112007297726000e1_wp, & 3.173836778626027646833156112007297726000e1_wp, & -6.347673557252055293666312224014595500000e1_wp, & -1.408401653709690307884869527168714351650e4_wp, & 7.783861483683196951282356057244934082800e3_wp, & 1.023696380261861982074552448790615600052e5_wp, & -2.808856262725180126416959436145273050140e4_wp, & -1.945663328138133045593670938904445416121e5_wp ] b(11,:) = [ 2.293828327398878395231483560344797028000e1_wp, & -2.293828327398878395231483560344797028000e1_wp, & 4.587656654797756790462967120689594010000e1_wp, & 1.017894691809977833541568779743218681630e4_wp, & -5.625633330631149943201655806994276135700e3_wp, & -7.398564953037362656904975446516125224370e4_wp, & 2.030045812822167651338439901004609940560e4_wp, & 1.406190037156383022701488158207833280861e5_wp ] b(12,:) = [ -2.361324633071542145259900641263517579000e-1_wp, & 2.361324633071542145259900641263517579000e-1_wp, & -4.722649266143084290519801282527035280000e-1_wp, & -1.047846423785882160944152250822090481310e2_wp, & 5.791168589896686562832131791261995638500e1_wp, & 7.616269039975525878738694882178573744430e2_wp, & -2.089780271183806579187902254943500565490e2_wp, & -1.447567428588892382130036646632729629570e3_wp ] b(14,:) = [ -2.800000000000000000000000000000000000000e-37_wp, & 2.800000000000000000000000000000000000000e-37_wp, & -9.999999999999999999999999999999999985000e-1_wp, & -1.386361337326743122469601010648237947730e1_wp, & 1.523445924570790606643731098536786401900e1_wp, & 1.018019076567058926442727278417486740185e2_wp, & -4.493783698283162426574924394147145608250e1_wp, & -1.929796145225588132959328212730088960570e2_wp ] b(15,:) = [ -1.720000000000000000000000000000000000000e-36_wp, & 1.720000000000000000000000000000000000000e-36_wp, & 6.100000000000000000000000000000000000000e-36_wp, & 1.434968575290546223276673100484047073210e1_wp, & -6.823976739116491511183354682495270228000e0_wp, & -1.147715862032391499084617729429976110027e2_wp, & 2.729590695646596604473341872998108089670e1_wp, & 2.298729377727072096359824945955196848017e2_wp ] b(16,:) = [ 2.200000000000000000000000000000000000000e-36_wp, & -2.200000000000000000000000000000000000000e-36_wp, & -5.000000000000000000000000000000000000000e-36_wp, & -1.025452470111040085560664290210906322490e2_wp, & 4.876527558584640178619311834243828117000e1_wp, & 8.201768916566261752269372741604835638680e2_wp, & -1.950611023433856071447724733697531246760e2_wp, & -1.642710341604349689799450723704711058784e3_wp ] b(17,:) = [ 4.000000000000000000000000000000000000000e-37_wp, & -4.000000000000000000000000000000000000000e-37_wp, & -4.400000000000000000000000000000000000000e-36_wp, & -3.813206313286473398334122725888547021350e1_wp, & 1.813365925316544375842493172333461372500e1_wp, & 3.049876803103379564611422864078234159030e2_wp, & -7.253463701266177503369972689333845488000e1_wp, & -6.108516609091004863949139257772330194915e2_wp ] b(18,:) = [ -3.800000000000000000000000000000000000000e-37_wp, & 3.800000000000000000000000000000000000000e-37_wp, & 1.200000000000000000000000000000000000000e-36_wp, & -6.638279583069588062871084016403504860100e1_wp, & 5.399983561289169607183189822084630570200e1_wp, & 4.090643411138268706944841483866826765826e2_wp, & -2.159993424515667842873275928833852228116e2_wp, & -5.711622272434449401356158886201861909946e2_wp ] b(19,:) = [ -1.640000000000000000000000000000000000000e-36_wp, & 1.640000000000000000000000000000000000000e-36_wp, & 5.600000000000000000000000000000000000000e-36_wp, & -9.041887573173058787343992868450872086300e1_wp, & 3.029342771398165999553907388317154648800e1_wp, & 7.208192092273725349071999068697808498760e2_wp, & -5.717371085592663998215629553268618596400e1_wp, & -1.592261308015418013416409177206823360972e3_wp ] b(20,:) = [ 4.900000000000000000000000000000000000000e-37_wp, & -4.900000000000000000000000000000000000000e-37_wp, & -3.000000000000000000000000000000000000000e-36_wp, & -5.973884363038871206457816967313835076550e1_wp, & 4.983159781045732460927805546730901932800e1_wp, & 3.755777296179377040175517228331666569468e2_wp, & -1.993263912418292984371122218692360772978e2_wp, & -6.049687555793670790184208565420961249773e2_wp ] b(21,:) = [ 2.500000000000000000000000000000000000000e-37_wp, & -2.500000000000000000000000000000000000000e-37_wp, & -1.500000000000000000000000000000000000000e-36_wp, & -5.920053764683937384859682230934791521200e1_wp, & 1.873876677738784897131139418885832105800e1_wp, & 4.346261735507422047601327250649118076426e2_wp, & -7.495506710955139588524557675543328422660e1_wp, & -8.334379054819676018284720384103746063216e2_wp ] end procedure rkv87e_dense_coefficients !***************************************************************************************** !***************************************************************************************** !> ! The additional function evaluations for the dense output of [[rkv87e]]. ! These are stored in `me%funcs(:,14:21)`. Stage 14 is `f` at the end of the step, ! which is cached for the first stage of the next step (see [[set_fsal_cache]]). module procedure rkv87e_dense_extra_stages real(wp),parameter :: c15 = .3110177634953863863927417318829099695921_wp real(wp),parameter :: a1501 = .4620700646754963101730413150238116432863e-1_wp real(wp),parameter :: a1506 = .4503904160842480866828520384400679697151e-1_wp real(wp),parameter :: a1507 = .2336816697713424410788701065340221126565_wp real(wp),parameter :: a1508 = 37.83901368421067410780338220861855254153_wp real(wp),parameter :: a1509 = -15.94911328945424610266139490307397370835_wp real(wp),parameter :: a1510 = 23.02836835181610285142510596329590091940_wp real(wp),parameter :: a1511 = -44.85578507769412524816130998016948002745_wp real(wp),parameter :: a1512 = -.6379858768647444009509067402330140781326e-1_wp real(wp),parameter :: a1514 = -.1259503554386166268241032464519842162533e-1_wp real(wp),parameter :: c16 = .1725_wp real(wp),parameter :: a1601 = .5037946855482040993065158747220696112586e-1_wp real(wp),parameter :: a1606 = .4109836131046079339916530614028848248545e-1_wp real(wp),parameter :: a1607 = .1718054153348195783296309209549424619697_wp real(wp),parameter :: a1608 = 4.61410531998151886974342237185977124648_wp real(wp),parameter :: a1609 = -1.791667883085396449712744996746836471721_wp real(wp),parameter :: a1610 = 2.531658930485041408462243518792913614971_wp real(wp),parameter :: a1611 = -5.32497786020573071925718815977276269909_wp real(wp),parameter :: a1612 = -.3065532595385634734924449496356513113607e-1_wp real(wp),parameter :: a1614 = -.5254479979429613570549519094377878106127e-2_wp real(wp),parameter :: a1615 = -.8399194644224792997538653464258058697156e-1_wp real(wp),parameter :: c17 = .7846_wp real(wp),parameter :: a1701 = .4082897132997079620207118756242653796386e-1_wp real(wp),parameter :: a1706 = .4244479514247632218892086657732332485609_wp real(wp),parameter :: a1707 = .2326091531275234539465100096964845486081_wp real(wp),parameter :: a1708 = 2.677982520711806062780528871014035962908_wp real(wp),parameter :: a1709 = .7420826657338945216477607044022963622057_wp real(wp),parameter :: a1710 = .1460377847941461193920992339731312296021_wp real(wp),parameter :: a1711 = -3.579344509890565218033356743825917680543_wp real(wp),parameter :: a1712 = .1138844389600173704531638716149985665239_wp real(wp),parameter :: a1714 = .1267790651033190047378693537615687232109e-1_wp real(wp),parameter :: a1715 = -.7443436349946674429752785032561552478382e-1_wp real(wp),parameter :: a1716 = .4782748079757851554575511473876987663388e-1_wp real(wp),parameter :: c18 = .37_wp real(wp),parameter :: a1801 = .5212682393668413629928136927994514676607e-1_wp real(wp),parameter :: a1806 = .5392508396744797718209106862347065628649e-1_wp real(wp),parameter :: a1807 = .1660758097434640828541930599928251901718e-1_wp real(wp),parameter :: a1808 = -4.454485757926779655418936993298463071587_wp real(wp),parameter :: a1809 = 6.835218278632146381711296817968152631469_wp real(wp),parameter :: a1810 = -8.711334822181993739847172734848837971169_wp real(wp),parameter :: a1811 = 6.491635839232917053651267142703105653517_wp real(wp),parameter :: a1812 = -.7072551809844346422069985227700294651922e-1_wp real(wp),parameter :: a1814 = -.1854031491993216429111842937941202966440e-1_wp real(wp),parameter :: a1815 = .2350402105435384645116542087045962190647e-1_wp real(wp),parameter :: a1816 = .2344795103407822090556377813402774776461_wp real(wp),parameter :: a1817 = -.8241072501152898885823089698097768766651e-1_wp real(wp),parameter :: c19 = .5_wp real(wp),parameter :: a1901 = .5020102870355713598699964419977883461362e-1_wp real(wp),parameter :: a1906 = .1552209034795498114932226104700567642339_wp real(wp),parameter :: a1907 = .1264268424089234914713091134864747506300_wp real(wp),parameter :: a1908 = -5.14920630353984701704917414605721854951_wp real(wp),parameter :: a1909 = 8.46834099903692926607453176331494311551_wp real(wp),parameter :: a1910 = -10.66213068108149527544209836207095498430_wp real(wp),parameter :: a1911 = 7.54183322495972836290996201569018333903_wp real(wp),parameter :: a1912 = -.743696811383214243944066492459357053774e-1_wp real(wp),parameter :: a1914 = -.2055887686618382619339821759221121764364e-1_wp real(wp),parameter :: a1915 = .775379526471029807261782993777862395844e-1_wp real(wp),parameter :: a1916 = .1046259220352544296313761971333987587377_wp real(wp),parameter :: a1917 = -.1179213306451979352145022687063013455111_wp real(wp),parameter :: c20 = .7_wp real(wp),parameter :: a2001 = .3737341446457825692757506548800094134977e-1_wp real(wp),parameter :: a2006 = .3504930705338316406767087468339071089224_wp real(wp),parameter :: a2007 = .4922652819373025433298989824173484805373_wp real(wp),parameter :: a2008 = 8.553695439359312242284304421725315855379_wp real(wp),parameter :: a2009 = -10.35317299030591348532574006719207803272_wp real(wp),parameter :: a2010 = 13.83320427252914990351082875460544773493_wp real(wp),parameter :: a2011 = -12.28092433078461863729523583784519048012_wp real(wp),parameter :: a2012 = .1719151595656509762746810113378644307112_wp real(wp),parameter :: a2014 = .3641583114314496380113822384214528216140e-1_wp real(wp),parameter :: a2015 = .2961920580288763054890146412520723429115e-1_wp real(wp),parameter :: a2016 = -.2651793938627067002647615623738425030047_wp real(wp),parameter :: a2017 = .942950396173806655317007970358739475630e-1_wp real(wp),parameter :: c21 = .9_wp real(wp),parameter :: a2101 = .3939058345528250943410670634923521987132e-1_wp real(wp),parameter :: a2106 = .3558516141234424183136697322755323715063_wp real(wp),parameter :: a2107 = .4197382225952610029372225526720065366258_wp real(wp),parameter :: a2108 = .872044977807194166293172525204036071060_wp real(wp),parameter :: a2109 = .898952083487659486126627160171417043611_wp real(wp),parameter :: a2110 = -.630580616105988359023456649527853470403_wp real(wp),parameter :: a2111 = -1.121887220595483550736681645425215081433_wp real(wp),parameter :: a2112 = .4298219512400197176967511031829197714867e-1_wp real(wp),parameter :: a2114 = .1332557566873915707013495891889190564164e-1_wp real(wp),parameter :: a2115 = .1876227053964148034446101291928097773800e-1_wp real(wp),parameter :: a2116 = -.1859411132922105570515379368592596513699_wp real(wp),parameter :: a2117 = .1773614271924602745226064729836361000042_wp associate (t => me%dense_t, h => me%dense_h, x => me%dense_x, k => me%funcs) ! f at the end of the step (this is also the first stage of the next step, ! so it is saved in the function evaluation cache): call me%set_fsal_cache(t+h, me%dense_xf, k(:,14)) if (me%stopped) return call me%f(t+c15*h, x+h*(a1501*k(:,1) + a1506*k(:,6) + a1507*k(:,7) + a1508*k(:,8) + & a1509*k(:,9) + a1510*k(:,10) + a1511*k(:,11) + a1512*k(:,12) + & a1514*k(:,14)), k(:,15)) if (me%stopped) return call me%f(t+c16*h, x+h*(a1601*k(:,1) + a1606*k(:,6) + a1607*k(:,7) + a1608*k(:,8) + & a1609*k(:,9) + a1610*k(:,10) + a1611*k(:,11) + a1612*k(:,12) + & a1614*k(:,14) + a1615*k(:,15)), k(:,16)) if (me%stopped) return call me%f(t+c17*h, x+h*(a1701*k(:,1) + a1706*k(:,6) + a1707*k(:,7) + a1708*k(:,8) + & a1709*k(:,9) + a1710*k(:,10) + a1711*k(:,11) + a1712*k(:,12) + & a1714*k(:,14) + a1715*k(:,15) + a1716*k(:,16)), k(:,17)) if (me%stopped) return call me%f(t+c18*h, x+h*(a1801*k(:,1) + a1806*k(:,6) + a1807*k(:,7) + a1808*k(:,8) + & a1809*k(:,9) + a1810*k(:,10) + a1811*k(:,11) + a1812*k(:,12) + & a1814*k(:,14) + a1815*k(:,15) + a1816*k(:,16) + a1817*k(:,17)), k(:,18)) if (me%stopped) return call me%f(t+c19*h, x+h*(a1901*k(:,1) + a1906*k(:,6) + a1907*k(:,7) + a1908*k(:,8) + & a1909*k(:,9) + a1910*k(:,10) + a1911*k(:,11) + a1912*k(:,12) + & a1914*k(:,14) + a1915*k(:,15) + a1916*k(:,16) + a1917*k(:,17)), k(:,19)) if (me%stopped) return call me%f(t+c20*h, x+h*(a2001*k(:,1) + a2006*k(:,6) + a2007*k(:,7) + a2008*k(:,8) + & a2009*k(:,9) + a2010*k(:,10) + a2011*k(:,11) + a2012*k(:,12) + & a2014*k(:,14) + a2015*k(:,15) + a2016*k(:,16) + a2017*k(:,17)), k(:,20)) if (me%stopped) return call me%f(t+c21*h, x+h*(a2101*k(:,1) + a2106*k(:,6) + a2107*k(:,7) + a2108*k(:,8) + & a2109*k(:,9) + a2110*k(:,10) + a2111*k(:,11) + a2112*k(:,12) + & a2114*k(:,14) + a2115*k(:,15) + a2116*k(:,16) + a2117*k(:,17)), k(:,21)) end associate end procedure rkv87e_dense_extra_stages !***************************************************************************************** !***************************************************************************************** !> ! Dense output for Verner's "most robust" 8(7) method [[rkv87r]]. ! This is the interpolant of order 8 (`bi8`) from the coefficient file, which requires 8 additional ! function evaluations (see [[rkv87r_dense_extra_stages]]). ! !@note The interpolant is published in the polynomial form. Here it has been ! converted (exactly, using high precision arithmetic) ! to the form used in `DOPRI5` and `DOP853` (see [[dense_coefficients]]), which has ! much smaller roundoff error for these coefficients. This code is generated by ! `scripts/generate_dense.py`. ! !### References ! * J.H. Verner, "Numerically optimal Runge--Kutta pairs with interpolants", ! Numerical Algorithms, 53, (2010) pp. 383-396. [doi](https://doi.org/10.1007/s11075-009-9290-3) ! * Coefficients: ! * [RKV87.IIa.Robust.00000754677.081208.CoeffsOnlyFLOAT](https://web.archive.org/web/https://www.sfu.ca/~jverner/RKV87.IIa.Robust.00000754677.081208.CoeffsOnlyFLOAT) module procedure rkv87r_dense_coefficients form = dense_form_dopri allocate(b(21,8)) b = 0.0_wp b(1,:) = [ 4.472956466669571420301584042904938255000e-2_wp, & 9.552704353333042857969841595709506174500e-1_wp, & -9.105408706666085715939683191419012354200e-1_wp, & -9.691232641845086332345535507614873258800e0_wp, & 8.020974285630397028779403348796943323900e0_wp, & 3.787005946021333916730459065720655049316e1_wp, & -1.641062691938850782359503419698494761476e1_wp, & -5.713436418098746726741860275346451549916e1_wp ] b(6,:) = [ 1.569103352770819981336869801072664509900e-1_wp, & -1.569103352770819981336869801072664509900e-1_wp, & 3.138206705541639962673739602145329174000e-1_wp, & -4.335995671088398387805626043337848836650e1_wp, & 3.120422366997099655227140957766452403900e1_wp, & 3.079372796893748326011330997461193143516e2_wp, & -5.441475104486891781041242393781865817660e1_wp, & -7.231283387936238598986478466688726493811e2_wp ] b(7,:) = [ 1.846097340815163774070245187352627800500e-1_wp, & -1.846097340815163774070245187352627800500e-1_wp, & 3.692194681630327548140490374705255548000e-1_wp, & -5.101429465463324049304001418153526419300e1_wp, & 3.671270871839687325892758634464185699500e1_wp, & 3.622974816595408571692416676501795405396e2_wp, & -6.402059305249850189817638807278158619260e1_wp, & -8.507822642515111225480921303415510063376e2_wp ] b(8,:) = [ 2.251638060208699104247941940035072190900e-1_wp, & -2.251638060208699104247941940035072190900e-1_wp, & 4.503276120417398208495883880070144394000e-1_wp, & -6.222084010388814188347070939275397438800e1_wp, & 4.477755880802963922110241162454513172000e1_wp, & 4.418850408300662187361683656495595096820e2_wp, & -7.808429207231659798021361520992182790000e1_wp, & -1.037677529123876112401217700900292618796e3_wp ] b(9,:) = [ 1.479461565197023468700517988544914183600e-1_wp, & -1.479461565197023468700517988544914183600e-1_wp, & 2.958923130394046937401035977089828308000e-1_wp, & -4.088283242087312568043685586449240575050e1_wp, & 2.942154798790754312177190829253203567500e1_wp, & 2.903450362191007028698435935429060761310e2_wp, & -5.130607401258065144337214398038054267500e1_wp, & -6.818165177333586028419626797962516386195e2_wp ] b(10,:) = [ 7.605554244495582526979836191033649145000e-2_wp, & -7.605554244495582526979836191033649145000e-2_wp, & 1.521110848899116505395967238206729804000e-1_wp, & -2.101687579860616621660785317536682386570e1_wp, & 1.512490654998942570342542625068576530000e1_wp, & 1.492593639828912361803202355012922501018e2_wp, & -2.637521231738263633777869820748419663180e1_wp, & -3.505053887441520158806914746272446155388e2_wp ] b(11,:) = [ 1.227729023501861961082434631592143749600e-1_wp, & -1.227729023501861961082434631592143749600e-1_wp, & 2.455458047003723922164869263184287442000e-1_wp, & -3.392655889602951431669989043738191700250e1_wp, & 2.441542871463273670603836349008726242600e1_wp, & 2.409424051164302742311458671525281903098e2_wp, & -4.257627073859899515498257965481217790280e1_wp, & -5.658044434650382773588442714671675697003e2_wp ] b(12,:) = [ 4.181195863899163158338484280087188227000e-2_wp, & -4.181195863899163158338484280087188227000e-2_wp, & 8.362391727798326316676968560174376440000e-2_wp, & -1.155414468640643246615618317240418312970e1_wp, & 8.315001731063372403744186453126939249200e0_wp, & 8.205616780462199246812245029287127400990e1_wp, & -1.449992007232295459021804971750014848870e1_wp, & -1.926922923141433811127686430865745560445e2_wp ] b(14,:) = [ -2.200000000000000000000000000000000000000e-38_wp, & 2.200000000000000000000000000000000000000e-38_wp, & -9.999999999999999999999999999999999998900e-1_wp, & 1.280089315230128997622478115178850386850e0_wp, & 1.214526658294784720699270859503774880420e1_wp, & -1.633768081413412302265308771151546083612e1_wp, & -3.126073065386297197044006830249323867748e1_wp, & 4.159290247459829082789858665770259971247e1_wp ] b(15,:) = [ -1.290000000000000000000000000000000000000e-35_wp, & 1.290000000000000000000000000000000000000e-35_wp, & 4.600000000000000000000000000000000000000e-35_wp, & -3.787724238112221117869223214531251760350e2_wp, & 2.923048945259776813119373129404403050400e2_wp, & 2.444376034006293946505843863551276680274e3_wp, & -6.303204814666199773256223628619521683540e2_wp, & -5.205692776382276693190581775232023008049e3_wp ] b(16,:) = [ -3.800000000000000000000000000000000000000e-37_wp, & 3.800000000000000000000000000000000000000e-37_wp, & 2.200000000000000000000000000000000000000e-36_wp, & 4.050201956693941882963694376959127659450e1_wp, & -3.125607307543530818046684278491619417100e1_wp, & -2.613763825838144257157743234144348393888e2_wp, & 6.739997652661045583621718069646876732080e1_wp, & 5.566431382913769754094230327362926809033e2_wp ] b(17,:) = [ -1.460000000000000000000000000000000000000e-36_wp, & 1.460000000000000000000000000000000000000e-36_wp, & 2.400000000000000000000000000000000000000e-36_wp, & 8.287874270864547488242848049989674796000e1_wp, & -6.395888566051017670305033696122239160500e1_wp, & -5.348510072807066900977094370708424840840e2_wp, & 1.379196734593797153737594683090522991920e2_wp, & 1.139051433293012791190246990806004069021e3_wp ] b(18,:) = [ -8.900000000000000000000000000000000000000e-36_wp, & 8.900000000000000000000000000000000000000e-36_wp, & 4.800000000000000000000000000000000000000e-35_wp, & 6.696598644458484065321915951171526447600e2_wp, & -5.091332373841825361638643865010116060100e2_wp, & -4.002367935374949534827244213737313659466e3_wp, & 6.824515520353815526691641610864065213660e2_wp, & 8.955559150961236338683195435020456911796e3_wp ] b(19,:) = [ -3.200000000000000000000000000000000000000e-36_wp, & 3.200000000000000000000000000000000000000e-36_wp, & 1.200000000000000000000000000000000000000e-35_wp, & -2.440795068338021105727923077508888854030e2_wp, & 1.877569738326881532377272300727173370300e2_wp, & 1.185275785880848271417648290848137881656e3_wp, & -1.991885609601217401712182030748473548600e1_wp, & -2.761779565981876362182063486724758808276e3_wp ] b(20,:) = [ -9.000000000000000000000000000000000000000e-37_wp, & 9.000000000000000000000000000000000000000e-37_wp, & 3.400000000000000000000000000000000000000e-36_wp, & 6.811738686016824405368089974720930253800e1_wp, & -4.410018495578055594559885870373932497900e1_wp, & -5.047663575101844172639638093991741964740e2_wp, & 4.589112206770594183809840002291297460100e1_wp, & 1.225463681207094529961063380448631794824e3_wp ] b(21,:) = [ 2.640000000000000000000000000000000000000e-36_wp, & -2.640000000000000000000000000000000000000e-36_wp, & -9.700000000000000000000000000000000000000e-36_wp, & 3.408056366135824033096753411991316907850e1_wp, & -4.175110433132608875973752203938633291200e1_wp, & -2.225452910855924804194271532587966272754e2_wp, & 9.552548435737522063469397433477366562440e1_wp, & 5.087031747435249686104611859291129299849e2_wp ] end procedure rkv87r_dense_coefficients !***************************************************************************************** !***************************************************************************************** !> ! The additional function evaluations for the dense output of [[rkv87r]]. ! These are stored in `me%funcs(:,14:21)`. Stage 14 is `f` at the end of the step, ! which is cached for the first stage of the next step (see [[set_fsal_cache]]). module procedure rkv87r_dense_extra_stages real(wp),parameter :: c15 = .3110177634953863863927417318829099695921_wp real(wp),parameter :: a1501 = .4584909962591746145269856107990713535990e-1_wp real(wp),parameter :: a1506 = .781870674437554054886342907851634313370e-2_wp real(wp),parameter :: a1507 = .6155870622618810523423270651592335069895e-2_wp real(wp),parameter :: a1508 = .2195571449623681030222448202857462664973_wp real(wp),parameter :: a1509 = -.1222581732697346808194498111979519895319e-1_wp real(wp),parameter :: a1510 = .3261300492958830346952381858768844254700e-1_wp real(wp),parameter :: a1511 = .1254797805520437910933978712641756586385e-1_wp real(wp),parameter :: a1512 = .1129681142614891903100335083803550169907e-1_wp real(wp),parameter :: a1514 = -.1259503554386166268241032464519842162533e-1_wp real(wp),parameter :: c16 = .536_wp real(wp),parameter :: a1601 = .4167847284138612914768499355650240017377e-1_wp real(wp),parameter :: a1606 = -.4073849103699890037905839225012252772752_wp real(wp),parameter :: a1607 = -.2302256719981637874461745430790814930259e-1_wp real(wp),parameter :: a1608 = .2482392242653016164318596667033368515282_wp real(wp),parameter :: a1609 = .7400026817560596139676184623176512480681e-1_wp real(wp),parameter :: a1610 = -.1930141806259700048703973978096116667950_wp real(wp),parameter :: a1611 = .6690280284382524116505423177469790757302_wp real(wp),parameter :: a1612 = -.1099229564510122581232500661965260824893e-1_wp real(wp),parameter :: a1614 = .1480389368968139520045852342144953350362e-1_wp real(wp),parameter :: a1615 = .1226640664306490993906164335783647158793_wp real(wp),parameter :: c17 = .132_wp real(wp),parameter :: a1701 = .5646609456796632256456589489425447704092e-1_wp real(wp),parameter :: a1706 = -.3522910440496784686626294493945401506422e-1_wp real(wp),parameter :: a1707 = .6568161065081292073153281369860458994523e-1_wp real(wp),parameter :: a1708 = .7383545165330150178251163230257161341400e-1_wp real(wp),parameter :: a1709 = -.1925857315058937001203418153699886798201_wp real(wp),parameter :: a1710 = .3765563615652640794139928459441174753686e-1_wp real(wp),parameter :: a1711 = -.5186009042220349899736698512682662022351_wp real(wp),parameter :: a1712 = .9190279569743070828256110335276636243207e-1_wp real(wp),parameter :: a1714 = -.1036134627781298571789386786962991910815_wp real(wp),parameter :: a1715 = .1051965910775021413022577292180480274178_wp real(wp),parameter :: a1716 = .5512910231074863915343848322133512704141_wp real(wp),parameter :: c18 = .334_wp real(wp),parameter :: a1801 = .3934466647027139074996525141009606759203e-1_wp real(wp),parameter :: a1806 = -.1291667299515080090913774283193174596360e-1_wp real(wp),parameter :: a1807 = .1945891642767310407515289952358063928465e-2_wp real(wp),parameter :: a1808 = -.3950219424001926036159581502570781977097e-1_wp real(wp),parameter :: a1809 = -.4501175554198611234137687445953085443754e-2_wp real(wp),parameter :: a1810 = -.998931335985195283267955477426531672585e-2_wp real(wp),parameter :: a1811 = -.1198539244102761951947005140289460975181e-1_wp real(wp),parameter :: a1812 = .1883152441766023211187579107594459282194e-2_wp real(wp),parameter :: a1814 = -.2288623796054366081491763462894660549547e-2_wp real(wp),parameter :: a1815 = .1239520837697704951665852820403595312797_wp real(wp),parameter :: a1816 = .1852238959793020024702401439703066043180e-1_wp real(wp),parameter :: a1817 = .2295351884637971911562351980362084556913_wp real(wp),parameter :: c19 = .35_wp real(wp),parameter :: a1901 = .3884096081168368537927354267981728648614e-1_wp real(wp),parameter :: a1906 = -.2427245527708823693410501445101278168941e-1_wp real(wp),parameter :: a1907 = .3128428375929178417630603820449548457213e-2_wp real(wp),parameter :: a1908 = -.7358478051063769623956621836398499737217e-1_wp real(wp),parameter :: a1909 = -.8694978890870701525550785891871065419813e-2_wp real(wp),parameter :: a1910 = -.1865660321520561173012739484485694591980e-1_wp real(wp),parameter :: a1911 = -.2246454061884046921515692264772804996290e-1_wp real(wp),parameter :: a1912 = .3374661498794678376921336894109603882993e-2_wp real(wp),parameter :: a1914 = -.3968398066426385004345096186031114335021e-2_wp real(wp),parameter :: a1915 = .1587081125691119989004468338247673365291_wp real(wp),parameter :: a1916 = .3991436969512407771883001041914220455046e-1_wp real(wp),parameter :: a1917 = .2576752236284254818557491047471989747932_wp real(wp),parameter :: c20 = .755_wp real(wp),parameter :: a2001 = .7198494848014507493654612298114036727341e-1_wp real(wp),parameter :: a2006 = .5076533962269455568306708785536594617654_wp real(wp),parameter :: a2007 = -.1572363214423407410943664004339326120598e-1_wp real(wp),parameter :: a2008 = 1.478215037037996550905747792359144808839_wp real(wp),parameter :: a2009 = .2041154554617372630941975148106197833715_wp real(wp),parameter :: a2010 = .3793878839474901527036840682829095687643_wp real(wp),parameter :: a2011 = .4643932887538339042346694891381803095294_wp real(wp),parameter :: a2012 = -.5513989032591925973677454467001044028380e-1_wp real(wp),parameter :: a2014 = .5665552314958980596114050893400471586229e-1_wp real(wp),parameter :: a2015 = -.6894099744391873790684724028801677847224_wp real(wp),parameter :: a2016 = -.5600773486771496063800966221681643192010_wp real(wp),parameter :: a2017 = -1.087054687471247989371876165297923209992_wp real(wp),parameter :: c21 = .95_wp real(wp),parameter :: a2101 = .5696655426838587836617455691105627859732e-1_wp real(wp),parameter :: a2106 = .3386438327349868685266036706413028229066_wp real(wp),parameter :: a2107 = .1470519205552199854380086863315794888932_wp real(wp),parameter :: a2108 = .7933978111737925583755765349543256757167_wp real(wp),parameter :: a2109 = .2067171591185578546324388929571733817983_wp real(wp),parameter :: a2110 = .2188160879075981870381792015062558329913_wp real(wp),parameter :: a2111 = .2925192926438936110272980692420782663877_wp real(wp),parameter :: a2112 = .1215446220886579630956353870439258394935e-1_wp real(wp),parameter :: a2114 = -.229214343240334151600210819062643523147e-2_wp real(wp),parameter :: a2115 = -.3225954701993023482175429441563279152604_wp real(wp),parameter :: a2116 = -.2993120292318899024302980359438239792811_wp real(wp),parameter :: a2117 = -.4920674777477051475500000629573860014675_wp associate (t => me%dense_t, h => me%dense_h, x => me%dense_x, k => me%funcs) ! f at the end of the step (this is also the first stage of the next step, ! so it is saved in the function evaluation cache): call me%set_fsal_cache(t+h, me%dense_xf, k(:,14)) if (me%stopped) return call me%f(t+c15*h, x+h*(a1501*k(:,1) + a1506*k(:,6) + a1507*k(:,7) + a1508*k(:,8) + & a1509*k(:,9) + a1510*k(:,10) + a1511*k(:,11) + a1512*k(:,12) + & a1514*k(:,14)), k(:,15)) if (me%stopped) return call me%f(t+c16*h, x+h*(a1601*k(:,1) + a1606*k(:,6) + a1607*k(:,7) + a1608*k(:,8) + & a1609*k(:,9) + a1610*k(:,10) + a1611*k(:,11) + a1612*k(:,12) + & a1614*k(:,14) + a1615*k(:,15)), k(:,16)) if (me%stopped) return call me%f(t+c17*h, x+h*(a1701*k(:,1) + a1706*k(:,6) + a1707*k(:,7) + a1708*k(:,8) + & a1709*k(:,9) + a1710*k(:,10) + a1711*k(:,11) + a1712*k(:,12) + & a1714*k(:,14) + a1715*k(:,15) + a1716*k(:,16)), k(:,17)) if (me%stopped) return call me%f(t+c18*h, x+h*(a1801*k(:,1) + a1806*k(:,6) + a1807*k(:,7) + a1808*k(:,8) + & a1809*k(:,9) + a1810*k(:,10) + a1811*k(:,11) + a1812*k(:,12) + & a1814*k(:,14) + a1815*k(:,15) + a1816*k(:,16) + a1817*k(:,17)), k(:,18)) if (me%stopped) return call me%f(t+c19*h, x+h*(a1901*k(:,1) + a1906*k(:,6) + a1907*k(:,7) + a1908*k(:,8) + & a1909*k(:,9) + a1910*k(:,10) + a1911*k(:,11) + a1912*k(:,12) + & a1914*k(:,14) + a1915*k(:,15) + a1916*k(:,16) + a1917*k(:,17)), k(:,19)) if (me%stopped) return call me%f(t+c20*h, x+h*(a2001*k(:,1) + a2006*k(:,6) + a2007*k(:,7) + a2008*k(:,8) + & a2009*k(:,9) + a2010*k(:,10) + a2011*k(:,11) + a2012*k(:,12) + & a2014*k(:,14) + a2015*k(:,15) + a2016*k(:,16) + a2017*k(:,17)), k(:,20)) if (me%stopped) return call me%f(t+c21*h, x+h*(a2101*k(:,1) + a2106*k(:,6) + a2107*k(:,7) + a2108*k(:,8) + & a2109*k(:,9) + a2110*k(:,10) + a2111*k(:,11) + a2112*k(:,12) + & a2114*k(:,14) + a2115*k(:,15) + a2116*k(:,16) + a2117*k(:,17)), k(:,21)) end associate end procedure rkv87r_dense_extra_stages !***************************************************************************************** !***************************************************************************************** !> ! Dense output for Verner's "most efficient" 9(8) method [[rkv98e]]. ! This is the interpolant of order 9 (`bi9`) from the coefficient file, which requires 10 additional ! function evaluations (see [[rkv98e_dense_extra_stages]]). ! !@note The interpolant is published in the polynomial form. Here it has been ! converted (exactly, using high precision arithmetic) ! to the form used in `DOPRI5` and `DOP853` (see [[dense_coefficients]]), which has ! much smaller roundoff error for these coefficients. This code is generated by ! `scripts/generate_dense.py`. ! !### References ! * J.H. Verner, "Numerically optimal Runge--Kutta pairs with interpolants", ! Numerical Algorithms, 53, (2010) pp. 383-396. [doi](https://doi.org/10.1007/s11075-009-9290-3) ! * Coefficients: ! * [RKV98.IIa.Efficient.000000349.081209.FLOAT6040OnWeb](https://web.archive.org/web/https://www.sfu.ca/~jverner/RKV98.IIa.Efficient.000000349.081209.FLOAT6040OnWeb) ! * [RKV98.IIa.Efficient.000000349.081210.CoeffsOnlyRADandFLOATS](https://web.archive.org/web/https://www.sfu.ca/~jverner/RKV98.IIa.Efficient.000000349.081210.CoeffsOnlyRADandFLOATS) module procedure rkv98e_dense_coefficients form = dense_form_dopri allocate(b(26,9)) b = 0.0_wp b(1,:) = [ 1.461197685842315252051541915018784740000e-2_wp, & 9.853880231415768474794845808498121526000e-1_wp, & -9.707760462831536949589691616996243076000e-1_wp, & -2.637432463119266769046921927403938348680e1_wp, & 2.662373586934323228524553986598284507950e1_wp, & 1.773103840298460763640663504147635792890e2_wp, & -1.803659001326322052043358708639748570580e2_wp, & -3.602356071514077449712577435776812132872e2_wp, & 3.687958613829998340610814211036270246107e2_wp ] b(8,:) = [ -3.915211862331339089410228267288242031190e-1_wp, & 3.915211862331339089410228267288242031190e-1_wp, & -7.830423724662678178820456534576484053200e-1_wp, & 3.824219802469492939508450383737897284700e0_wp, & -1.089055975957825383000248358158210906500e1_wp, & -7.854717108610908774516648907685957990100e1_wp, & 2.023877746331953159677827195710041818110e2_wp, & 4.056749190361309957310643884604588000482e2_wp, & -9.192934944890586676320942978986329899642e2_wp ] b(9,:) = [ 2.310932500289506415909675644868993670340e-1_wp, & -2.310932500289506415909675644868993670340e-1_wp, & 4.621865000579012831819351289737987333600e-1_wp, & -2.257224932015584705420024199011715147000e0_wp, & 6.428093645938330909398683444200459790500e0_wp, & 4.636204038281700096488870120648811430100e1_wp, & -1.194582828482291707753994953095320113990e2_wp, & -2.394474138098560440235722472353826804126e2_wp, & 5.426079835318221405169412532400889768401e2_wp ] b(10,:) = [ 1.274766769992852382560589467488989178900e-1_wp, & -1.274766769992852382560589467488989178900e-1_wp, & 2.549533539985704765121178934977978334400e-1_wp, & -1.245140364494578097833513717320252273000e0_wp, & 3.545893345313125885208623162493497216500e0_wp, & 2.557443303154802282059111897138225875000e1_wp, & -6.589610443241086147117209679198861537100e1_wp, & -1.320850376405596611114884974599119629404e2_wp, & 2.993158071265785138462868993727082901209e2_wp ] b(11,:) = [ 2.246434176204157731566981937082069679690e-1_wp, & -2.246434176204157731566981937082069679690e-1_wp, & 4.492868352408315463133963874164139412200e-1_wp, & -2.194225590762462158525288157285090730100e0_wp, & 6.248685001516756011431575754085912984500e0_wp, & 4.506807186340142222223740477328845454400e1_wp, & -1.161241919386681626333128081302600267300e2_wp, & -2.327644159744056704155810877594301364774e2_wp, & 5.274637482204278644179968961638568925209e2_wp ] b(12,:) = [ 5.684352689748512932705226972873692149210e-1_wp, & -5.684352689748512932705226972873692149210e-1_wp, & 1.136870537949702586541045394574738414900e0_wp, & -5.552244651050077280327639912126807771000e0_wp, & 1.581160479662098110835663977563675363900e1_wp, & 1.140397605379126189528653276941578598950e2_wp, & -2.938393965795200616457824748216184383800e2_wp, & -5.889845551840504285132965527892994927790e2_wp, & 1.334688551172190921099749059976639173619e3_wp ] b(13,:) = [ 5.825871557215827200814768021863420863830e-2_wp, & -5.825871557215827200814768021863420863830e-2_wp, & 1.165174311443165440162953604372684198100e-1_wp, & -5.690474528364889614553757119684679897000e-1_wp, & 1.620525391126758704572308723626127960200e0_wp, & 1.168789189502107698674012676467766010480e1_wp, & -3.011548854119085473634893984513174406880e1_wp, & -6.036480413811171971813297190679316339960e1_wp, & 1.367917244480423273434147193694336909663e2_wp ] b(14,:) = [ 1.364317403482215641609022744494239841501e-1_wp, & -1.364317403482215641609022744494239841501e-1_wp, & 2.728634806964431283218045488988479684000e-1_wp, & -1.332609783252874831920549764292705279400e0_wp, & 3.794987534801833972494696485716918801000e0_wp, & 2.737100220248684757459214589310884501500e1_wp, & -7.052521623176904359940982517912716128800e1_wp, & -1.413638320628874645660359292753972443304e2_wp, & 3.203423358787481875842587982911148385364e2_wp ] b(15,:) = [ 3.057013983082797397721005067920369640230e-2_wp, & -3.057013983082797397721005067920369640230e-2_wp, & 6.114027966165594795442010135840739318600e-2_wp, & -2.985966997854893020569951149206441786900e-1_wp, & 8.503395126312590577968392073601393201000e-1_wp, & 6.132996343111081058218285676277561234300e0_wp, & -1.580252305146719597804528565242054186030e1_wp, & -3.167526927498059401319033134641697060458e1_wp, & 7.177882490212657594681492031347005327988e1_wp ] b(17,:) = [ -9.050000000000000000000000000000000000000e-38_wp, & 9.050000000000000000000000000000000000000e-38_wp, & -9.999999999999999999999999999999999998200e-1_wp, & 5.527558093255900558295661824035176973000e-1_wp, & 1.813638202994577983855179752567625208370e1_wp, & -1.584023651189561810655026292535406483000e0_wp, & -1.138727071799036179942913514628938258960e2_wp, & -9.567966115228650651917704209071222150400e0_wp, & 2.486243327770222987362193390543305737239e2_wp ] b(18,:) = [ -4.790000000000000000000000000000000000000e-38_wp, & 4.790000000000000000000000000000000000000e-38_wp, & 1.700000000000000000000000000000000000000e-37_wp, & -1.507876007899797948720901584434839157500e-1_wp, & 8.741438023337320473924434808369258824500e-1_wp, & 4.351609820910940598179738328896931649500e0_wp, & -1.956907482678494674590176993148583098680e1_wp, & -1.681535793587097508424107764872161423792e1_wp, & 6.970457078142274038253812108643441743987e1_wp ] b(19,:) = [ 5.915000000000000000000000000000000000000e-37_wp, & -5.915000000000000000000000000000000000000e-37_wp, & -2.910000000000000000000000000000000000000e-36_wp, & -6.413652207435296452288504944964177514000e-1_wp, & 3.718113623455295804540722809450322061300e0_wp, & 1.850928842130449090763628813845471638000e1_wp, & -8.323578285132855661128262757472927635200e1_wp, & -7.152302774180107756124674782038207395020e1_wp, & 2.964838434449778148523985255750527153802e2_wp ] b(20,:) = [ -5.950000000000000000000000000000000000000e-37_wp, & 5.950000000000000000000000000000000000000e-37_wp, & 1.540000000000000000000000000000000000000e-36_wp, & 1.810729313444845732964058528284532355100e0_wp, & -1.049713503470644128498706300479128801500e1_wp, & -5.225620291135724096128737190590253896700e1_wp, & 2.349947690673105785282725197801339672630e2_wp, & 2.019268253558706214778963156521776929512e2_wp, & -8.370456783831301740195014698000522807852e2_wp ] b(21,:) = [ 3.087000000000000000000000000000000000000e-38_wp, & -3.087000000000000000000000000000000000000e-38_wp, & -1.510000000000000000000000000000000000000e-37_wp, & 5.176008760353717918864555990277480376000e-2_wp, & -3.000628669057631822942504016123986812800e-1_wp, & -1.493754820467511215734523897132065955300e0_wp, & 6.717376110820595274696522515117143918900e0_wp, & 5.772121814298099976139910906751945603100e0_wp, & -2.392713108446217690295957956014097092250e1_wp ] b(22,:) = [ -1.000000000000000000000000000000000000000e-37_wp, & 1.000000000000000000000000000000000000000e-37_wp, & 0.0_wp, & 3.132178255668799909977422939838912846080e1_wp, & -3.392420649401873037663618247998439760100e1_wp, & -2.590893142265158773925207906443653480920e2_wp, & 2.866823287602136863278753834362081748040e2_wp, & 5.895648155050292008953915349445446995122e2_wp, & -6.785938956640391428503413103061359428182e2_wp ] b(23,:) = [ 7.330000000000000000000000000000000000000e-37_wp, & -7.330000000000000000000000000000000000000e-37_wp, & -2.610000000000000000000000000000000000000e-36_wp, & -2.719607334185924760747802644504744091040e0_wp, & 1.500803987485908225866556647903159336750e0_wp, & 7.970043750272488226418632681947221879890e1_wp, & -5.099535355799443025608158655435623819350e1_wp, & -1.583747217185762625252436285453161829049e2_wp, & -1.233784294340547057337599296127606178639e1_wp ] b(24,:) = [ -4.970000000000000000000000000000000000000e-37_wp, & 4.970000000000000000000000000000000000000e-37_wp, & 1.300000000000000000000000000000000000000e-36_wp, & 3.177214801432923432265738869490200555600e0_wp, & -7.520518797062100203150956724560663824000e0_wp, & -9.993492137158635125376097072828656105200e1_wp, & 1.939160363419817684647692048571334260330e2_wp, & 4.661545453302695554186806153605924696298e2_wp, & -8.214378043648253954175634277881875971878e2_wp ] b(25,:) = [ 5.408000000000000000000000000000000000000e-37_wp, & -5.408000000000000000000000000000000000000e-37_wp, & -2.470000000000000000000000000000000000000e-36_wp, & 8.772843083465530069477626269697233862000e-1_wp, & 2.953487199177759004908979283138970972500e0_wp, & -3.271702282236840403462928162309868480400e1_wp, & -3.690270909644872824643482573797112761000e1_wp, & 2.169042458546745945296385613610223960816e2_wp, & -1.411134269128985501802080532710905715931e2_wp ] b(26,:) = [ -1.124000000000000000000000000000000000000e-36_wp, & 1.124000000000000000000000000000000000000e-36_wp, & 5.650000000000000000000000000000000000000e-36_wp, & 1.719427581798715782378897599231938077600e0_wp, & -2.897431278741946397869446997357742792750e1_wp, & -3.048550514149042630044736051278801471500e1_wp, & 2.720044463548258913344026076958928013620e2_wp, & 1.572045358514632251263931928882559536476e2_wp, & -7.828483098245396412116558219612402319811e2_wp ] end procedure rkv98e_dense_coefficients !***************************************************************************************** !***************************************************************************************** !> ! The additional function evaluations for the dense output of [[rkv98e]]. ! These are stored in `me%funcs(:,17:26)`. Stage 17 is `f` at the end of the step, ! which is cached for the first stage of the next step (see [[set_fsal_cache]]). module procedure rkv98e_dense_extra_stages real(wp),parameter :: c18 = .7404185470631561083004100761798676215811_wp real(wp),parameter :: a1801 = .1549973668189559302279946863304789372788e-1_wp real(wp),parameter :: a1808 = .3355153219059635054403439303177105512242_wp real(wp),parameter :: a1809 = .2003613944191860651552622660712101217322_wp real(wp),parameter :: a1810 = .1252060659283549312946162355194540994211_wp real(wp),parameter :: a1811 = .2298676393184206750544046308957155868736_wp real(wp),parameter :: a1812 = -.2020250653476181447824906889122391003637_wp real(wp),parameter :: a1813 = .5917103230665456601422111997583025339897e-1_wp real(wp),parameter :: a1814 = -.2651834783047638681693835956996437528251e-1_wp real(wp),parameter :: a1815 = -.2384094602130971415278110567256446033405e-1_wp real(wp),parameter :: a1817 = .2718171570208501807097257892166705118335e-1_wp real(wp),parameter :: c19 = .888_wp real(wp),parameter :: a1901 = .1302453943114338366054520296881099431474e-1_wp real(wp),parameter :: a1908 = -.7452850902413112085299330666038981625179_wp real(wp),parameter :: a1909 = .2643867896429300961465132150322749722129_wp real(wp),parameter :: a1910 = .1313961758372753932588328082078842388890_wp real(wp),parameter :: a1911 = .2167253815122927263092467187957410643315_wp real(wp),parameter :: a1912 = .8734117564076052559016338094938888451419_wp real(wp),parameter :: a1913 = .1185905643935776688228545787724340848142e-1_wp real(wp),parameter :: a1914 = .5876002941689550612992712203494447529933e-1_wp real(wp),parameter :: a1915 = .3266518630202087866399279690939423159022e-2_wp real(wp),parameter :: a1917 = -.8959308648417929824525368306101792182274e-2_wp real(wp),parameter :: a1918 = .6941415157202692219907482080827253287034e-1_wp real(wp),parameter :: c20 = .696_wp real(wp),parameter :: a2001 = .1397089996925942721283716334050740168797e-1_wp real(wp),parameter :: a2008 = -.4665765335957674596054673402956853940520_wp real(wp),parameter :: a2009 = .2416372787216257077935214889875485248580_wp real(wp),parameter :: a2010 = .1290363341345674735721677437066933999929_wp real(wp),parameter :: a2011 = .2216700671735105311080225734522323922813_wp real(wp),parameter :: a2012 = .6257275123364644931771253383573999863003_wp real(wp),parameter :: a2013 = .4355312415679284117869124964829805160429e-1_wp real(wp),parameter :: a2014 = .1011962491667290833450024852274278874501_wp real(wp),parameter :: a2015 = .1808582254679721049279369742685497400353e-1_wp real(wp),parameter :: a2017 = -.2079875587689169691156509689282083267654e-1_wp real(wp),parameter :: a2018 = -.9022232517086218976198252891464664868640e-1_wp real(wp),parameter :: a2019 = -.1212796735622254216011467740438097427634_wp real(wp),parameter :: c21 = .487_wp real(wp),parameter :: a2101 = .1604638888318112738641232352800290501904e-1_wp real(wp),parameter :: a2108 = .9517712399458336651642257453589397190702e-1_wp real(wp),parameter :: a2109 = .1359187264655317806136927180199100622471_wp real(wp),parameter :: a2110 = .1237765280959854006935081364365637515893_wp real(wp),parameter :: a2111 = .2335656264102966047058755123098072346246_wp real(wp),parameter :: a2112 = -.9051508172625873314662090873741762206189e-1_wp real(wp),parameter :: a2113 = -.2537574270006131028513276914038326155331e-1_wp real(wp),parameter :: a2114 = -.1359631696887162048002744757083947500478_wp real(wp),parameter :: a2115 = -.4679214284145113075088049469061349990847e-1_wp real(wp),parameter :: a2117 = .5177958859391748239949773879090325427473e-1_wp real(wp),parameter :: a2118 = .9672595677476773313884172931875718705561e-1_wp real(wp),parameter :: a2119 = .1477312690340742769720989417101989769314_wp real(wp),parameter :: a2120 = -.1150750712958503934434410263732282100773_wp real(wp),parameter :: c22 = .025_wp real(wp),parameter :: a2201 = .1802918623893620731908165792176564180038e-1_wp real(wp),parameter :: a2208 = .6983601042028873702545973390560096201728e-1_wp real(wp),parameter :: a2209 = -.2541247660791663512384395986842781657182e-1_wp real(wp),parameter :: a2210 = .8487827035463274491721441398893680307535e-2_wp real(wp),parameter :: a2211 = -.2427525516089801645451101966852425715128e-2_wp real(wp),parameter :: a2212 = -.1047839752893819879012607694745789515746_wp real(wp),parameter :: a2213 = -.1473147795248041942353840372690095884761e-1_wp real(wp),parameter :: a2214 = -.3916338390816177165706892282751065537530e-1_wp real(wp),parameter :: a2215 = -.1005657343293959419073236542225421561652e-1_wp real(wp),parameter :: a2217 = .1102510392204834322538452331445716455061e-1_wp real(wp),parameter :: a2218 = .5092830749095398308703438556315975226108e-2_wp real(wp),parameter :: a2219 = .4759715599420644505591133410826632557391e-1_wp real(wp),parameter :: a2220 = .3386307003288382751110965442296681690349e-1_wp real(wp),parameter :: a2221 = .2764422831404797700452373965825845732168e-1_wp real(wp),parameter :: c23 = .15_wp real(wp),parameter :: a2301 = .1677431640522778042988664067637191163626e-1_wp real(wp),parameter :: a2308 = .6220437408820475326702539861577894278533_wp real(wp),parameter :: a2309 = -.2060859809768841878234097076241307428139_wp real(wp),parameter :: a2310 = .1156394989766058889629372195583391792474_wp real(wp),parameter :: a2311 = .2664101793378358946544219293685167025971e-1_wp real(wp),parameter :: a2312 = -.9376810793418770527505892794460093668860_wp real(wp),parameter :: a2313 = -.1367806466702160302637074581619101741312_wp real(wp),parameter :: a2314 = -.3678480995268296672182605288991379118419_wp real(wp),parameter :: a2315 = -.9547871314402478902820445838193201497337e-1_wp real(wp),parameter :: a2317 = .1013492018422369748729008873270013785313_wp real(wp),parameter :: a2318 = -.8911323084568593396468400926074881389560e-1_wp real(wp),parameter :: a2319 = .4664140988974760478895528270623735057521_wp real(wp),parameter :: a2320 = .4502736292354579812232681662308722738519_wp real(wp),parameter :: a2321 = .1838522463326818655346135218242696774099_wp real(wp),parameter :: c24 = .32_wp real(wp),parameter :: a2401 = .1071149731491444187554380927165768658192e-1_wp real(wp),parameter :: a2408 = -.7094336118221108191937165464264324417735e-1_wp real(wp),parameter :: a2409 = .1002164900340091596740582334112699697590_wp real(wp),parameter :: a2410 = .1383453980468025108839271214703390659581_wp real(wp),parameter :: a2411 = .1796330633578163411338104055485109917477_wp real(wp),parameter :: a2412 = .9048246545576180974879274948815422276563e-1_wp real(wp),parameter :: a2413 = -.5460662294523338383345981122023862069115e-2_wp real(wp),parameter :: a2414 = -.3000457905119619782973021046143166498567e-1_wp real(wp),parameter :: a2415 = -.1145192026962799093665613252151017277867e-1_wp real(wp),parameter :: a2417 = .1003394686109385076849515422360600302176e-1_wp real(wp),parameter :: a2418 = -.9506485282809046129031027932806241113157e-1_wp real(wp),parameter :: a2419 = .4853358804093591445756711642658478691640e-1_wp real(wp),parameter :: a2420 = .8013325919783924638483373011297347396327e-1_wp real(wp),parameter :: a2421 = -.1251643326835242045676140618774248455713_wp real(wp),parameter :: c25 = .78_wp real(wp),parameter :: a2501 = .1410172088869221367153586187761994182069e-1_wp real(wp),parameter :: a2508 = -.3713379753704491105936205420001801316029_wp real(wp),parameter :: a2509 = .2231265548117180273161442520179150684520_wp real(wp),parameter :: a2510 = .1287005345918120122888629169443916280865_wp real(wp),parameter :: a2511 = .2224600659675494761192249831098918110654_wp real(wp),parameter :: a2512 = .5382853042550701952740528638168708946100_wp real(wp),parameter :: a2513 = .5417202616988763101781128062036252796548e-1_wp real(wp),parameter :: a2514 = .1256968791308743925752109039299467082975_wp real(wp),parameter :: a2515 = .2784492789002054061504430663197543089132e-1_wp real(wp),parameter :: a2517 = -.3077409246205059733390460511525401688205e-1_wp real(wp),parameter :: a2518 = .8569805293689777608077303071761466118035e-2_wp real(wp),parameter :: a2519 = -.1535174690587044615794997685221990516897_wp real(wp),parameter :: a2520 = -.2179957030548196497189489878038029238243e-1_wp real(wp),parameter :: a2521 = .1447128819737186799295514239727801525027e-1_wp real(wp),parameter :: c26 = .96_wp real(wp),parameter :: a2601 = .1424600411735646609296566581447532773183e-1_wp real(wp),parameter :: a2608 = -.3767107393295407091303982522049390741260_wp real(wp),parameter :: a2609 = .2252399780730421480874737297000189000070_wp real(wp),parameter :: a2610 = .1283603076292529988314451246143633426068_wp real(wp),parameter :: a2611 = .2230238705261692544876826347415151339678_wp real(wp),parameter :: a2612 = .5463127827750747224899202176094949607118_wp real(wp),parameter :: a2613 = .5526190791375779994553849469706124289752e-1_wp real(wp),parameter :: a2614 = .1285613508749982456581494397108686240388_wp real(wp),parameter :: a2615 = .2857250681296406482698934635829147899039e-1_wp real(wp),parameter :: a2617 = -.2398761886357108720930416967644499057175e-1_wp real(wp),parameter :: a2618 = .5556224458910509454379297181908734648749e-1_wp real(wp),parameter :: a2619 = -.1740675650762838674257930398070760254668e-1_wp real(wp),parameter :: a2620 = -.3815462365996979065575121886854199471011e-1_wp real(wp),parameter :: a2621 = .1111878504898917877407531966545730451506e-1_wp associate (t => me%dense_t, h => me%dense_h, x => me%dense_x, k => me%funcs) ! f at the end of the step (this is also the first stage of the next step, ! so it is saved in the function evaluation cache): call me%set_fsal_cache(t+h, me%dense_xf, k(:,17)) if (me%stopped) return call me%f(t+c18*h, x+h*(a1801*k(:,1) + a1808*k(:,8) + a1809*k(:,9) + a1810*k(:,10) + & a1811*k(:,11) + a1812*k(:,12) + a1813*k(:,13) + a1814*k(:,14) + & a1815*k(:,15) + a1817*k(:,17)), k(:,18)) if (me%stopped) return call me%f(t+c19*h, x+h*(a1901*k(:,1) + a1908*k(:,8) + a1909*k(:,9) + a1910*k(:,10) + & a1911*k(:,11) + a1912*k(:,12) + a1913*k(:,13) + a1914*k(:,14) + & a1915*k(:,15) + a1917*k(:,17) + a1918*k(:,18)), k(:,19)) if (me%stopped) return call me%f(t+c20*h, x+h*(a2001*k(:,1) + a2008*k(:,8) + a2009*k(:,9) + a2010*k(:,10) + & a2011*k(:,11) + a2012*k(:,12) + a2013*k(:,13) + a2014*k(:,14) + & a2015*k(:,15) + a2017*k(:,17) + a2018*k(:,18) + a2019*k(:,19)), k(:,20)) if (me%stopped) return call me%f(t+c21*h, x+h*(a2101*k(:,1) + a2108*k(:,8) + a2109*k(:,9) + a2110*k(:,10) + & a2111*k(:,11) + a2112*k(:,12) + a2113*k(:,13) + a2114*k(:,14) + & a2115*k(:,15) + a2117*k(:,17) + a2118*k(:,18) + a2119*k(:,19) + & a2120*k(:,20)), k(:,21)) if (me%stopped) return call me%f(t+c22*h, x+h*(a2201*k(:,1) + a2208*k(:,8) + a2209*k(:,9) + a2210*k(:,10) + & a2211*k(:,11) + a2212*k(:,12) + a2213*k(:,13) + a2214*k(:,14) + & a2215*k(:,15) + a2217*k(:,17) + a2218*k(:,18) + a2219*k(:,19) + & a2220*k(:,20) + a2221*k(:,21)), k(:,22)) if (me%stopped) return call me%f(t+c23*h, x+h*(a2301*k(:,1) + a2308*k(:,8) + a2309*k(:,9) + a2310*k(:,10) + & a2311*k(:,11) + a2312*k(:,12) + a2313*k(:,13) + a2314*k(:,14) + & a2315*k(:,15) + a2317*k(:,17) + a2318*k(:,18) + a2319*k(:,19) + & a2320*k(:,20) + a2321*k(:,21)), k(:,23)) if (me%stopped) return call me%f(t+c24*h, x+h*(a2401*k(:,1) + a2408*k(:,8) + a2409*k(:,9) + a2410*k(:,10) + & a2411*k(:,11) + a2412*k(:,12) + a2413*k(:,13) + a2414*k(:,14) + & a2415*k(:,15) + a2417*k(:,17) + a2418*k(:,18) + a2419*k(:,19) + & a2420*k(:,20) + a2421*k(:,21)), k(:,24)) if (me%stopped) return call me%f(t+c25*h, x+h*(a2501*k(:,1) + a2508*k(:,8) + a2509*k(:,9) + a2510*k(:,10) + & a2511*k(:,11) + a2512*k(:,12) + a2513*k(:,13) + a2514*k(:,14) + & a2515*k(:,15) + a2517*k(:,17) + a2518*k(:,18) + a2519*k(:,19) + & a2520*k(:,20) + a2521*k(:,21)), k(:,25)) if (me%stopped) return call me%f(t+c26*h, x+h*(a2601*k(:,1) + a2608*k(:,8) + a2609*k(:,9) + a2610*k(:,10) + & a2611*k(:,11) + a2612*k(:,12) + a2613*k(:,13) + a2614*k(:,14) + & a2615*k(:,15) + a2617*k(:,17) + a2618*k(:,18) + a2619*k(:,19) + & a2620*k(:,20) + a2621*k(:,21)), k(:,26)) end associate end procedure rkv98e_dense_extra_stages !***************************************************************************************** !***************************************************************************************** !> ! Dense output for Verner's "most robust" 9(8) method [[rkv98r]]. ! This is the interpolant of order 9 (`bi9`) from the coefficient file, which requires 10 additional ! function evaluations (see [[rkv98r_dense_extra_stages]]). ! !@note The interpolant is published in the polynomial form. Here it has been ! converted (exactly, using high precision arithmetic) ! to the form used in `DOPRI5` and `DOP853` (see [[dense_coefficients]]), which has ! much smaller roundoff error for these coefficients. This code is generated by ! `scripts/generate_dense.py`. ! !### References ! * J.H. Verner, "Numerically optimal Runge--Kutta pairs with interpolants", ! Numerical Algorithms, 53, (2010) pp. 383-396. [doi](https://doi.org/10.1007/s11075-009-9290-3) ! * Coefficients: ! * [RKV98.IIa.Robust.000000351.081209.CoeffsOnlyFLOAT6040](https://web.archive.org/web/https://www.sfu.ca/~jverner/RKV98.IIa.Robust.000000351.081209.CoeffsOnlyFLOAT6040) ! * [RKV98.IIa.Robust.000000351.081209.CoeffsOnlyRADandFLOATS40](https://web.archive.org/web/https://www.sfu.ca/~jverner/RKV98.IIa.Robust.000000351.081209.CoeffsOnlyRADandFLOATS40) module procedure rkv98r_dense_coefficients form = dense_form_dopri allocate(b(26,9)) b = 0.0_wp b(1,:) = [ 1.458885278405539719101539582255752906000e-2_wp, & 9.854111472159446028089846041774424709400e-1_wp, & -9.708222944318892056179692083548849402000e-1_wp, & -2.637066763385014715848211085903324897380e1_wp, & 2.660320174559264463682575087921035465450e1_wp, & 1.772067251706265541061912112509357845920e2_wp, & -1.799086577126680615986898704592441623500e2_wp, & -3.598290182379035162462470562409222601692e2_wp, & 3.671554574315021425912203561515071665187e2_wp ] b(8,:) = [ 2.024197887889332650566666683195656101220e-3_wp, & -2.024197887889332650566666683195656101220e-3_wp, & 4.048395775778665301133333366391312167700e-3_wp, & -1.960136618337276162769723133916471559300e-2_wp, & 5.531857983304611009285027556921979028500e-2_wp, & 4.011844228456285693801130630206043463000e-1_wp, & -1.024276754469481646618879849564428451780e0_wp, & -2.078396401843996195208891754419523644080e0_wp, & 4.674158442717667575904712851691470968655e0_wp ] b(9,:) = [ 2.178047084569716646796256135839225743450e-1_wp, & -2.178047084569716646796256135839225743450e-1_wp, & 4.356094169139433293592512271678451494000e-1_wp, & -2.109116837079348887399766793715451622800e0_wp, & 5.952306948286395162126274481842078105000e0_wp, & 4.316764520808936940164061738905736695400e1_wp, & -1.102126927516463968944139861508143575190e2_wp, & -2.236364957547074442393208735387683032966e2_wp, & 5.029417938773550941496020410671334197216e2_wp ] b(10,:) = [ 1.274895340854389692868677968654808675533e-1_wp, & -1.274895340854389692868677968654808675533e-1_wp, & 2.549790681708779385737355937309617299200e-1_wp, & -1.234547796491375740392284556721378495900e0_wp, & 3.484115862079580772279676232951059583000e0_wp, & 2.526769514825283360713661763226758842900e1_wp, & -6.451175894567435092846231016374016783300e1_wp, & -1.309031051268592972887356736771996492388e2_wp, & 2.943913170095037265437261593264399863548e2_wp ] b(11,:) = [ 2.244617745463131861258531547137348033760e-1_wp, & -2.244617745463131861258531547137348033760e-1_wp, & 4.489235490926263722517063094274696059100e-1_wp, & -2.173580687627161778645801801047976361300e0_wp, & 6.134235525586262681747094992931877088000e0_wp, & 4.448703756240235718898809565260904499400e1_wp, & -1.135812754821567698287476571360806256760e2_wp, & -2.304718068128333713234406086507888172788e2_wp, & 5.183138984780918128529461223609697359788e2_wp ] b(12,:) = [ 1.787254491259903095100090833796054441670e-1_wp, & -1.787254491259903095100090833796054441670e-1_wp, & 3.574508982519806190200181667592108911600e-1_wp, & -1.730691942505290489911333018731866402000e0_wp, & 4.884323852339505668617782888843418892000e0_wp, & 3.542235993053079998015926093615899385900e1_wp, & -9.043795770518971042229120948659313689200e1_wp, & -1.835108773721431312385479901153920605670e2_wp, & 4.127022718276975364384031516482241464250e2_wp ] b(13,:) = [ 7.594344758096557172908303416513173071020e-2_wp, & -7.594344758096557172908303416513173071020e-2_wp, & 1.518868951619311434581660683302634615300e-1_wp, & -7.354000980677168332946589478840249153000e-1_wp, & 2.075431306859499574776849039395681721400e0_wp, & 1.505155615908963731350887337540990917380e1_wp, & -3.842860842650351340235082956830033617980e1_wp, & -7.797685648239174193825861324061456007220e1_wp, & 1.753641325304376264979817562401272782969e2_wp ] b(14,:) = [ 1.294845879197561516869001434704642294505e-1_wp, & -1.294845879197561516869001434704642294505e-1_wp, & 2.589691758395123033738002869409284523100e-1_wp, & -1.253866945570603997411573624060682600200e0_wp, & 3.538637974500109895677536309218684514000e0_wp, & 2.566310338667416323519405481173356503000e1_wp, & -6.552128833933312965777341363663141797700e1_wp, & -1.329515771342347930809400620793492718190e2_wp, & 2.989981777216614016801117759308092190290e2_wp ] b(15,:) = [ 2.947744761261941714007911131590716599580e-2_wp, & -2.947744761261941714007911131590716599580e-2_wp, & 5.895489522523883428015822263181433231500e-2_wp, & -2.854455329012426521536430973146036065400e-1_wp, & 8.055786189627093137627906185284474536000e-1_wp, & 5.842261212791819782825976489021398248500e0_wp, & -1.491606356836032223347328803975834579720e1_wp, & -3.026671523577966975300765540210259557502e1_wp, & 6.806758442573064372308655760592173752592e1_wp ] b(17,:) = [ 1.131000000000000000000000000000000000000e-37_wp, & -1.131000000000000000000000000000000000000e-37_wp, & -1.000000000000000000000000000000000000970e0_wp, & 5.458150170288863269034998372331305350000e-1_wp, & 1.817661909544495089602219720966551185990e1_wp, & -1.383717924690821375893089544366403067000e0_wp, & -1.147734767597614978790533863269211013350e2_wp, & -1.034198139378041982072302882124553707820e1_wp, & 2.518328522135679646449465513799755444867e2_wp ] b(18,:) = [ -7.070000000000000000000000000000000000000e-38_wp, & 7.070000000000000000000000000000000000000e-38_wp, & 4.220000000000000000000000000000000000000e-37_wp, & -1.642163800536330788290065204813354325100e-1_wp, & 9.519929365114165055101380496084025410500e-1_wp, & 4.739153673458544025901423505009778408500e0_wp, & -2.131184933122730491680094714445253547520e1_wp, & -1.831289306990792069265916156113594432898e1_wp, & 7.591229137507532480130986827669750176353e1_wp ] b(19,:) = [ 3.900000000000000000000000000000000000000e-37_wp, & -3.900000000000000000000000000000000000000e-37_wp, & -1.790000000000000000000000000000000000000e-36_wp, & -6.748844193125564745929227754653454946000e-1_wp, & 3.912430659702321225670397405943369852100e0_wp, & 1.947662574160048447368276443195596868200e1_wp, & -8.758587331960754878214405437565889419200e1_wp, & -7.526098310887908074617101188121033899940e1_wp, & 3.119787603807909928458652100409945704424e2_wp ] b(20,:) = [ 7.730000000000000000000000000000000000000e-37_wp, & -7.730000000000000000000000000000000000000e-37_wp, & -2.740000000000000000000000000000000000000e-36_wp, & 1.778873344420799933681245486496450244300e0_wp, & -1.031246004986806282727141824566679162800e1_wp, & -5.133686506836873201794597274110383607500e1_wp, & 2.308605309850987048297641643879070668530e2_wp, & 1.983743481049099109017199049274608500936e2_wp, & -8.223196224208751607207370309035750144516e2_wp ] b(21,:) = [ 4.316000000000000000000000000000000000000e-38_wp, & -4.316000000000000000000000000000000000000e-38_wp, & -1.980000000000000000000000000000000000000e-37_wp, & 4.134694024574781779455481610866606195000e-2_wp, & -2.396959124750886878829081649079368525300e-1_wp, & -1.193239698061238269798491979757585653620e0_wp, & 5.365967515158004398451414547972962168220e0_wp, & 4.610880444697151389389213911963751684110e0_wp, & -1.911344638322878584148388848904209912604e1_wp ] b(22,:) = [ 1.040000000000000000000000000000000000000e-36_wp, & -1.040000000000000000000000000000000000000e-36_wp, & -5.500000000000000000000000000000000000000e-36_wp, & 3.131247017428404500084714290987346868600e1_wp, & -3.387022087901636525344456435703796891050e1_wp, & -2.588205663021983641434655081348870539980e2_wp, & 2.854737764182181819679920389377461350840e2_wp, & 5.885263279733103419273222746605193404750e2_wp, & -6.742890614799327434029117465686431697155e2_wp ] b(23,:) = [ -2.550000000000000000000000000000000000000e-37_wp, & 2.550000000000000000000000000000000000000e-37_wp, & 6.720000000000000000000000000000000000000e-37_wp, & -2.730958729620607116715051590434636563256e0_wp, & 1.566610140744408727858992403426037489610e0_wp, & 8.002802971584994113128834556971356162016e1_wp, & -5.246852711792871340459703830423998954264e1_wp, & -1.596405935562192937597599467349257817650e2_wp, & -7.090434387268304770405112504542156157026e0_wp ] b(24,:) = [ 4.250000000000000000000000000000000000000e-37_wp, & -4.250000000000000000000000000000000000000e-37_wp, & -1.300000000000000000000000000000000000000e-36_wp, & 3.159813308843413258842148106168074298100e0_wp, & -7.419639103912373270945802130100008608000e0_wp, & -9.943272818292694588113489712183004393500e1_wp, & 1.916576867907068770537324986336022165370e2_wp, & 4.642139857286974302377727707301602292418e2_wp, & -8.133936179360131427508167163638763606878e2_wp ] b(25,:) = [ 8.500000000000000000000000000000000000000e-38_wp, & -8.500000000000000000000000000000000000000e-38_wp, & -7.000000000000000000000000000000000000000e-38_wp, & 9.128020953320920409453895709005633217000e-1_wp, & 2.747583973691318169360530714969178335000e0_wp, & -3.374203780936702542849738237632098313140e1_wp, & -3.229324367913034554974039429876197988960e1_wp, & 2.208650775021236000229375410529858952982e2_wp, & -1.575322311379308302075254133968369612868e2_wp ] b(26,:) = [ -1.078000000000000000000000000000000000000e-36_wp, & 1.078000000000000000000000000000000000000e-36_wp, & 5.390000000000000000000000000000000000000e-36_wp, & 1.731857489108072590441870089449362046100e0_wp, & -2.904637127486227930078416860439061586550e1_wp, & -3.084422234659900569916201220862765848500e1_wp, & 2.736175881844753788952171484335330984650e2_wp, & 1.585906799337452418438788684149845770398e2_wp, & -7.885942819688829666512243546539760160873e2_wp ] end procedure rkv98r_dense_coefficients !***************************************************************************************** !***************************************************************************************** !> ! The additional function evaluations for the dense output of [[rkv98r]]. ! These are stored in `me%funcs(:,17:26)`. Stage 17 is `f` at the end of the step, ! which is cached for the first stage of the next step (see [[set_fsal_cache]]). module procedure rkv98r_dense_extra_stages real(wp),parameter :: c18 = .7421010083583088006243877166286162723673_wp real(wp),parameter :: a1801 = .1560140526108861651166696161636643678715e-1_wp real(wp),parameter :: a1808 = .2681164393327584462067257298507705692190_wp real(wp),parameter :: a1809 = .1883053124587791093781990103116496723062_wp real(wp),parameter :: a1810 = .1249199137461030776325034278436425745187_wp real(wp),parameter :: a1811 = .2302302127814521984877267202569592992854_wp real(wp),parameter :: a1812 = -.1360312216132798635363504089171114829643_wp real(wp),parameter :: a1813 = .7488659971306952790649243500419796404442e-1_wp real(wp),parameter :: a1814 = -.2812840029795629049644800354887547182869e-1_wp real(wp),parameter :: a1815 = -.2314455726481949553795911074427483913356e-1_wp real(wp),parameter :: a1817 = .2734530424111347407183095495529155013292e-1_wp real(wp),parameter :: c19 = .888_wp real(wp),parameter :: a1901 = .1311195721844068401556584238625026007821e-1_wp real(wp),parameter :: a1908 = -.1464024265969826989710171144006401375029_wp real(wp),parameter :: a1909 = .2471264389666795959943268183777606168285_wp real(wp),parameter :: a1910 = .1311375203080032288323212687159488912255_wp real(wp),parameter :: a1911 = .2170560346982582529677346605154232633905_wp real(wp),parameter :: a1912 = .2867536713760320279685175119456430108611_wp real(wp),parameter :: a1913 = .2323311339149421811667180525022488323943e-1_wp real(wp),parameter :: a1914 = .5250677264199395668210993035986718117205e-1_wp real(wp),parameter :: a1915 = .2833951586009950740445688183974089075836e-2_wp real(wp),parameter :: a1917 = -.8502403851995712467765683685024083445736e-2_wp real(wp),parameter :: a1918 = .6914537026206649612108927235057202507752e-1_wp real(wp),parameter :: c20 = .696_wp real(wp),parameter :: a2001 = .1398921213361768395172663025684236914687e-1_wp real(wp),parameter :: a2008 = -.3157406517950500155656872013463001170279e-1_wp real(wp),parameter :: a2009 = .2271812513272158193884043729886141068401_wp real(wp),parameter :: a2010 = .1289486410996786560165007139967355200862_wp real(wp),parameter :: a2011 = .2216682589135277018738330104776507206003_wp real(wp),parameter :: a2012 = .1948368236542480784645408517950426740505_wp real(wp),parameter :: a2013 = .5740088404417652948912844658374513373184e-1_wp real(wp),parameter :: a2014 = .9008366542675954700472151691971957675586e-1_wp real(wp),parameter :: a2015 = .1579153208844212165580559397387558558087e-1_wp real(wp),parameter :: a2017 = -.1899131505909185814243914464042645862057e-1_wp real(wp),parameter :: a2018 = -.8830926811918835070775190595441414629725e-1_wp real(wp),parameter :: a2019 = -.1150256203298809274379013662627550701719_wp real(wp),parameter :: c21 = .487_wp real(wp),parameter :: a2101 = .1615147291900762304051270768727488699116e-1_wp real(wp),parameter :: a2108 = .8098685003242905691737114344809608208533e-1_wp real(wp),parameter :: a2109 = .1276916294306930449548063682447180310200_wp real(wp),parameter :: a2110 = .1234814359383480575498127887662911829596_wp real(wp),parameter :: a2111 = .2339851259140109824702007335541675692697_wp real(wp),parameter :: a2112 = -.6595995683357367659329999201079842018819e-1_wp real(wp),parameter :: a2113 = -.2565276859406432791053055820612070574951e-1_wp real(wp),parameter :: a2114 = -.1258973463819247076247947215937788425689_wp real(wp),parameter :: a2115 = -.4307672490364843615025132534342781501874e-1_wp real(wp),parameter :: a2117 = .4973042479196705159702504356466546825472e-1_wp real(wp),parameter :: a2118 = .1000473540179392626782905455570668634250_wp real(wp),parameter :: a2119 = .1378658806763623213069940499298119480365_wp real(wp),parameter :: a2120 = -.1223533770075462522361367835979662485167_wp real(wp),parameter :: c22 = .25e-1_wp real(wp),parameter :: a2201 = .1802729411106818689763249275362776542125e-1_wp real(wp),parameter :: a2208 = -.3566948654464566871345751863957855326182e-3_wp real(wp),parameter :: a2209 = -.2359965598171250974225291789331305531887e-1_wp real(wp),parameter :: a2210 = .8489537818582623161952771938798430467577e-2_wp real(wp),parameter :: a2211 = -.2424017121415795431658865084807451176669e-2_wp real(wp),parameter :: a2212 = -.3430595122831044180368089698215097859892e-1_wp real(wp),parameter :: a2213 = -.1890658549297046191788725638037661418072e-1_wp real(wp),parameter :: a2214 = -.3742593046527093615091178254251506793186e-1_wp real(wp),parameter :: a2215 = -.9695346568972416886241067894611986071147e-2_wp real(wp),parameter :: a2217 = .1102424057658581840569553621781637903330e-1_wp real(wp),parameter :: a2218 = .4923260682022719983443039400621144750980e-2_wp real(wp),parameter :: a2219 = .4757764286275925144559532694675125937637e-1_wp real(wp),parameter :: a2220 = .3403995941160684337366745033612222756612e-1_wp real(wp),parameter :: a2221 = .2763224626147357535178074437043373219519e-1_wp real(wp),parameter :: c23 = .15_wp real(wp),parameter :: a2301 = .1674115470096526424930400903558890871308e-1_wp real(wp),parameter :: a2308 = -.3173492444543278603351394583591742321103e-2_wp real(wp),parameter :: a2309 = -.1907932208630995154539013764218262860213_wp real(wp),parameter :: a2310 = .1156811286400889022914973296917837873123_wp real(wp),parameter :: a2311 = .2659243297776452866445546052304774936345e-1_wp real(wp),parameter :: a2312 = -.3088649490220827144035155968652816431108_wp real(wp),parameter :: a2313 = -.1753758456547022422433014122559418385679_wp real(wp),parameter :: a2314 = -.3519244231945119122948228740531438774071_wp real(wp),parameter :: a2315 = -.9209686132127586046007074704921651582851e-1_wp real(wp),parameter :: a2317 = .1013643084710043289094345867037510169960_wp real(wp),parameter :: a2318 = -.8614613115672610072132806092797279629388e-1_wp real(wp),parameter :: a2319 = .4667555353741319512606533809977766283817_wp real(wp),parameter :: a2320 = .4471784578431838563500095104120691993359_wp real(wp),parameter :: a2321 = .1840619056498027924549371847929574094481_wp real(wp),parameter :: c24 = .32_wp real(wp),parameter :: a2401 = .1067816157945957352545935273312729770712e-1_wp real(wp),parameter :: a2408 = .3813266550650219306405894185406227733203e-3_wp real(wp),parameter :: a2409 = .9563838511927889653937246525395111512962e-1_wp real(wp),parameter :: a2410 = .1383796710757440233248935289235401276023_wp real(wp),parameter :: a2411 = .1794589653248393403412396094337145578219_wp real(wp),parameter :: a2412 = .2328104103984205159658342243398415303992e-1_wp real(wp),parameter :: a2413 = -.6102337776291304739559151883232786325123e-2_wp real(wp),parameter :: a2414 = -.2957864564982467043648509518504953052119e-1_wp real(wp),parameter :: a2415 = -.1107556819372404441098997386306135386770e-1_wp real(wp),parameter :: a2417 = .1005006241909649985189552106106942758409e-1_wp real(wp),parameter :: a2418 = -.9189958889836397344227776088639708677529e-1_wp real(wp),parameter :: a2419 = .4889782809929091409586999965373170403839e-1_wp real(wp),parameter :: a2420 = .7683137001604329651351035116374049373690e-1_wp real(wp),parameter :: a2421 = -.1249406708104556246901528582576587419443_wp real(wp),parameter :: c25 = .78_wp real(wp),parameter :: a2501 = .1407873162514538134880663877417425640142e-1_wp real(wp),parameter :: a2508 = .1921038236500748458986816137745896655262e-2_wp real(wp),parameter :: a2509 = .2103911064595803591331678030844204745297_wp real(wp),parameter :: a2510 = .1287126282116406147526805293017325463897_wp real(wp),parameter :: a2511 = .2222827890207636907592258288998389315181_wp real(wp),parameter :: a2512 = .1689157537655836922893273799046879855404_wp real(wp),parameter :: a2513 = .7069536651338852930255328632439555742574e-1_wp real(wp),parameter :: a2514 = .1192421401475245249188274865443088811860_wp real(wp),parameter :: a2515 = .2685261833359796204499227534872501897484e-1_wp real(wp),parameter :: a2517 = -.3077554523024171431258004671769555940844e-1_wp real(wp),parameter :: a2518 = .8284466445798771733365157887905394914006e-2_wp real(wp),parameter :: a2519 = -.1535503041842782550897093173814614623206_wp real(wp),parameter :: a2520 = -.2150191510960740019354801525582536987010e-1_wp real(wp),parameter :: a2521 = .1445112576460309485390417714704744806405e-1_wp real(wp),parameter :: c26 = .96_wp real(wp),parameter :: a2601 = .1422598489194661552790199000685205084824e-1_wp real(wp),parameter :: a2608 = .1949172492119241778367670693682637328188e-2_wp real(wp),parameter :: a2609 = .2124063537210806078797518674148478107397_wp real(wp),parameter :: a2610 = .1283651739828476111356431808186160040103_wp real(wp),parameter :: a2611 = .2228603986302258483512559751598074395656_wp real(wp),parameter :: a2612 = .1715923679711209630782228726709623464662_wp real(wp),parameter :: a2613 = .7212912939761323676316514664066037355886e-1_wp real(wp),parameter :: a2614 = .1220420686370054673857430981768668520159_wp real(wp),parameter :: a2615 = .2757049488622773377977068820958406578502e-1_wp real(wp),parameter :: a2617 = -.2399703787121476720950542355073551377435e-1_wp real(wp),parameter :: a2618 = .5371225309992069259478306291407462056681e-1_wp real(wp),parameter :: a2619 = -.1761964270699028520222472394692440648034e-1_wp real(wp),parameter :: a2620 = -.3622477925515253053764132865960991436725e-1_wp real(wp),parameter :: a2621 = .1098806212324956467476592345131563373716e-1_wp associate (t => me%dense_t, h => me%dense_h, x => me%dense_x, k => me%funcs) ! f at the end of the step (this is also the first stage of the next step, ! so it is saved in the function evaluation cache): call me%set_fsal_cache(t+h, me%dense_xf, k(:,17)) if (me%stopped) return call me%f(t+c18*h, x+h*(a1801*k(:,1) + a1808*k(:,8) + a1809*k(:,9) + a1810*k(:,10) + & a1811*k(:,11) + a1812*k(:,12) + a1813*k(:,13) + a1814*k(:,14) + & a1815*k(:,15) + a1817*k(:,17)), k(:,18)) if (me%stopped) return call me%f(t+c19*h, x+h*(a1901*k(:,1) + a1908*k(:,8) + a1909*k(:,9) + a1910*k(:,10) + & a1911*k(:,11) + a1912*k(:,12) + a1913*k(:,13) + a1914*k(:,14) + & a1915*k(:,15) + a1917*k(:,17) + a1918*k(:,18)), k(:,19)) if (me%stopped) return call me%f(t+c20*h, x+h*(a2001*k(:,1) + a2008*k(:,8) + a2009*k(:,9) + a2010*k(:,10) + & a2011*k(:,11) + a2012*k(:,12) + a2013*k(:,13) + a2014*k(:,14) + & a2015*k(:,15) + a2017*k(:,17) + a2018*k(:,18) + a2019*k(:,19)), k(:,20)) if (me%stopped) return call me%f(t+c21*h, x+h*(a2101*k(:,1) + a2108*k(:,8) + a2109*k(:,9) + a2110*k(:,10) + & a2111*k(:,11) + a2112*k(:,12) + a2113*k(:,13) + a2114*k(:,14) + & a2115*k(:,15) + a2117*k(:,17) + a2118*k(:,18) + a2119*k(:,19) + & a2120*k(:,20)), k(:,21)) if (me%stopped) return call me%f(t+c22*h, x+h*(a2201*k(:,1) + a2208*k(:,8) + a2209*k(:,9) + a2210*k(:,10) + & a2211*k(:,11) + a2212*k(:,12) + a2213*k(:,13) + a2214*k(:,14) + & a2215*k(:,15) + a2217*k(:,17) + a2218*k(:,18) + a2219*k(:,19) + & a2220*k(:,20) + a2221*k(:,21)), k(:,22)) if (me%stopped) return call me%f(t+c23*h, x+h*(a2301*k(:,1) + a2308*k(:,8) + a2309*k(:,9) + a2310*k(:,10) + & a2311*k(:,11) + a2312*k(:,12) + a2313*k(:,13) + a2314*k(:,14) + & a2315*k(:,15) + a2317*k(:,17) + a2318*k(:,18) + a2319*k(:,19) + & a2320*k(:,20) + a2321*k(:,21)), k(:,23)) if (me%stopped) return call me%f(t+c24*h, x+h*(a2401*k(:,1) + a2408*k(:,8) + a2409*k(:,9) + a2410*k(:,10) + & a2411*k(:,11) + a2412*k(:,12) + a2413*k(:,13) + a2414*k(:,14) + & a2415*k(:,15) + a2417*k(:,17) + a2418*k(:,18) + a2419*k(:,19) + & a2420*k(:,20) + a2421*k(:,21)), k(:,24)) if (me%stopped) return call me%f(t+c25*h, x+h*(a2501*k(:,1) + a2508*k(:,8) + a2509*k(:,9) + a2510*k(:,10) + & a2511*k(:,11) + a2512*k(:,12) + a2513*k(:,13) + a2514*k(:,14) + & a2515*k(:,15) + a2517*k(:,17) + a2518*k(:,18) + a2519*k(:,19) + & a2520*k(:,20) + a2521*k(:,21)), k(:,25)) if (me%stopped) return call me%f(t+c26*h, x+h*(a2601*k(:,1) + a2608*k(:,8) + a2609*k(:,9) + a2610*k(:,10) + & a2611*k(:,11) + a2612*k(:,12) + a2613*k(:,13) + a2614*k(:,14) + & a2615*k(:,15) + a2617*k(:,17) + a2618*k(:,18) + a2619*k(:,19) + & a2620*k(:,20) + a2621*k(:,21)), k(:,26)) end associate end procedure rkv98r_dense_extra_stages !***************************************************************************************** !***************************************************************************************** !> ! Dense output for the Prince-Dormand 8(7) (RK8(7)13M) method [[rkdp87]]. ! This is an interpolant of order 8 constructed by P. Stone, which requires 8 additional ! function evaluations (see [[rkdp87_dense_extra_stages]]). ! !@note The interpolant is published in the polynomial form. Here it has been ! converted (exactly, using high precision arithmetic) ! to the form used in `DOPRI5` and `DOP853` (see [[dense_coefficients]]), which has ! much smaller roundoff error for these coefficients. This code is generated by ! `scripts/generate_dense.py`. ! !### References ! * P. Stone, [Maple worksheets on interpolation for Runge-Kutta schemes](http://www.peterstone.name/Maplepgs/RKinterp.html) ! * Coefficients: ! * [RKinterp8c.mws](http://www.peterstone.name/Maplepgs/Maple/nmthds/RKinterp/RKinterp8c.mws) module procedure rkdp87_dense_coefficients form = dense_form_dopri allocate(b(21,8)) b = 0.0_wp b(1,:) = [ 4.174749114153024622208592846850711513419e-2_wp, & 9.582525088584697537779140715314928848658e-1_wp, & -9.165050177169395075558281430629857697316e-1_wp, & -1.067510547241701717731728121412964567601e1_wp, & 5.691786028593493178164148738407631635299e0_wp, & 4.391912066338726224652799717492210092833e1_wp, & -1.017913579379592465489590917619992278487e1_wp, & -6.312660810334737675623449298204836161855e1_wp ] b(6,:) = [ -5.545232861123930896152189465471671889359e-2_wp, & 5.545232861123930896152189465471671889359e-2_wp, & -1.109046572224786179230437893094334377872e-1_wp, & 3.449518777036079206627691527553235890416e-1_wp, & 4.338661596866255964339784301132702041197e0_wp, & -6.463813607794334373215951712190746917385e0_wp, & -1.176736972872198138487472967397777592899e1_wp, & 1.618183016733174273109322680926209644941e1_wp ] b(7,:) = [ 2.393128072011800970467473542487569696603e-1_wp, & -2.393128072011800970467473542487569696603e-1_wp, & 4.786256144023601940934947084975139393206e-1_wp, & -1.488691354718631521455725986235172737537e0_wp, & -1.872414220007294806896720176025009142402e1_wp, & 2.789555314351439420032820337417064200403e1_wp, & 5.078384179134131409008911601834228131304e1_wp, & -6.983510521525696384845573954192793856731e1_wp ] b(8,:) = [ 7.035106694034430230580464108897021513664e-1_wp, & -7.035106694034430230580464108897021513664e-1_wp, & 1.407021338806886046116092821779404302733e0_wp, & -4.376323456073095709295351966539791207022e0_wp, & -5.504358068937323180233713560418213287751e1_wp, & 8.200488513293557254437746560142317967795e1_wp, & 1.492898560312776123727914157616806048263e2_wp, & -2.052950788235249679998092356640622140626e2_wp ] b(9,:) = [ -7.597596138144609298844876770850584076554e-1_wp, & 7.597596138144609298844876770850584076554e-1_wp, & -1.519519227628921859768975354170116815311e0_wp, & 4.726230835607266833542837137549102782132e0_wp, & 5.944457053222149283315986829128732491297e1_wp, & -8.856155644708327385857151489902701752935e1_wp, & -1.612262731721192057608006152143643898497e2_wp, & 2.217093735582331955708414259793933941582e2_wp ] b(10,:) = [ 6.605630309222863414613785948378206399404e-1_wp, & -6.605630309222863414613785948378206399404e-1_wp, & 1.321126061844572682922757189675641279881e0_wp, & -4.109159408898925537704689622948532167116e0_wp, & -5.168330215065515345282938146817937363275e1_wp, & 7.699868364438589400402303681372127594555e1_wp, & 1.401760685016980925129405039305637625478e2_wp, & -1.927623068120502484076542206283109961234e2_wp ] b(11,:) = [ 1.581874825101233355296148386006854439729e-1_wp, & -1.581874825101233355296148386006854439729e-1_wp, & 3.163749650202466710592296772013708879457e-1_wp, & -9.840356660876780719349460025702418826601e-1_wp, & -1.237679233063835761691372370078031298946e1_wp, & 1.843916076455660441626278493459138111024e1_wp, & 3.356848377283612665559626372654010864371e1_wp, & -4.616150557936628828786353650482791038667e1_wp ] b(12,:) = [ -2.381095387528628044718635553056971935251e-1_wp, & 2.381095387528628044718635553056971935251e-1_wp, & -4.762190775057256089437271106113943870503e-1_wp, & 1.481206191858501807837682005759356370452e0_wp, & 1.862999692722002483251108221071441697523e1_wp, & -2.775529387641327426196780595241034538870e1_wp, & -5.052849985947183238644455534612397732624e1_wp, & 6.948397324002675250289610412902008999206e1_wp ] b(13,:) = [ 2.500000000000000000000000000000000000000e-1_wp, & -2.500000000000000000000000000000000000000e-1_wp, & 5.000000000000000000000000000000000000000e-1_wp, & -1.555173093460008625853261033667484280365e0_wp, & -1.956032192662003899264263284937968363182e1_wp, & 2.914130826277068960086479182976538490826e1_wp, & 5.305173841850584552705788953878580092085e1_wp, & -7.295379009589482768757034169130799756939e1_wp ] b(14,:) = [ 0.0_wp, & 0.0_wp, & -1.000000000000000000000000000000000000000e0_wp, & 1.248095493190467055679458993877508147821e0_wp, & 8.962759864628559365773153140423005072906e0_wp, & -1.124005230834643269771495028030076614127e1_wp, & -1.951065477078916606492101697633305894155e1_wp, & 2.420108837245054261940413358615466057839e1_wp ] b(15,:) = [ 0.0_wp, & 0.0_wp, & 0.0_wp, & -1.152381207848544894098573138707010290294e2_wp, & -8.161678854828192264319460627035862542498e-3_wp, & 7.043462141070517784109568359798449690394e2_wp, & 1.024840158775703415943290820036273208816e1_wp, & -1.226301165556452677636627649159032906282e3_wp ] b(16,:) = [ 0.0_wp, & 0.0_wp, & 0.0_wp, & 1.072227416468638830508933092707742151603e2_wp, & 7.593993873691578364746608574312709989164e-3_wp, & -6.553554642403781657605036293623812584235e2_wp, & -9.535574758190635465566845864870725099858e0_wp, & 1.141005876876366910230959221355828940894e3_wp ] b(17,:) = [ 0.0_wp, & 0.0_wp, & 0.0_wp, & 1.774521499088860140238969411753059848836e1_wp, & 1.256795450837918899229110230004098454672e-3_wp, & -1.084604201476251091257368789195346690825e2_wp, & -1.578124393639141734281636143575494277425e0_wp, & 1.888348896890075762452210450748404568992e2_wp ] b(18,:) = [ 0.0_wp, & 0.0_wp, & 0.0_wp, & -2.499153732923508241320113753205026447962e2_wp, & -2.085560865731488099393454396013556178237e2_wp, & 3.393477035209722648264943891290367210121e3_wp, & 1.656538945060443836168205178805513388287e3_wp, & -1.076617611968236409826793469998838358377e4_wp ] b(19,:) = [ 0.0_wp, & 0.0_wp, & 0.0_wp, & 1.912933570385397538519982898913037415943e2_wp, & 2.213258557752666454057813519133363975438e2_wp, & -3.216323190099662606833949730704788310879e3_wp, & -1.676047958969323894411391067286864400552e3_wp, & 1.074978832006861187030893711193027141323e4_wp ] b(20,:) = [ 0.0_wp, & 0.0_wp, & 0.0_wp, & 6.098424439680860344443292585461926346588e1_wp, & 1.276599286562775583408554699059797423503e1_wp, & -1.848696209409996464590679748517406237792e2_wp, & -2.160892130174864101287960772426191300084e1_wp, & -8.453163096986787687019747170166512119314e0_wp ] b(21,:) = [ 0.0_wp, & 0.0_wp, & 0.0_wp, & 3.295940057399984817992978593125432178041e0_wp, & 3.478391316961461115422092314005047901693e1_wp, & -7.719254926002200031755657031643240559417e1_wp, & -1.316748224160594386100572925752170208657e2_wp, & 2.398594909932156463698173944652973683008e2_wp ] end procedure rkdp87_dense_coefficients !***************************************************************************************** !***************************************************************************************** !> ! The additional function evaluations for the dense output of [[rkdp87]]. ! These are stored in `me%funcs(:,14:21)`. Stage 14 is `f` at the end of the step, ! which is cached for the first stage of the next step (see [[set_fsal_cache]]). module procedure rkdp87_dense_extra_stages real(wp),parameter :: c15 = 3.110177634953863863927417318829099695921e-1_wp real(wp),parameter :: a1501 = 4.435331117999704827201756397132859707678e-2_wp real(wp),parameter :: a1506 = 4.364686831550274979708965624817361214833e-2_wp real(wp),parameter :: a1507 = 2.214612909096201333992579232656568309507e-1_wp real(wp),parameter :: a1508 = 2.408257402439057850132725800478257151750e-2_wp real(wp),parameter :: a1509 = -3.180313138118299747698448283626399728695e-2_wp real(wp),parameter :: a1510 = 8.864429685322726493007988822309060325983e-3_wp real(wp),parameter :: a1511 = 9.993963109955243576776602485537222732890e-4_wp real(wp),parameter :: a1512 = 1.200805999460228573175848880356799421183e-2_wp real(wp),parameter :: a1514 = -1.259503554386166268241032464519842162533e-2_wp real(wp),parameter :: c16 = 3.900000000000000000000000000000000000000e-1_wp real(wp),parameter :: a1601 = 4.460675166606287059090925020917971268785e-2_wp real(wp),parameter :: a1606 = -1.001717952276942783566249927518240212486e-2_wp real(wp),parameter :: a1607 = 2.180133337532843995237192727488683885895e-1_wp real(wp),parameter :: a1608 = 1.316435515076679550847863952857179920145e-1_wp real(wp),parameter :: a1609 = -1.054077497187771551607663761052289530586e-1_wp real(wp),parameter :: a1610 = 3.481381867357951093264146787230837884826e-2_wp real(wp),parameter :: a1611 = -7.424646830363417175728177729584917674680e-4_wp real(wp),parameter :: a1612 = 8.834418216000000000000000000000000000000e-3_wp real(wp),parameter :: a1614 = -8.750783908552131428208782515955331903073e-3_wp real(wp),parameter :: a1615 = 7.700630401654032001015408955325070671395e-2_wp real(wp),parameter :: c17 = 7.500000000000000000000000000000000000000e-2_wp real(wp),parameter :: a1701 = 3.993020695653906563383468777519318548746e-2_wp real(wp),parameter :: a1706 = 3.107501798295948025417866480313858185119e-1_wp real(wp),parameter :: a1707 = 8.046434929536454903337676393017676713636e-2_wp real(wp),parameter :: a1708 = -6.157526423016600959408280648540208463029e-1_wp real(wp),parameter :: a1709 = 3.903010728140840942288433275705286336872e-1_wp real(wp),parameter :: a1710 = -1.422085411509405491050826041202067568781e-1_wp real(wp),parameter :: a1711 = 1.178022481753467227246517905742710665540e-2_wp real(wp),parameter :: a1712 = 2.689680149305807469351023586518679842111e-2_wp real(wp),parameter :: a1714 = -3.169098279840923653089425306993753388912e-2_wp real(wp),parameter :: a1715 = -3.145101260582300176675067993363288100508e-1_wp real(wp),parameter :: a1716 = 3.190394571030646408404948791505956372215e-1_wp real(wp),parameter :: c18 = 4.100000000000000000000000000000000000000e-1_wp real(wp),parameter :: a1801 = 1.450557047543664677045600699283886558262e-2_wp real(wp),parameter :: a1806 = 3.144536140911054963706342485304596004512e-2_wp real(wp),parameter :: a1807 = 3.715783248657690545268699075721576544859e-2_wp real(wp),parameter :: a1808 = -2.359517157110257865599846716538230421773e-2_wp real(wp),parameter :: a1809 = 3.920705621082224068192042184136974996456e-2_wp real(wp),parameter :: a1810 = -1.098622165458479216169190236380897121577e-2_wp real(wp),parameter :: a1811 = -8.096822446841331912908730838423549719449e-4_wp real(wp),parameter :: a1812 = -2.163269701613257487387923735849554667750e-3_wp real(wp),parameter :: a1814 = 2.654521616460219475066757083951961376948e-3_wp real(wp),parameter :: a1815 = 2.511562219208545046642238771986242777615e-1_wp real(wp),parameter :: a1816 = -6.498542101960180660630822424200730178886e-2_wp real(wp),parameter :: a1817 = 1.364132020723255014212599118638439066828e-1_wp real(wp),parameter :: c19 = 4.000000000000000000000000000000000000000e-1_wp real(wp),parameter :: a1901 = 1.446336618471592605799637355230438657926e-2_wp real(wp),parameter :: a1906 = 3.104751193922325663097227478535346108288e-2_wp real(wp),parameter :: a1907 = 3.661908805398219670465188920153386551595e-2_wp real(wp),parameter :: a1908 = -2.339478958466827786855261195290634759347e-2_wp real(wp),parameter :: a1909 = 3.879150732565195534576788590033483049488e-2_wp real(wp),parameter :: a1910 = -1.090859397382524594034743083379617423642e-2_wp real(wp),parameter :: a1911 = -8.025266815048055198422264011084742961702e-4_wp real(wp),parameter :: a1912 = -2.135231945628697281987752929223855801666e-3_wp real(wp),parameter :: a1914 = 2.622796236389624857351788846906102046025e-3_wp real(wp),parameter :: a1915 = 2.534442274803626704614239908457133748025e-1_wp real(wp),parameter :: a1916 = -7.643511000478099461439026575399944229462e-2_wp real(wp),parameter :: a1917 = 1.366877549700823911669560847388882737008e-1_wp real(wp),parameter :: c20 = 3.000000000000000000000000000000000000000e-1_wp real(wp),parameter :: a2001 = 1.454820927990212743631839742362023495019e-2_wp real(wp),parameter :: a2006 = 3.295833574477986148501962572781101218370e-2_wp real(wp),parameter :: a2007 = 3.908142787627776032656206482454472464491e-2_wp real(wp),parameter :: a2008 = -2.453625117065831858258113695738948397469e-2_wp real(wp),parameter :: a2009 = 4.093420165107695002931056233872486867656e-2_wp real(wp),parameter :: a2010 = -1.139339116697368269416077973418873431685e-2_wp real(wp),parameter :: a2011 = -8.425286304615995039446748728711864702110e-4_wp real(wp),parameter :: a2012 = -2.268675533978017976128644496915957259603e-3_wp real(wp),parameter :: a2014 = 2.778258636762732394270017683630499064869e-3_wp real(wp),parameter :: a2015 = 2.001138381588676427760578503039957001874e-1_wp real(wp),parameter :: a2016 = -1.271353237985016727288823264896836477597e-1_wp real(wp),parameter :: a2017 = 1.357618989529062170381590442487219700734e-1_wp real(wp),parameter :: c21 = 6.700000000000000000000000000000000000000e-1_wp real(wp),parameter :: a2101 = 6.617091212153454764237535387701483714178e-2_wp real(wp),parameter :: a2106 = 1.799282440576909948036063074652655951235e-1_wp real(wp),parameter :: a2107 = 3.624950045318951547013441289753072027999e-1_wp real(wp),parameter :: a2108 = 7.936089314877400356470244958243950579864e-2_wp real(wp),parameter :: a2109 = 4.850720370946592974218363233076530576227e-2_wp real(wp),parameter :: a2110 = 7.118400389001109880604226157877318948773e-2_wp real(wp),parameter :: a2111 = 2.113188147405872759613334579630890266022e-3_wp real(wp),parameter :: a2112 = -1.383686274944340964500317887611596828088e-2_wp real(wp),parameter :: a2114 = 1.258226100918847618852170065178443441806e-2_wp real(wp),parameter :: a2115 = -1.802109810437443834350265465118022589859e-1_wp real(wp),parameter :: a2116 = 1.433584739547597910202038973213568278251e-1_wp real(wp),parameter :: a2117 = -1.016523407775380761485633409744195613562e-1_wp associate (t => me%dense_t, h => me%dense_h, x => me%dense_x, k => me%funcs) ! f at the end of the step (this is also the first stage of the next step, ! so it is saved in the function evaluation cache): call me%set_fsal_cache(t+h, me%dense_xf, k(:,14)) if (me%stopped) return call me%f(t+c15*h, x+h*(a1501*k(:,1) + a1506*k(:,6) + a1507*k(:,7) + a1508*k(:,8) + & a1509*k(:,9) + a1510*k(:,10) + a1511*k(:,11) + a1512*k(:,12) + & a1514*k(:,14)), k(:,15)) if (me%stopped) return call me%f(t+c16*h, x+h*(a1601*k(:,1) + a1606*k(:,6) + a1607*k(:,7) + a1608*k(:,8) + & a1609*k(:,9) + a1610*k(:,10) + a1611*k(:,11) + a1612*k(:,12) + & a1614*k(:,14) + a1615*k(:,15)), k(:,16)) if (me%stopped) return call me%f(t+c17*h, x+h*(a1701*k(:,1) + a1706*k(:,6) + a1707*k(:,7) + a1708*k(:,8) + & a1709*k(:,9) + a1710*k(:,10) + a1711*k(:,11) + a1712*k(:,12) + & a1714*k(:,14) + a1715*k(:,15) + a1716*k(:,16)), k(:,17)) if (me%stopped) return call me%f(t+c18*h, x+h*(a1801*k(:,1) + a1806*k(:,6) + a1807*k(:,7) + a1808*k(:,8) + & a1809*k(:,9) + a1810*k(:,10) + a1811*k(:,11) + a1812*k(:,12) + & a1814*k(:,14) + a1815*k(:,15) + a1816*k(:,16) + a1817*k(:,17)), k(:,18)) if (me%stopped) return call me%f(t+c19*h, x+h*(a1901*k(:,1) + a1906*k(:,6) + a1907*k(:,7) + a1908*k(:,8) + & a1909*k(:,9) + a1910*k(:,10) + a1911*k(:,11) + a1912*k(:,12) + & a1914*k(:,14) + a1915*k(:,15) + a1916*k(:,16) + a1917*k(:,17)), k(:,19)) if (me%stopped) return call me%f(t+c20*h, x+h*(a2001*k(:,1) + a2006*k(:,6) + a2007*k(:,7) + a2008*k(:,8) + & a2009*k(:,9) + a2010*k(:,10) + a2011*k(:,11) + a2012*k(:,12) + & a2014*k(:,14) + a2015*k(:,15) + a2016*k(:,16) + a2017*k(:,17)), k(:,20)) if (me%stopped) return call me%f(t+c21*h, x+h*(a2101*k(:,1) + a2106*k(:,6) + a2107*k(:,7) + a2108*k(:,8) + & a2109*k(:,9) + a2110*k(:,10) + a2111*k(:,11) + a2112*k(:,12) + & a2114*k(:,14) + a2115*k(:,15) + a2116*k(:,16) + a2117*k(:,17)), k(:,21)) end associate end procedure rkdp87_dense_extra_stages !***************************************************************************************** !***************************************************************************************** !> ! Dense output for the Enright-Verner 8(7) method [[rkev87]]. ! This is an interpolant of order 8 constructed by P. Stone, which requires 8 additional ! function evaluations (see [[rkev87_dense_extra_stages]]). ! !@note The interpolant is published in the polynomial form. Here it has been ! converted (exactly, using high precision arithmetic) ! to the form used in `DOPRI5` and `DOP853` (see [[dense_coefficients]]), which has ! much smaller roundoff error for these coefficients. This code is generated by ! `scripts/generate_dense.py`. ! !### References ! * P. Stone, [Maple worksheets on interpolation for Runge-Kutta schemes](http://www.peterstone.name/Maplepgs/RKinterp.html) ! * Coefficients: ! * [RKinterp8b.mws](http://www.peterstone.name/Maplepgs/Maple/nmthds/RKinterp/RKinterp8b.mws) module procedure rkev87_dense_coefficients form = dense_form_dopri allocate(b(21,8)) b = 0.0_wp b(1,:) = [ 4.391770364439902647673901682456129382042e-2_wp, & 9.560822963556009735232609831754387061796e-1_wp, & -9.121645927112019470465219663508774123592e-1_wp, & -9.226068447724062761577407304039914142130e0_wp, & 5.598045717439134379318726813131900122754e0_wp, & 3.512103616983960144470662346810604936278e1_wp, & -1.025820796200233312492454603401904783481e1_wp, & -4.760577062696842016269814659493753203233e1_wp ] b(6,:) = [ 3.510246253011980590037119470671975224711e-1_wp, & -3.510246253011980590037119470671975224711e-1_wp, & 7.020492506023961180074238941343950449422e-1_wp, & -2.386212040023204584488861182995857139752e0_wp, & -2.849774847168109309335439075045528029279e1_wp, & 4.403889539371758169864774544999594799587e1_wp, & 7.885817390551639300969897476128707791243e1_wp, & -1.120786845632716657475330061879462015662e2_wp ] b(7,:) = [ 2.461428263549239789274962176113575612945e-1_wp, & -2.461428263549239789274962176113575612945e-1_wp, & 4.922856527098479578549924352227151225890e-1_wp, & -1.673241514920736564488663170200079000965e0_wp, & -1.998297511906028297539220361752425578749e1_wp, & 3.088062033385066725878472329907825448961e1_wp, & 5.529633081906083171660114632809091417263e1_wp, & -7.859096543119796808522095407071565620779e1_wp ] b(8,:) = [ 9.003244930529127822924857436630861679802e-1_wp, & -9.003244930529127822924857436630861679802e-1_wp, & 1.800648986105825564584971487326172335960e0_wp, & -6.120269036416562910916801156023926341804e0_wp, & -7.309236759073646446129677954298249017693e1_wp, & 1.129530332407243601263088284551277801784e2_wp, & 2.022591588371967230789133878855458133382e2_wp, & -2.874646893359150443113642867952308621369e2_wp ] b(9,:) = [ 4.549418727254746872059927459954873596623e0_wp, & -4.549418727254746872059927459954873596623e0_wp, & 9.098837454509493744119854919909747193245e0_wp, & -3.092625690510285981340028036736195085958e1_wp, & -3.693421521934997210686072931028055228916e2_wp, & 5.707615961697253868138814117428456224196e2_wp, & 1.022033291410918639945732308373736261728e3_wp, & -1.452584319520916942161729927734819696928e3_wp ] b(10,:) = [ 4.802501519237061312514354065953544905625e-3_wp, & -4.802501519237061312514354065953544905625e-3_wp, & 9.605003038474122625028708131907089811249e-3_wp, & -3.264667525134480850672911831119496330729e-2_wp, & -3.898885447499599745013492647927533261575e-1_wp, & 6.025128916593102094076203595025241112233e-1_wp, & 1.078888695231976050774468291546015106765e0_wp, & -1.533391147208971761445002780707571149517e0_wp ] b(11,:) = [ -4.741054352120046004585409588596098810320e0_wp, & 4.741054352120046004585409588596098810320e0_wp, & -9.482108704240092009170819177192197620640e0_wp, & 3.222896675048357438978769668158568445128e1_wp, & 3.848999890003583264638834028294779183413e2_wp, & -5.948038445729852261399932434009665032115e2_wp, & -1.065084503043587200187373623123000494432e3_wp, & 1.513771675627155382782458718369879581767e3_wp ] b(12,:) = [ -3.545765250073717754874651505909308767743e-1_wp, & 3.545765250073717754874651505909308767743e-1_wp, & -7.091530500147435509749303011818617535486e-1_wp, & 2.410357314265870879555490251578106649710e0_wp, & 2.878610756995325377680442601483898625719e1_wp, & -4.448451011227170657363804529121904559592e1_wp, & -7.965611315118629938961956670143901227218e1_wp, & 1.132127709437557604642739023063636240156e2_wp ] b(14,:) = [ 0.0_wp, & 0.0_wp, & -1.000000000000000000000000000000000000000e0_wp, & 1.290564266715399738082418597337832257123e0_wp, & 9.154258302172844046720482707479283583241e0_wp, & -1.185046187603745209794222346365454571997e1_wp, & -2.032875139086708579896907289611208924741e1_wp, & 2.606009301180272428529884952500840836902e1_wp ] b(15,:) = [ 0.0_wp, & 0.0_wp, & 0.0_wp, & 4.231470307474816277003914433432315531246e1_wp, & -1.365909559879506945681075473158055923938e-1_wp, & -2.587527159683363886327380231471380299930e2_wp, & -2.242081699568125039450868637819902574330e0_wp, & 4.492522422300282664511058083618432794105e2_wp ] b(16,:) = [ 0.0_wp, & 0.0_wp, & 0.0_wp, & -5.381990108699555475713723369386504921274e1_wp, & 1.737294889595162571935015793088652540454e-1_wp, & 3.291065414025760813220995135133517997745e2_wp, & 2.851694719127780067552326058315612506541e0_wp, & -5.714021246284023202861480742839820204234e2_wp ] b(17,:) = [ 0.0_wp, & 0.0_wp, & 0.0_wp, & 2.057792964324803569396609715006582088157e1_wp, & -6.642511651940184856130026382404824051241e-2_wp, & -1.258332162886731592441692070591278668710e2_wp, & -1.090339671925050958430284899255856226265e0_wp, & 2.184744394011316284826723830925058283502e2_wp ] b(18,:) = [ 0.0_wp, & 0.0_wp, & 0.0_wp, & 1.186424657408506895604566858468413103974e2_wp, & 1.998034884356286106313307838203909604296e1_wp, & -1.101022354780388947346820039638589292049e3_wp, & -1.668848919591184186069128554238276960415e2_wp, & 2.737591391703951336030846574911265631682e3_wp ] b(19,:) = [ 0.0_wp, & 0.0_wp, & 0.0_wp, & -1.195553183477307236415538581084511096651e1_wp, & -1.031725451800117816836910206812571966188e1_wp, & 2.248007018675025913750944329499202214279e2_wp, & 1.493879142910939668624544064592706851096e2_wp, & -8.538000374460980257653173793472978320075e2_wp ] b(20,:) = [ 0.0_wp, & 0.0_wp, & 0.0_wp, & -1.046952830989712272749686961517673313314e2_wp, & 1.960307395153888292922862065260665531446e1_wp, & 8.650688719889626051763068244353878261590e2_wp, & -3.729954100468316473836072900687708193664e1_wp, & -1.892755582554576129625018454853954953334e3_wp ] b(21,:) = [ 0.0_wp, & 0.0_wp, & 0.0_wp, & 3.370423849866892807752525093678504008655e0_wp, & 3.362984963625123336836828717894317105382e1_wp, & -7.658670585986530538993694167262074247840e1_wp, & -1.289210227952086328876854714354411993081e2_wp, & 2.394529523367303894098189960827259721917e2_wp ] end procedure rkev87_dense_coefficients !***************************************************************************************** !***************************************************************************************** !> ! The additional function evaluations for the dense output of [[rkev87]]. ! These are stored in `me%funcs(:,14:21)`. Stage 14 is `f` at the end of the step, ! which is cached for the first stage of the next step (see [[set_fsal_cache]]). module procedure rkev87_dense_extra_stages real(wp),parameter :: c15 = 3.110177634953863863927417318829099695921e-1_wp real(wp),parameter :: a1501 = 4.606118229177062378812330141508378303081e-2_wp real(wp),parameter :: a1506 = 5.563879265037094774988904887741689318457e-2_wp real(wp),parameter :: a1507 = 2.313778722397200205998710838964445918652e-1_wp real(wp),parameter :: a1508 = 2.774634770079758598408348133791100318753e-1_wp real(wp),parameter :: a1509 = 1.367534288907514327272420159245793697941e0_wp real(wp),parameter :: a1510 = -8.852053288949034134775680261299163362189e-2_wp real(wp),parameter :: a1511 = -1.470142317006875734859984437074622346224e0_wp real(wp),parameter :: a1512 = -9.579996416173765396824511059812662683433e-2_wp real(wp),parameter :: a1514 = -1.259503554386166268241032464519842162533e-2_wp real(wp),parameter :: c16 = 3.970000000000000000000000000000000000000e-1_wp real(wp),parameter :: a1601 = 4.744923181087010987698361118302172988759e-2_wp real(wp),parameter :: a1606 = 6.110556556028856636616048593054552986600e-2_wp real(wp),parameter :: a1607 = 2.204313596130638110210525058989630030785e-1_wp real(wp),parameter :: a1608 = 1.448749862745162472671772888868297968950e-1_wp real(wp),parameter :: a1609 = 8.256944261409224879929084716709217044802e-1_wp real(wp),parameter :: a1610 = -4.819911723269977121603368554900968489895e-2_wp real(wp),parameter :: a1611 = -8.635144982398826686314898674660474217939e-1_wp real(wp),parameter :: a1612 = -6.570900319522628662585304753751592460333e-2_wp real(wp),parameter :: a1614 = -8.030823824924940691410543552334783893484e-3_wp real(wp),parameter :: a1615 = 8.289787309307244464050478053462605098246e-2_wp real(wp),parameter :: c17 = 1.000000000000000000000000000000000000000e-1_wp real(wp),parameter :: a1701 = 4.890422422590243735643355928726905013253e-2_wp real(wp),parameter :: a1706 = 2.989337764810148368025189802500811497661e-2_wp real(wp),parameter :: a1707 = 8.010975565601287534176031826276361729309e-2_wp real(wp),parameter :: a1708 = -1.855925041794188436352098711408322328964e-1_wp real(wp),parameter :: a1709 = -5.672695518352900753386893208551888258026e-1_wp real(wp),parameter :: a1710 = 5.591411903664724506777966059129269423962e-2_wp real(wp),parameter :: a1711 = 6.837022871826681925244868001102237923926e-1_wp real(wp),parameter :: a1712 = 1.623428693609802296096003887418509596122e-2_wp real(wp),parameter :: a1714 = 3.488973286951633835819828802381050928671e-3_wp real(wp),parameter :: a1715 = -2.372012865046302651545935537589887956462e-2_wp real(wp),parameter :: a1716 = -4.166483930720994527813355658120347766068e-2_wp real(wp),parameter :: c18 = 3.630000000000000000000000000000000000000e-1_wp real(wp),parameter :: a1801 = 2.830610977564639786013880171591726169634e-2_wp real(wp),parameter :: a1806 = 1.808140586894492018873685554113617873645e-2_wp real(wp),parameter :: a1807 = 2.373848424564159475407909712024445350880e-2_wp real(wp),parameter :: a1808 = -5.633044498884479334388672069580263943199e-3_wp real(wp),parameter :: a1809 = -9.052690834203304458197301691952427041874e-2_wp real(wp),parameter :: a1810 = 6.455492970119906188733114964871342318416e-5_wp real(wp),parameter :: a1811 = 8.532476541091549683466744998774272980888e-2_wp real(wp),parameter :: a1812 = 9.617354401809802968006874345767582430442e-3_wp real(wp),parameter :: a1814 = 1.445000263984237186278121995021307415604e-3_wp real(wp),parameter :: a1815 = 2.030393801641215936234049750722068729254e-1_wp real(wp),parameter :: a1816 = -6.208488038564734995553585386027329942649e-2_wp real(wp),parameter :: a1817 = 1.516277781657996313946980359216927338434e-1_wp real(wp),parameter :: c19 = 4.500000000000000000000000000000000000000e-1_wp real(wp),parameter :: a1901 = 2.856290781663390425643250364579208858791e-2_wp real(wp),parameter :: a1906 = 2.316257839791545293498901433471915365768e-2_wp real(wp),parameter :: a1907 = 3.028654858021448917336337908830409463426e-2_wp real(wp),parameter :: a1908 = -6.638350434785318341883271931088849237251e-3_wp real(wp),parameter :: a1909 = -1.123580918926519086504223022775848960256e-1_wp real(wp),parameter :: a1910 = 8.472657524428672350650240464287981879645e-5_wp real(wp),parameter :: a1911 = 1.056422011075090480444438440399884595596e-1_wp real(wp),parameter :: a1912 = 1.201030302483773382796917358695352514706e-2_wp real(wp),parameter :: a1914 = 1.777469590426180785109808804304014856394e-3_wp real(wp),parameter :: a1915 = 1.913048458729429321787330256483826603242e-1_wp real(wp),parameter :: a1916 = 2.823832173555990836270193336222456579414e-2_wp real(wp),parameter :: a1917 = 1.479265396261532907050563892933623028828e-1_wp real(wp),parameter :: c20 = 3.030000000000000000000000000000000000000e-1_wp real(wp),parameter :: a2001 = 2.829184679280092973425764200351982457640e-2_wp real(wp),parameter :: a2006 = 2.010976600316258331802409429117206044189e-2_wp real(wp),parameter :: a2007 = 2.638412296846209109225112925829174916105e-2_wp real(wp),parameter :: a2008 = -6.183465183025303299354154885240659615577e-3_wp real(wp),parameter :: a2009 = -1.001731352792896258893736304854090196607e-1_wp real(wp),parameter :: a2010 = 7.208311650232223117809837798999916128577e-5_wp real(wp),parameter :: a2011 = 9.438010364821978038925968927657732098424e-2_wp real(wp),parameter :: a2012 = 1.065253633263355242041179222164602436088e-2_wp real(wp),parameter :: a2014 = 1.596590333948665651173908987745709127765e-3_wp real(wp),parameter :: a2015 = 1.560277402099908222947533850106655013358e-1_wp real(wp),parameter :: a2016 = -7.882660737994427001882445556948162584414e-2_wp real(wp),parameter :: a2017 = 1.506684184365384520762425015125231159711e-1_wp real(wp),parameter :: c21 = 6.800000000000000000000000000000000000000e-1_wp real(wp),parameter :: a2101 = 5.860772551157834546429464002350799906243e-2_wp real(wp),parameter :: a2106 = 3.477925728033368984466843926118366543751e-1_wp real(wp),parameter :: a2107 = 4.307844812588595385347441086505540713238e-1_wp real(wp),parameter :: a2108 = 1.307584239578842166335279530686311373068e-2_wp real(wp),parameter :: a2109 = -9.827916500818530921904784134255636922003e-1_wp real(wp),parameter :: a2110 = 1.668543705589422508453335285587334832077e-3_wp real(wp),parameter :: a2111 = 8.718284019203299084528640634435306456894e-1_wp real(wp),parameter :: a2112 = 1.198921810070596855514963392810281132805e-1_wp real(wp),parameter :: a2114 = 1.362109857187728652190264565068064962825e-2_wp real(wp),parameter :: a2115 = -1.533654673545670480346334054504936107439e-1_wp real(wp),parameter :: a2116 = 9.141568504729135645965614521691441835292e-2_wp real(wp),parameter :: a2117 = -1.325294147852907233783366465944456973310e-1_wp associate (t => me%dense_t, h => me%dense_h, x => me%dense_x, k => me%funcs) ! f at the end of the step (this is also the first stage of the next step, ! so it is saved in the function evaluation cache): call me%set_fsal_cache(t+h, me%dense_xf, k(:,14)) if (me%stopped) return call me%f(t+c15*h, x+h*(a1501*k(:,1) + a1506*k(:,6) + a1507*k(:,7) + a1508*k(:,8) + & a1509*k(:,9) + a1510*k(:,10) + a1511*k(:,11) + a1512*k(:,12) + & a1514*k(:,14)), k(:,15)) if (me%stopped) return call me%f(t+c16*h, x+h*(a1601*k(:,1) + a1606*k(:,6) + a1607*k(:,7) + a1608*k(:,8) + & a1609*k(:,9) + a1610*k(:,10) + a1611*k(:,11) + a1612*k(:,12) + & a1614*k(:,14) + a1615*k(:,15)), k(:,16)) if (me%stopped) return call me%f(t+c17*h, x+h*(a1701*k(:,1) + a1706*k(:,6) + a1707*k(:,7) + a1708*k(:,8) + & a1709*k(:,9) + a1710*k(:,10) + a1711*k(:,11) + a1712*k(:,12) + & a1714*k(:,14) + a1715*k(:,15) + a1716*k(:,16)), k(:,17)) if (me%stopped) return call me%f(t+c18*h, x+h*(a1801*k(:,1) + a1806*k(:,6) + a1807*k(:,7) + a1808*k(:,8) + & a1809*k(:,9) + a1810*k(:,10) + a1811*k(:,11) + a1812*k(:,12) + & a1814*k(:,14) + a1815*k(:,15) + a1816*k(:,16) + a1817*k(:,17)), k(:,18)) if (me%stopped) return call me%f(t+c19*h, x+h*(a1901*k(:,1) + a1906*k(:,6) + a1907*k(:,7) + a1908*k(:,8) + & a1909*k(:,9) + a1910*k(:,10) + a1911*k(:,11) + a1912*k(:,12) + & a1914*k(:,14) + a1915*k(:,15) + a1916*k(:,16) + a1917*k(:,17)), k(:,19)) if (me%stopped) return call me%f(t+c20*h, x+h*(a2001*k(:,1) + a2006*k(:,6) + a2007*k(:,7) + a2008*k(:,8) + & a2009*k(:,9) + a2010*k(:,10) + a2011*k(:,11) + a2012*k(:,12) + & a2014*k(:,14) + a2015*k(:,15) + a2016*k(:,16) + a2017*k(:,17)), k(:,20)) if (me%stopped) return call me%f(t+c21*h, x+h*(a2101*k(:,1) + a2106*k(:,6) + a2107*k(:,7) + a2108*k(:,8) + & a2109*k(:,9) + a2110*k(:,10) + a2111*k(:,11) + a2112*k(:,12) + & a2114*k(:,14) + a2115*k(:,15) + a2116*k(:,16) + a2117*k(:,17)), k(:,21)) end associate end procedure rkev87_dense_extra_stages !***************************************************************************************** !***************************************************************************************** !> ! Dense output for the Verner 8(7) (Maple's `dverk78`) method [[dverk78]]. ! This is an interpolant of order 8 constructed by P. Stone, which requires 9 additional ! function evaluations (see [[dverk78_dense_extra_stages]]). ! !@note The interpolant is published in the polynomial form. Here it has been ! converted (exactly, using high precision arithmetic) ! to the form used in `DOPRI5` and `DOP853` (see [[dense_coefficients]]), which has ! much smaller roundoff error for these coefficients. This code is generated by ! `scripts/generate_dense.py`. ! !### References ! * P. Stone, [Maple worksheets on interpolation for Runge-Kutta schemes](http://www.peterstone.name/Maplepgs/RKinterp.html) ! * Coefficients: ! * [RKinterp8a.mws](http://www.peterstone.name/Maplepgs/Maple/nmthds/RKinterp/RKinterp8a.mws) module procedure dverk78_dense_coefficients form = dense_form_dopri allocate(b(22,8)) b = 0.0_wp b(1,:) = [ 4.462884896530702174156131408819200715231e-2_wp, & 9.553711510346929782584386859118079928477e-1_wp, & -9.107423020693859565168773718236159856954e-1_wp, & -1.082197010114616359188532900231586553513e1_wp, & 7.736658447990972068809008947424705735139e0_wp, & 4.941422065330257962747888219179071047121e1_wp, & -1.805004097551440852641518391865110295381e1_wp, & -7.753378150895846850967062594800837606696e1_wp ] b(6,:) = [ 3.538496728080569550060787149620213259966e-1_wp, & -3.538496728080569550060787149620213259966e-1_wp, & 7.076993456161139100121574299240426519932e-1_wp, & -4.844350275725777128764188699117446501218e1_wp, & 1.589357941562891986628742302887865676704e1_wp, & 4.602814690655731127702895267743343053226e2_wp, & -7.625473002874366873046596745356060447554e1_wp, & -1.012486931337643552946247671045000438788e3_wp ] b(7,:) = [ 2.492141378240346531574135224332454745769e-1_wp, & -2.492141378240346531574135224332454745769e-1_wp, & 4.984282756480693063148270448664909491538e-1_wp, & -3.411845961879705334653287668727635111075e1_wp, & 1.119374976263535558572290489125791483859e1_wp, & 3.241733941966358940617213735287871331907e2_wp, & -5.370573511714484738197485429932833535046e1_wp, & -7.130882887329551056028527695282593843953e2_wp ] b(8,:) = [ -5.228338470980250040474048151553375184780e0_wp, & 5.228338470980250040474048151553375184780e0_wp, & -1.045667694196050008094809630310675036956e1_wp, & 7.157814422289939277115560801444405027685e2_wp, & -2.348370482891129334826385624189881631000e2_wp, & -6.800931291238562901286914722666618978384e3_wp, & 1.126708795403514045922626648762031995776e3_wp, & 1.496009402091160531384367705108919569505e4_wp ] b(9,:) = [ 5.151165492368612233259287600784658428717e0_wp, & -5.151165492368612233259287600784658428717e0_wp, & 1.030233098473722446651857520156931685743e1_wp, & -7.052161381197921431272629967524427540471e2_wp, & 2.313707320577848552342049255578059508221e2_wp, & 6.700546029650942455438126156199804036109e3_wp, & -1.110077991133315899967494986587494604682e3_wp, & -1.473927529956982416497053282646521667825e4_wp ] b(10,:) = [ -9.509792872164367300522532985563737470781e0_wp, & 9.509792872164367300522532985563737470781e0_wp, & -1.901958574432873460104506597112747494156e1_wp, & 1.301930488073508350211590372813365899873e3_wp, & -4.271436710411018781842358280680703099136e2_wp, & -1.237017233610245403917521565552938754496e4_wp, & 2.049363737830911592433216818761072707923e3_wp, & 2.721082352962136029406774986974110196408e4_wp ] b(11,:) = [ 1.005963413734314886385260052473208136948e1_wp, & -1.005963413734314886385260052473208136948e1_wp, & 2.011926827468629772770520104946416273896e1_wp, & -1.377206060986616648408629046120447012167e3_wp, & 4.518404462133350294998469999243452919740e2_wp, & 1.308539624257373452372027857895946061037e4_wp, & -2.167854725549344120159092828700124450662e3_wp, & -2.878411054409220189892573937955400446630e4_wp ] b(12,:) = [ -1.203609461645423860203605398830859503626e-1_wp, & 1.203609461645423860203605398830859503626e-1_wp, & -2.407218923290847720407210797661719007252e-1_wp, & 1.647791781497841025124209203233325264189e1_wp, & -5.406155221864696634920207316710844556330e0_wp, & -1.565634148509985640917918245116563187525e2_wp, & 2.593782660005429990129218800947025635179e1_wp, & 3.443945110022384281279311761500223718941e2_wp ] b(14,:) = [ 0.0_wp, & 0.0_wp, & -1.000000000000000000000000000000000000000e0_wp, & 5.088704054024093087799387668553443594206e0_wp, & 5.994582015478069564188573832816516408376e0_wp, & -4.667692215748528524596123726804303132185e1_wp, & -9.498726354527230029495717627276903990294e0_wp, & 1.024315875585863189658699845689645976399e2_wp ] b(18,:) = [ 0.0_wp, & 0.0_wp, & 0.0_wp, & 3.131789498273708573049730375632215461001e1_wp, & 1.846213917235432944565002829297699439341e1_wp, & -3.734468736000289229193968548497487348188e2_wp, & -1.391073488646631981199410207629536854981e2_wp, & 1.169734934755895686119858763016208356237e3_wp ] b(19,:) = [ 0.0_wp, & 0.0_wp, & 0.0_wp, & 1.308173556520692942353495482071898398753e1_wp, & -3.339310337035622457016455777318149517476e1_wp, & -6.488497100548994914179896811147627319391e1_wp, & 2.367615505507669069308954819364420005148e2_wp, & -2.080479427566821734337347805616399683058e2_wp ] b(20,:) = [ 0.0_wp, & 0.0_wp, & 0.0_wp, & 3.291747252196333540086079847020008036281e1_wp, & 8.110865450027851548852633741504771718145e-1_wp, & -3.365583044113196321287789454150256348773e2_wp, & -1.355501505795424590164958408048044037892e1_wp, & 7.935585352516262851477429979844373387833e2_wp ] b(21,:) = [ 0.0_wp, & 0.0_wp, & 0.0_wp, & 1.943585266543131996319679094947528043748e1_wp, & -2.112087474603707162670707481870812558225e1_wp, & -1.732207076857579426268930799391711977986e2_wp, & 7.456987635802592241461547505058943153922e1_wp, & 3.700613874286637534556070828092263962593e2_wp ] b(22,:) = [ 0.0_wp, & 0.0_wp, & 0.0_wp, & 3.977462367676632798167435489824684959602e1_wp, & -2.140212096173751192092889745399756978355e1_wp, & -2.973565350880913290011432293630490813576e2_wp, & 7.476252633793485121388353091026373588585e1_wp, & 5.834442814682892846603411277429725921595e2_wp ] end procedure dverk78_dense_coefficients !***************************************************************************************** !***************************************************************************************** !> ! The additional function evaluations for the dense output of [[dverk78]]. ! These are stored in `me%funcs(:,14:22)`. Stage 14 is `f` at the end of the step, ! which is cached for the first stage of the next step (see [[set_fsal_cache]]). module procedure dverk78_dense_extra_stages real(wp),parameter :: c15 = 3.110177634953863863927417318829099695921e-1_wp real(wp),parameter :: a1501 = 4.584388949209796115333574616962239112076e-2_wp real(wp),parameter :: a1506 = 3.630388961716748777014844642836338085588e-2_wp real(wp),parameter :: a1507 = 2.368382446163510984791462973731356347065e-1_wp real(wp),parameter :: a1508 = 4.485575209861175949930620094528135096786e-2_wp real(wp),parameter :: a1509 = -7.321174242409605925313660599505542337546e-2_wp real(wp),parameter :: a1510 = 1.116796392388900052858067478262523693343e-2_wp real(wp),parameter :: a1511 = 5.433408462262251996481696132791993361756e-2_wp real(wp),parameter :: a1512 = -3.251928290739571906704566450378411360910e-2_wp real(wp),parameter :: a1514 = -1.259503554386166268241032464519842162533e-2_wp real(wp),parameter :: c16 = 4.748000000000000000000000000000000000000e-1_wp real(wp),parameter :: a1601 = 5.013938388149897679433589683705420340682e-2_wp real(wp),parameter :: a1606 = 9.141586214121856346416973884712596700136e-2_wp real(wp),parameter :: a1607 = 2.091392497721726472202048080499861142508e-1_wp real(wp),parameter :: a1608 = -2.400012712906264187646533945836691022068e0_wp real(wp),parameter :: a1609 = 8.900134418750022717562231888853548439473e-1_wp real(wp),parameter :: a1610 = -2.245534956015498921270900517229655156380e0_wp real(wp),parameter :: a1611 = 3.746489998397286616147902470749466758309e0_wp real(wp),parameter :: a1612 = 4.031671151895576459870187466249515023204e-3_wp real(wp),parameter :: a1614 = 3.432993897133524674027806705287275818332e-3_wp real(wp),parameter :: a1615 = 1.256850678055549324007003655258215006911e-1_wp real(wp),parameter :: c17 = 1.108000000000000000000000000000000000000e-1_wp real(wp),parameter :: a1701 = 4.922008100217028714063404492454142677427e-2_wp real(wp),parameter :: a1706 = 1.284497392104008890590815978453353497457e-2_wp real(wp),parameter :: a1707 = 1.038723293704250078931641665575346450846e-1_wp real(wp),parameter :: a1708 = -4.351647910071053869843388605538477801359e-1_wp real(wp),parameter :: a1709 = 1.369971819296699751972180408519310459054e-1_wp real(wp),parameter :: a1710 = -4.247661738893296126410657536853715193403e-1_wp real(wp),parameter :: a1711 = 7.366727838070516174761904151140862540873e-1_wp real(wp),parameter :: a1712 = -1.358157893841445902942840010049297027193e-2_wp real(wp),parameter :: a1714 = -4.716733873244034091777040578517197856614e-3_wp real(wp),parameter :: a1715 = -6.901019526658074709850194004986309652592e-2_wp real(wp),parameter :: a1716 = 1.843212294431726323199716773546565730450e-2_wp real(wp),parameter :: c18 = 3.972000000000000000000000000000000000000e-1_wp real(wp),parameter :: a1801 = 3.290668231387221106499148568686740782308e-2_wp real(wp),parameter :: a1806 = 1.288934041390745904508465192821521112783e-2_wp real(wp),parameter :: a1807 = 1.698357377896843092639739824390708045600e-2_wp real(wp),parameter :: a1808 = 7.292864046310047944624850925409755668153e-2_wp real(wp),parameter :: a1809 = -3.622539061772744485854767044823592754575e-2_wp real(wp),parameter :: a1810 = 9.057043432999283980529409146494162871781e-2_wp real(wp),parameter :: a1811 = -1.300470488380912846496984279896433146661e-1_wp real(wp),parameter :: a1812 = 2.170488422343271834358917939739442800179e-3_wp real(wp),parameter :: a1814 = 9.123635805634085083277263532800829734771e-4_wp real(wp),parameter :: a1815 = 1.805148604892410096659535546269709142551e-1_wp real(wp),parameter :: a1816 = -8.769595306265018517208922617874235190155e-3_wp real(wp),parameter :: a1817 = 1.623656509700946377287986855577341525670e-1_wp real(wp),parameter :: c19 = 4.530000000000000000000000000000000000000e-1_wp real(wp),parameter :: a1901 = 3.199928466558771236920747252814410784540e-2_wp real(wp),parameter :: a1906 = -6.661246984901002866428709970539443432525e-3_wp real(wp),parameter :: a1907 = -8.450986710613110627019614990764027078245e-3_wp real(wp),parameter :: a1908 = -2.682339678507348581157101040131459940270e-2_wp real(wp),parameter :: a1909 = 9.485330779997651557596630551518064183624e-3_wp real(wp),parameter :: a1910 = -2.877844401616755913542024429172333753057e-2_wp real(wp),parameter :: a1911 = 4.672501620339446429477702462622362090388e-2_wp real(wp),parameter :: a1912 = -8.512794192876897223347342801936860355256e-4_wp real(wp),parameter :: a1914 = -2.942053482228113873286779224488495461837e-4_wp real(wp),parameter :: a1915 = 2.130071642261069082282907135462123603660e-1_wp real(wp),parameter :: a1916 = 4.512850648175037408233536825426407144800e-2_wp real(wp),parameter :: a1917 = 1.785142569074285490178957823506217182789e-1_wp real(wp),parameter :: c20 = 7.116000000000000000000000000000000000000e-1_wp real(wp),parameter :: a2001 = 6.229218721188033740588795753796336349035e-2_wp real(wp),parameter :: a2006 = 4.200503869104953475235476119758027627861e-1_wp real(wp),parameter :: a2007 = 5.218421464442877145342989770677193715290e-1_wp real(wp),parameter :: a2008 = 1.322753273348260953269394879203970060808e0_wp real(wp),parameter :: a2009 = -2.847511151223089892095021724646754517499e-1_wp real(wp),parameter :: a2010 = 1.203017397955593527042485075749566065335e0_wp real(wp),parameter :: a2011 = -2.251406269982176583122823927442719901192e0_wp real(wp),parameter :: a2012 = 4.450471608613795047051990303470629205984e-2_wp real(wp),parameter :: a2014 = 1.387836911995496720875196385462012995343e-2_wp real(wp),parameter :: a2015 = -1.466949613936094135244062240083955466537e-1_wp real(wp),parameter :: a2016 = 1.313624120726726666229319137973471404051e-2_wp real(wp),parameter :: a2017 = -2.070223717857830782604472358882918604075e-1_wp real(wp),parameter :: c21 = 9.029000000000000000000000000000000000000e-1_wp real(wp),parameter :: a2101 = 6.793530908700785243117501353377251857208e-2_wp real(wp),parameter :: a2106 = 5.739044183428677105248917072976695116637e-1_wp real(wp),parameter :: a2107 = 5.916896612490674638858939079760143934695e-1_wp real(wp),parameter :: a2108 = -2.233513897385083542818668596682321258585e0_wp real(wp),parameter :: a2109 = 3.045475241181545348485496484842096100731e0_wp real(wp),parameter :: a2110 = -5.060476001779407556941970943180120519311e0_wp real(wp),parameter :: a2111 = 4.541025377366928347946716727239641273976e0_wp real(wp),parameter :: a2112 = -3.975844072812268650718639766464814966597e-2_wp real(wp),parameter :: a2114 = -1.903374112575350720783879200966692374787e-2_wp real(wp),parameter :: a2115 = -1.834637916089901458527429412276356981432e-1_wp real(wp),parameter :: a2116 = -1.177630612929609298883176380810346351253e-1_wp real(wp),parameter :: a2117 = -2.631210733070983540574485320437666138334e-1_wp real(wp),parameter :: c22 = 1.012000000000000000000000000000000000000e-1_wp real(wp),parameter :: a2201 = 3.704521498946937261201622638053173421237e-2_wp real(wp),parameter :: a2206 = -3.706525050079791723549371980886930577599e-2_wp real(wp),parameter :: a2207 = -4.910283885838155162949354115462944620334e-2_wp real(wp),parameter :: a2208 = -2.185124579362087257242571915004179965439e-1_wp real(wp),parameter :: a2209 = 1.116471738453870695461006698427219035348e-1_wp real(wp),parameter :: a2210 = -2.750411933592176438972612878935368766766e-1_wp real(wp),parameter :: a2211 = 3.905492741078667527053830032795101430151e-1_wp real(wp),parameter :: a2212 = -6.460455853466138764727561016065122104754e-3_wp real(wp),parameter :: a2214 = -2.806524006504184875752356195713271277510e-3_wp real(wp),parameter :: a2215 = -1.137565900192644457212514851621504314417e-3_wp real(wp),parameter :: a2216 = 3.873574615952291475757000276215831243043e-2_wp real(wp),parameter :: a2217 = 1.133488773125226969631282701559314297038e-1_wp associate (t => me%dense_t, h => me%dense_h, x => me%dense_x, k => me%funcs) ! f at the end of the step (this is also the first stage of the next step, ! so it is saved in the function evaluation cache): call me%set_fsal_cache(t+h, me%dense_xf, k(:,14)) if (me%stopped) return call me%f(t+c15*h, x+h*(a1501*k(:,1) + a1506*k(:,6) + a1507*k(:,7) + a1508*k(:,8) + & a1509*k(:,9) + a1510*k(:,10) + a1511*k(:,11) + a1512*k(:,12) + & a1514*k(:,14)), k(:,15)) if (me%stopped) return call me%f(t+c16*h, x+h*(a1601*k(:,1) + a1606*k(:,6) + a1607*k(:,7) + a1608*k(:,8) + & a1609*k(:,9) + a1610*k(:,10) + a1611*k(:,11) + a1612*k(:,12) + & a1614*k(:,14) + a1615*k(:,15)), k(:,16)) if (me%stopped) return call me%f(t+c17*h, x+h*(a1701*k(:,1) + a1706*k(:,6) + a1707*k(:,7) + a1708*k(:,8) + & a1709*k(:,9) + a1710*k(:,10) + a1711*k(:,11) + a1712*k(:,12) + & a1714*k(:,14) + a1715*k(:,15) + a1716*k(:,16)), k(:,17)) if (me%stopped) return call me%f(t+c18*h, x+h*(a1801*k(:,1) + a1806*k(:,6) + a1807*k(:,7) + a1808*k(:,8) + & a1809*k(:,9) + a1810*k(:,10) + a1811*k(:,11) + a1812*k(:,12) + & a1814*k(:,14) + a1815*k(:,15) + a1816*k(:,16) + a1817*k(:,17)), k(:,18)) if (me%stopped) return call me%f(t+c19*h, x+h*(a1901*k(:,1) + a1906*k(:,6) + a1907*k(:,7) + a1908*k(:,8) + & a1909*k(:,9) + a1910*k(:,10) + a1911*k(:,11) + a1912*k(:,12) + & a1914*k(:,14) + a1915*k(:,15) + a1916*k(:,16) + a1917*k(:,17)), k(:,19)) if (me%stopped) return call me%f(t+c20*h, x+h*(a2001*k(:,1) + a2006*k(:,6) + a2007*k(:,7) + a2008*k(:,8) + & a2009*k(:,9) + a2010*k(:,10) + a2011*k(:,11) + a2012*k(:,12) + & a2014*k(:,14) + a2015*k(:,15) + a2016*k(:,16) + a2017*k(:,17)), k(:,20)) if (me%stopped) return call me%f(t+c21*h, x+h*(a2101*k(:,1) + a2106*k(:,6) + a2107*k(:,7) + a2108*k(:,8) + & a2109*k(:,9) + a2110*k(:,10) + a2111*k(:,11) + a2112*k(:,12) + & a2114*k(:,14) + a2115*k(:,15) + a2116*k(:,16) + a2117*k(:,17)), k(:,21)) if (me%stopped) return call me%f(t+c22*h, x+h*(a2201*k(:,1) + a2206*k(:,6) + a2207*k(:,7) + a2208*k(:,8) + & a2209*k(:,9) + a2210*k(:,10) + a2211*k(:,11) + a2212*k(:,12) + & a2214*k(:,14) + a2215*k(:,15) + a2216*k(:,16) + a2217*k(:,17)), k(:,22)) end associate end procedure dverk78_dense_extra_stages !***************************************************************************************** ! END GENERATED DENSE OUTPUT !***************************************************************************************** end submodule rklib_variable_dense !*****************************************************************************************