"""Explicit Runge--Kutta tableau definitions used by generic ERK steps.
Each tableau lists the literature source describing the coefficients
so integrators can reference the original derivations.
Published Classes
-----------------
:class:`ERKTableau`
Butcher tableau subclass for explicit Runge--Kutta methods.
>>> tab = ERKTableau(
... a=((0.0, 0.0), (1.0, 0.0)),
... b=(0.5, 0.5), c=(0.0, 1.0), order=2,
... )
>>> tab.stage_count
2
Module-Level Constants
----------------------
:data:`DEFAULT_ERK_TABLEAU`
Alias for :data:`DORMAND_PRINCE_54_TABLEAU`.
:data:`ERK_TABLEAU_REGISTRY`
Mapping from human-readable identifiers (e.g. ``'rk45'``,
``'tsit5'``) to :class:`ERKTableau` instances.
See Also
--------
:class:`~cubie.integrators.algorithms.base_algorithm_step.ButcherTableau`
Parent class providing typed accessors and FSAL detection.
:class:`~cubie.integrators.algorithms.generic_erk.ERKStep`
Step factory that consumes these tableaus.
"""
from typing import Dict
import attrs
from cubie.integrators.algorithms.base_algorithm_step import ButcherTableau
[docs]
@attrs.define
class ERKTableau(ButcherTableau):
"""Coefficient tableau describing an explicit Runge--Kutta scheme.
See Also
--------
:class:`~cubie.integrators.algorithms.base_algorithm_step.ButcherTableau`
Parent class defining ``a``, ``b``, ``b_hat``, ``c``, and
``order`` fields.
"""
#: Heun's improved Euler method offering second-order accuracy.
#: Heun, K. "Neue Methoden zur approximativen Integration der
#: Differentialgleichungen einer unabhängigen Veränderlichen." *Z.
#: Math. Phys.* 45 (1900).
HEUN_21_TABLEAU = ERKTableau(
a=((0.0, 0.0), (1.0, 0.0)),
b=(0.5, 0.5),
c=(0.0, 1.0),
order=2,
)
#: Ralston's third-order, three-stage explicit Runge--Kutta method.
#: Ralston, A. "Runge-Kutta methods with minimum error bounds." *Math. Comp.*
#: 16.80 (1962).
RALSTON_33_TABLEAU = ERKTableau(
a=((0.0, 0.0, 0.0), (0.5, 0.0, 0.0), (0.0, 0.75, 0.0)),
b=(2.0 / 9.0, 1.0 / 3.0, 4.0 / 9.0),
c=(0.0, 1.0 / 2.0, 3.0 / 4.0),
order=3,
)
#: Bogacki--Shampine 3(2) tableau with an embedded error estimate.
#: Bogacki, P. and Shampine, L. F. "An efficient Runge-Kutta (4,5) pair." *J.
#: Comput. Appl. Math.* 46.1 (1993).
BOGACKI_SHAMPINE_32_TABLEAU = ERKTableau(
a=(
(0.0, 0.0, 0.0, 0.0),
(0.5, 0.0, 0.0, 0.0),
(0.0, 0.75, 0.0, 0.0),
(2.0 / 9.0, 1.0 / 3.0, 4.0 / 9.0, 0.0),
),
b=(2.0 / 9.0, 1.0 / 3.0, 4.0 / 9.0, 0.0),
b_hat=(7.0 / 24.0, 1.0 / 4.0, 1.0 / 3.0, 1.0 / 8.0),
c=(0.0, 1.0 / 2.0, 3.0 / 4.0, 1.0),
order=3,
)
#: Dormand--Prince 5(4) tableau with an embedded error estimate.
#: Dormand, J. R. and Prince, P. J. "A family of embedded Runge-Kutta
#: formulae." *Journal of Computational and Applied Mathematics* 6.1 (1980).
DORMAND_PRINCE_54_TABLEAU = ERKTableau(
a=(
(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0),
(1.0 / 5.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0),
(3.0 / 40.0, 9.0 / 40.0, 0.0, 0.0, 0.0, 0.0, 0.0),
(
44.0 / 45.0,
-56.0 / 15.0,
32.0 / 9.0,
0.0,
0.0,
0.0,
0.0,
),
(
19372.0 / 6561.0,
-25360.0 / 2187.0,
64448.0 / 6561.0,
-212.0 / 729.0,
0.0,
0.0,
0.0,
),
(
9017.0 / 3168.0,
-355.0 / 33.0,
46732.0 / 5247.0,
49.0 / 176.0,
-5103.0 / 18656.0,
0.0,
0.0,
),
(
35.0 / 384.0,
0.0,
500.0 / 1113.0,
125.0 / 192.0,
-2187.0 / 6784.0,
11.0 / 84.0,
0.0,
),
),
b=(
35.0 / 384.0,
0.0,
500.0 / 1113.0,
125.0 / 192.0,
-2187.0 / 6784.0,
11.0 / 84.0,
0.0,
),
b_hat=(
5179.0 / 57600.0,
0.0,
7571.0 / 16695.0,
393.0 / 640.0,
-92097.0 / 339200.0,
187.0 / 2100.0,
1.0 / 40.0,
),
c=(
0.0,
1.0 / 5.0,
3.0 / 10.0,
4.0 / 5.0,
8.0 / 9.0,
1.0,
1.0,
),
order=5,
)
#: Classical four-stage Runge--Kutta method introduced by Kutta (1901).
#: Kutta, W. "Beitrag zur näherungsweisen Integration totaler
#: Differentialgleichungen." *Zeitschrift für Mathematik und Physik* 46 (1901).
CLASSICAL_RK4_TABLEAU = ERKTableau(
a=(
(0.0, 0.0, 0.0, 0.0),
(1.0 / 2.0, 0.0, 0.0, 0.0),
(0.0, 1.0 / 2.0, 0.0, 0.0),
(0.0, 0.0, 1.0, 0.0),
),
b=(
1.0 / 6.0,
1.0 / 3.0,
1.0 / 3.0,
1.0 / 6.0,
),
c=(0.0, 1.0 / 2.0, 1.0 / 2.0, 1.0),
order=4,
)
#: Cash--Karp 5(4) tableau with an embedded error estimate.
#: Cash, J. R. and Karp, A. H. "A variable order Runge-Kutta method for initial
#: value problems with rapidly varying right-hand sides." *ACM Transactions on
#: Mathematical Software* 16.3 (1990).
CASH_KARP_54_TABLEAU = ERKTableau(
a=(
(0.0, 0.0, 0.0, 0.0, 0.0, 0.0),
(1.0 / 5.0, 0.0, 0.0, 0.0, 0.0, 0.0),
(3.0 / 40.0, 9.0 / 40.0, 0.0, 0.0, 0.0, 0.0),
(
3.0 / 10.0,
-9.0 / 10.0,
6.0 / 5.0,
0.0,
0.0,
0.0,
),
(
-11.0 / 54.0,
5.0 / 2.0,
-70.0 / 27.0,
35.0 / 27.0,
0.0,
0.0,
),
(
1631.0 / 55296.0,
175.0 / 512.0,
575.0 / 13824.0,
44275.0 / 110592.0,
253.0 / 4096.0,
0.0,
),
),
b=(
37.0 / 378.0,
0.0,
250.0 / 621.0,
125.0 / 594.0,
0.0,
512.0 / 1771.0,
),
b_hat=(
2825.0 / 27648.0,
0.0,
18575.0 / 48384.0,
13525.0 / 55296.0,
277.0 / 14336.0,
1.0 / 4.0,
),
c=(
0.0,
1.0 / 5.0,
3.0 / 10.0,
3.0 / 5.0,
1.0,
7.0 / 8.0,
),
order=5,
)
#: Runge--Kutta--Fehlberg 5(4) tableau with an embedded error estimate.
#: Fehlberg, E. "Low-order classical Runge-Kutta formulas with stepsize control
#: and their application to some heat transfer problems." *NASA Technical
#: Report* 315 (1969).
FEHLBERG_45_TABLEAU = ERKTableau(
a=(
(0.0, 0.0, 0.0, 0.0, 0.0, 0.0),
(1.0 / 4.0, 0.0, 0.0, 0.0, 0.0, 0.0),
(3.0 / 32.0, 9.0 / 32.0, 0.0, 0.0, 0.0, 0.0),
(
1932.0 / 2197.0,
-7200.0 / 2197.0,
7296.0 / 2197.0,
0.0,
0.0,
0.0,
),
(
439.0 / 216.0,
-8.0,
3680.0 / 513.0,
-845.0 / 4104.0,
0.0,
0.0,
),
(
-8.0 / 27.0,
2.0,
-3544.0 / 2565.0,
1859.0 / 4104.0,
-11.0 / 40.0,
0.0,
),
),
b=(
16.0 / 135.0,
0.0,
6656.0 / 12825.0,
28561.0 / 56430.0,
-9.0 / 50.0,
2.0 / 55.0,
),
b_hat=(
25.0 / 216.0,
0.0,
1408.0 / 2565.0,
2197.0 / 4104.0,
-1.0 / 5.0,
0.0,
),
c=(
0.0,
1.0 / 4.0,
3.0 / 8.0,
12.0 / 13.0,
1.0,
1.0 / 2.0,
),
order=5,
)
#: Dormand--Prince 8(5,3) explicit Runge--Kutta method (DOP853).
#: Hairer, Nørsett and Wanner (1993), "Solving Ordinary Differential Equations I", p. 178.
DORMAND_PRINCE_853_TABLEAU = ERKTableau(
a=(
(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0),
(
0.05260015195876773,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
),
(
0.0197250569845379,
0.0591751709536137,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
),
(
0.02958758547680685,
0.0,
0.08876275643042055,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
),
(
0.2413651341592667,
0.0,
-0.8845494793282861,
0.924834003261792,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
),
(
0.037037037037037035,
0.0,
0.0,
0.17082860872947387,
0.12546768756682242,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
),
(
0.037109375,
0.0,
0.0,
0.17025221101954404,
0.06021653898045596,
-0.017578125,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
),
(
0.03709200011850479,
0.0,
0.0,
0.17038392571223999,
0.10726203044637328,
-0.015319437748624402,
0.008273789163814023,
0.0,
0.0,
0.0,
0.0,
0.0,
),
(
0.6241109587160757,
0.0,
0.0,
-3.360892629446941,
-0.868219346841726,
27.592099699446706,
20.154067550477893,
-43.48988418106996,
0.0,
0.0,
0.0,
0.0,
),
(
0.47766253643826435,
0.0,
0.0,
-2.4881146199716676,
-0.590290826836843,
21.230051448181194,
15.279233632882424,
-33.28821096898486,
-0.020331201708508627,
0.0,
0.0,
0.0,
),
(
-0.9371424300859873,
0.0,
0.0,
5.1863724288440636,
1.0914373489967296,
-8.149787010746926,
-18.52006565999696,
22.739487099350504,
2.4936055526796526,
-3.0467644718982197,
0.0,
0.0,
),
(
2.2733101475165383,
0.0,
0.0,
-10.53449546673725,
-2.0008720582248626,
-17.958931863118799,
27.94888452941996,
-2.8589982771350237,
-8.87285693353063,
12.360567175794303,
0.6433927460157635,
0.0,
),
),
b=(
0.05429373411656876,
0.0,
0.0,
0.0,
0.0,
4.450312892752409,
1.8915178993145004,
-5.801203960010585,
0.3111643669578199,
-0.15216094966251608,
0.20136540080403035,
0.04471061572777259,
),
b_hat=(
0.04117368912237387,
0.0,
0.0,
0.0,
0.0,
5.675469339128613,
2.3872768489717506,
-7.465581142465572,
0.6614932157077936,
-0.4863400683755336,
0.20136540080403035,
0.04471061572777259,
),
c=(
0.0,
0.05260015195876773,
0.0789002279381516,
0.1183503419072274,
0.2816496580927726,
0.3333333333333333,
0.25,
0.3076923076923077,
0.6512820512820513,
0.6,
0.8571428571428571,
1.0,
),
order=8,
)
#: Tsitouras 5(4) tableau (Tsit5) with an embedded error estimate.
#: Tsitouras, Ch. (2011). "Runge–Kutta pairs of order 5(4) satisfying only the
#: first column simplifying assumption." Applied Numerical Mathematics, 56(10–11).
# Note: the b_hat vector in the paper cited is not used in most libraries; it
# represents the d=(b - b_hat) vector, and the final entry is incorrect (b7
# = -1/66). I can't find explicit derivations for the fix - it is mentioned
# here https://arxiv.org/pdf/2108.12590#page=1.89 and used in:
# https://github.com/SciML/OrdinaryDiffEq.jl/blob/master/lib/OrdinaryDiffEqTsit5/src/tsit_tableaus.jl
# https://github.com/rtqichen/torchdiffeq/blob/master/torchdiffeq/_impl/tsit5.py
# https://github.com/patrick-kidger/diffrax/blob/main/diffrax/_solver/tsit5.py
tsitouras_b = (
0.09646076681806523,
0.01,
0.4798896504144996,
1.379008574103742,
-3.290069515436081,
2.324710524099774,
0.0,
)
tsitouras_d = (
-1.780011052225771443378550607539534775944678804333659557637450799792588061629796e-03,
-8.164344596567469032236360633546862401862537590159047610940604670770447527463931e-04,
7.880878010261996010314727672526304238628733777103128603258129604952959142646516e-03,
-1.44711007173262907537165147972635116720922712343167677619514233896760819649515e-01,
5.823571654525552250199376106520421794260781239567387797673045438803694038950012e-01,
-4.580821059291869466616365188325542974428047279788398179474684434732070620889539e-01,
1 / 66,
)
tsitouras_b_hat = tuple(b - d for b, d in zip(tsitouras_b, tsitouras_d))
TSITOURAS_54_TABLEAU = ERKTableau(
a=(
(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0),
(0.161, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0),
(-0.008480655492356989, 0.335480655492357, 0.0, 0.0, 0.0, 0.0, 0.0),
(
2.8971530571054935,
-6.359448489975075,
4.3622954328695815,
0.0,
0.0,
0.0,
0.0,
),
(
5.325864828439257,
-11.748883564062828,
7.4955393428898365,
-0.09249506636175525,
0.0,
0.0,
0.0,
),
(
5.86145544294642,
-12.92096931784711,
8.159367898576159,
-0.071584973281401,
-0.028269050394068383,
0.0,
0.0,
),
(
0.09646076681806523,
0.01,
0.4798896504144996,
1.379008574103742,
-3.290069515436081,
2.324710524099774,
0.0,
),
),
b=(
0.09646076681806523,
0.01,
0.4798896504144996,
1.379008574103742,
-3.290069515436081,
2.324710524099774,
0.0,
),
b_hat=tsitouras_b_hat,
c=(0.0, 0.161, 0.327, 0.9, 0.9800255409045097, 1.0, 1.0),
order=5,
)
#: Verner 7(6) explicit Runge--Kutta tableau (Vern7) with embedded error estimate.
#: Constants obtained from Verner's website:
# https://www.sfu.ca/~jverner/RKV76.IIa.Efficient.00001675585.240711.CoeffsOnlyRATandRAD
VERNER_76_TABLEAU = ERKTableau(
a=(
(0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
(1 / 200, 0, 0, 0, 0, 0, 0, 0, 0, 0),
(-4361 / 4050, 2401 / 2025, 0, 0, 0, 0, 0, 0, 0, 0),
(49 / 1200, 0, 49 / 400, 0, 0, 0, 0, 0, 0, 0),
(
2454451729 / 3841600000,
0,
-9433712007 / 3841600000,
4364554539 / 1920800000,
0,
0,
0,
0,
0,
0,
),
(
-6187101755456742839167388910402379177523537620
/ 2324599620333464857202963610201679332423082271,
0,
27569888999279458303270493567994248533230000
/ 2551701010245296220859455115479340650299761,
-37368161901278864592027018689858091583238040000
/ 4473131870960004275166624817435284159975481033,
1392547243220807196190880383038194667840000000
/ 1697219131380493083996999253929006193143549863,
0,
0,
0,
0,
0,
),
(
11272026205260557297236918526339 / 1857697188743815510261537500000,
0,
-48265918242888069 / 1953194276993750,
26726983360888651136155661781228 / 1308381343805114800955157615625,
-2090453318815827627666994432 / 1096684189897834170412307919,
1148577938985388929671582486744843844943428041509
/ 1141532118233823914568777901158338927629837500000,
0,
0,
0,
0,
),
(
1304457204588839386329181466225966641
/ 108211771565488329642169667802016000,
0,
-1990261989751005 / 40001418792832,
2392691599894847687194643439066780106875
/ 58155654089143548047476915856270826016,
-1870932273351008733802814881998561250
/ 419326053051486744762255151208232123,
1043329047173803328972823866240311074041739158858792987034783181
/ 510851127745017966999893975119259285040213723744255237522144000,
-311918858557595100410788125 / 3171569057622789618800376448,
0,
0,
0,
),
(
17579784273699839132265404100877911157
/ 1734023495717116205617154737841023480,
0,
-18539365951217471064750 / 434776548575709731377,
447448655912568142291911830292656995992000
/ 12511202807447096607487664209063950964109,
-65907597316483030274308429593905808000000
/ 15158061430635748897861852383197382130691,
273847823027445129865693702689010278588244606493753883568739168819449761
/ 136252034448398939768371761610231099586032870552034688235302796640584360,
694664732797172504668206847646718750
/ 1991875650119463976442052358853258111,
-19705319055289176355560129234220800
/ 72595753317320295604316217197876507,
0,
0,
),
(
-511858190895337044664743508805671
/ 11367030248263048398341724647960,
0,
2822037469238841750 / 15064746656776439,
-23523744880286194122061074624512868000
/ 152723005449262599342117017051789699,
10685036369693854448650967542704000000
/ 575558095977344459903303055137999707,
-6259648732772142303029374363607629515525848829303541906422993
/ 876479353814142962817551241844706205620792843316435566420120,
17380896627486168667542032602031250
/ 13279937889697320236613879977356033,
0,
0,
0,
),
),
c=(
0,
1 / 200,
49 / 450,
49 / 300,
911 / 2000,
3480084980 / 5709648941,
221 / 250,
37 / 40,
1,
1,
),
b=(
96762636172307789 / 2051985304794103980,
0,
0,
312188947591288252500000 / 1212357694274963646019729,
13550580884964304000000000000 / 51686919683339547115937980629,
72367769693133178898676076432831566019684378142853445230956642801
/ 475600216991873963561768100160364792981629064220601844848928537580,
1619421054120605468750 / 3278200730370057108183,
-66898316144057728000 / 227310933007074849597,
181081444637946577 / 2226845467039736466,
0,
),
b_hat=(
117807213929927 / 2640907728177740,
0,
0,
4758744518816629500000 / 17812069906509312711137,
1730775233574080000000000 / 7863520414322158392809673,
2682653613028767167314032381891560552585218935572349997
/ 12258338284789875762081637252125169126464880985167722660,
40977117022675781250 / 178949401077111131341,
0,
0,
2152106665253777 / 106040260335225546,
),
order=7,
)
DEFAULT_ERK_TABLEAU = DORMAND_PRINCE_54_TABLEAU
"""Default tableau used when constructing generic ERK integrators."""
ERK_TABLEAU_REGISTRY: Dict[str, ERKTableau] = {
"heun-21": HEUN_21_TABLEAU,
"ralston-33": RALSTON_33_TABLEAU,
"bogacki-shampine-32": BOGACKI_SHAMPINE_32_TABLEAU,
"rk23": BOGACKI_SHAMPINE_32_TABLEAU,
"ode23": BOGACKI_SHAMPINE_32_TABLEAU,
"dormand-prince-54": DORMAND_PRINCE_54_TABLEAU,
"dopri54": DORMAND_PRINCE_54_TABLEAU,
"rk45": DORMAND_PRINCE_54_TABLEAU,
"ode45": DORMAND_PRINCE_54_TABLEAU,
"classical-rk4": CLASSICAL_RK4_TABLEAU,
"rk4": CLASSICAL_RK4_TABLEAU,
"cash-karp-54": CASH_KARP_54_TABLEAU,
"fehlberg-45": FEHLBERG_45_TABLEAU,
"dormand-prince-853": DORMAND_PRINCE_853_TABLEAU,
"dop853": DORMAND_PRINCE_853_TABLEAU,
"tsit5": TSITOURAS_54_TABLEAU,
"Tsit5": TSITOURAS_54_TABLEAU,
"vern7": VERNER_76_TABLEAU,
"Vern7": VERNER_76_TABLEAU,
}
"""Mapping from human readable identifiers to ERK tableaus."""