########################################################################################
##
## EMBEDDED DIAGONALLY IMPLICIT RUNGE KUTTA METHOD
## (solvers/esdirk32.py)
##
## Milan Rother 2024
##
########################################################################################
# IMPORTS ==============================================================================
import numpy as np
from ._rungekutta import DiagonallyImplicitRungeKutta
# SOLVERS ==============================================================================
[docs]
class ESDIRK43(DiagonallyImplicitRungeKutta):
"""Six-stage, 4th order Embedded Singly Diagonally Implicit Runge-Kutta (ESDIRK) method.
Features an embedded 3rd order method for adaptive step size control. The first stage is
explicit. L-stable and stiffly accurate, making it excellent for stiff problems requiring
moderate-to-high accuracy with adaptive timestepping.
Characteristics
---------------
* Order: 4
* Embedded Order: 3
* Stages: 6 (1 Explicit, 5 Implicit)
* Implicit (ESDIRK)
* Adaptive timestep
* L-stable
* Stiffly accurate
When to Use
-----------
* **Stiff problems with adaptive stepping**: Excellent default for stiff ODEs
* **Moderate-to-high accuracy**: 4th order with good error control
* **Damping high frequencies**: L-stability damps spurious oscillations
* **General-purpose stiff solver**: Reliable choice for most stiff applications
**Recommended** as a default adaptive stiff solver. For very high accuracy, use ESDIRK54.
For non-stiff problems, RKDP54 is more efficient.
References
----------
.. [1] Kennedy, C. A., & Carpenter, M. H. (2019). "Diagonally implicit Runge-Kutta
methods for stiff ODEs". Applied Numerical Mathematics, 146, 221-244.
.. [2] 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 = 6
#order of scheme and embedded method
self.n = 4
self.m = 3
#flag adaptive timestep solver
self.is_adaptive = True
#intermediate evaluation times
self.eval_stages = [0.0, 1/2, (2-np.sqrt(2))/4, 2012122486997/3467029789466, 1.0, 1.0]
#butcher table
self.BT = {
0: None, # explicit first stage
1: [1/4, 1/4],
2: [-1356991263433/26208533697614, -1356991263433/26208533697614, 1/4],
3: [-1778551891173/14697912885533, -1778551891173/14697912885533,
7325038566068/12797657924939, 1/4],
4: [-24076725932807/39344244018142, -24076725932807/39344244018142,
9344023789330/6876721947151, 11302510524611/18374767399840, 1/4],
5: [657241292721/9909463049845, 657241292721/9909463049845,
1290772910128/5804808736437, 1103522341516/2197678446715, -3/28, 1/4]
}
#coefficients for truncation error estimate
_A1 = [
657241292721/9909463049845, 657241292721/9909463049845,
1290772910128/5804808736437, 1103522341516/2197678446715, -3/28, 1/4
]
_A2 = [
-71925161075/3900939759889, -71925161075/3900939759889,
2973346383745/8160025745289, 3972464885073/7694851252693,
-263368882881/4213126269514, 3295468053953/15064441987965
]
self.TR = [a1 - a2 for a1, a2 in zip(_A1, _A2)]