########################################################################################
##
## 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 L-stable Embedded Singly Diagonally Implicit Runge-Kutta method.
Features an embedded 5th order method for adaptive step size control. The first stage is
explicit. Designed for very stiff problems requiring very high accuracy. Computationally
expensive due to 16 stages, but can take very large timesteps with tight tolerances.
This is the ESDIRK(16,8)[2]SAL-[(16,5)] method.
Characteristics
---------------
* Order: 8
* Embedded Order: 5
* Stages: 16 (1 Explicit, 15 Implicit)
* Implicit (ESDIRK)
* Adaptive timestep
* L-stable, Stiffly accurate
When to Use
-----------
* **Extremely high accuracy on stiff problems**: When very tight tolerances are essential
* **Expensive right-hand sides**: When large timesteps justify the 16 stages
* **Benchmark computations**: Reference solutions for stiff problems
* **Specialized applications**: Only when 8th order accuracy is truly needed
**Warning**: Very expensive (16 implicit stages). Only use when extremely high accuracy
is essential and large timesteps are possible. For most applications, ESDIRK54 is more
practical.
References
----------
.. [1] Alamri, Y., & Ketcheson, D. I. (2019). "Very high-order A-stable stiffly accurate
diagonally implicit Runge-Kutta methods with error estimators". arXiv preprint
arXiv:1905.11370.
.. [2] Kennedy, C. A., & Carpenter, M. H. (2019). "Diagonally implicit Runge-Kutta
methods for stiff ODEs". Applied Numerical Mathematics, 146, 221-244.
.. [3] Hairer, E., & Wanner, G. (1996). "Solving Ordinary Differential Equations II:
Stiff and Differential-Algebraic Problems". Springer Series in Computational
Mathematics, Vol. 14.
"""
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)]