########################################################################################
##
## EMBEDDED DIAGONALLY IMPLICIT RUNGE KUTTA METHOD
## (solvers/esdirk85.py)
##
## Milan Rother 2024
##
########################################################################################
# IMPORTS ==============================================================================
from ._rungekutta import DiagonallyImplicitRungeKutta
# SOLVERS ==============================================================================
[docs]
class ESDIRK85(DiagonallyImplicitRungeKutta):
"""Sixteen-stage, 8th order ESDIRK method with embedded 5th order error
estimate. L-stable and stiffly accurate (ESDIRK(16,8)[2]SAL-[(16,5)]).
Characteristics
---------------
* Order: 8 (propagating) / 5 (embedded)
* Stages: 16 (1 explicit, 15 implicit)
* Adaptive timestep
* L-stable, stiffly accurate
* Stage order 2
Note
----
Fifteen implicit solves per step make this very expensive. It is only
justified when the right-hand side evaluation is itself costly (large
state dimension, expensive ``ODE`` blocks) and very tight tolerances are
required so that the 8th order convergence compensates through much
larger steps. For generating stiff reference solutions to validate other
solvers. In almost all practical block-diagram simulations, ``ESDIRK54``
is the better choice.
References
----------
.. [1] Alamri, Y., & Ketcheson, D. I. (2024). "Very high-order A-stable
stiffly accurate diagonally implicit Runge-Kutta methods with
error estimators". Journal of Scientific Computing, 100,
Article 84. :doi:`10.1007/s10915-024-02627-w`
.. [2] Kennedy, C. A., & Carpenter, M. H. (2019). "Diagonally implicit
Runge-Kutta methods for stiff ODEs". Applied Numerical
Mathematics, 146, 221-244.
:doi:`10.1016/j.apnum.2019.07.008`
"""
def __init__(self, *solver_args, **solver_kwargs):
super().__init__(*solver_args, **solver_kwargs)
#number of stages in RK scheme
self.s = 16
#order of scheme and embedded method
self.n = 8
self.m = 5
#flag adaptive timestep solver
self.is_adaptive = True
#intermediate evaluation times as ratios
self.eval_stages = [
0.0 , 0.234637638717043, 0.558545926594724,
0.562667638694992, 0.697898381329126, 0.956146958839776,
0.812903043340468, 0.148256733818785, 0.944650387704291,
0.428471803715736, 0.984131639774509, 0.320412672954752,
0.974077670791771, 0.852850433853921, 0.823320301074444,
1.0
]
#butcher table
self.BT = {
0:None, #explicit first stage
1:[0.117318819358521, 0.117318819358521],
2:[0.0557014605974616, 0.385525646638742, 0.117318819358521],
3:[0.063493276428895, 0.373556126263681, 0.0082994166438953,
0.117318819358521],
4:[0.0961351856230088, 0.335558324517178, 0.207077765910132,
-0.0581917140797146, 0.117318819358521],
5:[0.0497669214238319, 0.384288616546039, 0.0821728117583936,
0.120337007107103, 0.202262782645888, 0.117318819358521],
6:[0.00626710666809847, 0.496491452640725, -0.111303249827358,
0.170478821683603, 0.166517073971103, -0.0328669811542241,
0.117318819358521],
7:[0.0463439767281591, 0.00306724391019652, -0.00816305222386205,
-0.0353302599538294, 0.0139313601702569, -0.00992014507967429,
0.0210087909090165, 0.117318819358521],
8:[0.111574049232048, 0.467639166482209, 0.237773114804619,
0.0798895699267508, 0.109580615914593, 0.0307353103825936,
-0.0404391509541147, -0.16942110744293, 0.117318819358521],
9:[-0.0107072484863877, -0.231376703354252, 0.017541113036611,
0.144871527682418, -0.041855459769806, 0.0841832168332261,
-0.0850020937282192, 0.486170343825899, -0.0526717116822739,
0.117318819358521],
10:[-0.0142238262314935, 0.14752923682514, 0.238235830732566,
0.037950291904103, 0.252075123381518, 0.0474266904224567,
-0.00363139069342027, 0.274081442388563, -0.0599166970745255,
-0.0527138812389185, 0.117318819358521],
11:[-0.11837020183211, -0.635712481821264, 0.239738832602538,
0.330058936651707, -0.325784087988237, -0.0506514314589253,
-0.281914404487009, 0.852596345144291, 0.651444614298805,
-0.103476387303591, -0.354835880209975, 0.117318819358521],
12:[-0.00458164025442349, 0.296219694015248, 0.322146049419995,
0.15917778285238, 0.284864871688843, 0.185509526463076,
-0.0784621067883274, 0.166312223692047, -0.284152486083397,
-0.357125104338944, 0.078437074055306, 0.0884129667114481,
0.117318819358521],
13:[-0.0545561913848106, 0.675785423442753, 0.423066443201941,
-0.000165300126841193, 0.104252994793763, -0.105763019303021,
-0.15988308809318, 0.0515050001032011, 0.56013979290924,
-0.45781539708603, -0.255870699752664, 0.026960254296416,
-0.0721245985053681, 0.117318819358521],
14:[0.0649253995775223, -0.0216056457922249, -0.073738139377975,
0.0931033310077225, -0.0194339577299149, -0.0879623837313009,
0.057125517179467, 0.205120850488097, 0.132576503537441,
0.489416890627328, -0.1106765720501, -0.081038793996096,
0.0606031613503788, -0.00241467937442272, 0.117318819358521],
15:[0.0459979286336779, 0.0780075394482806, 0.015021874148058,
0.195180277284195, -0.00246643310153235, 0.0473977117068314,
-0.0682773558610363, 0.19568019123878, -0.0876765449323747,
0.177874852409192, -0.337519251582222, -0.0123255553640736,
0.311573291192553, 0.0458604327754991, 0.278352222645651,
0.117318819358521]
}
#coefficients for truncation error estimate (8th and 5th order solution)
_A1 = [
0.0459979286336779, 0.0780075394482806, 0.015021874148058,
0.195180277284195, -0.00246643310153235, 0.0473977117068314,
-0.0682773558610363, 0.19568019123878, -0.0876765449323747,
0.177874852409192, -0.337519251582222, -0.0123255553640736,
0.311573291192553, 0.0458604327754991, 0.278352222645651,
0.117318819358521
]
_A2 = [
0.0603373529853206, 0.175453809423998, 0.0537707777611352,
0.195309248607308, 0.0135893741970232, -0.0221160259296707,
-0.00726526156430691, 0.102961059369124, 0.000900215457460583,
0.0547959465692338, -0.334995726863153, 0.0464409662093384,
0.301388101652194, 0.00524851570622031, 0.229538601845236,
0.124643044573514
]
self.TR = [_a1 - _a2 for _a1, _a2 in zip(_A1, _A2)]