Runge-Kutta tables

Question: do the coefficients of a Runge-Kutta method have the order and the stability properties the method claims?

This is pure algebra on the Butcher table (A, b, c): no test problem, no step sizes, no tolerance games. It is the check to run on a library of methods, for example in continuous integration.

"""Check Runge-Kutta tables by algebra: order, linear order, stage order, stability."""
import math

from physverify.butcher import check_tableau, r_infinity, stage_order

RK4 = ([[0, 0, 0, 0], [.5, 0, 0, 0], [0, .5, 0, 0], [0, 0, 1, 0]], [1 / 6, 1 / 3, 1 / 3, 1 / 6], [0, .5, .5, 1])
print("RK4:", check_tableau(*RK4, claimed=4).message)

# labelled order 3, but only order 2; its linear order is 3, so a test on y' = -y looks fine
A, b, c = [[0, 0, 0], [1 / 3, 0, 0], [0, 1, 0]], [1 / 2, 0, 1 / 2], [0, 1 / 3, 1]
r = check_tableau(A, b, c, claimed=3)
print("labelled 3:", r.status, "-", r.message)
print("labelled 2, linear 3:", check_tableau(A, b, c, claimed=2, claimed_linear=3).status)

# a typo in the nodes c: the order conditions cannot see it, the row sums can
r = check_tableau(*RK4[:2], [0, .5, .4, 1], claimed=4)
print("typo in c:", r.status, "-", r.message.split(";")[0])

g = 1 - 1 / math.sqrt(2)  # two-stage SDIRK: L-stable and stiffly accurate
SDIRK = ([[g, 0], [1 - g, g]], [1 - g, g])
print("SDIRK:", check_tableau(*SDIRK, claimed=2, l_stable=True, stiffly_accurate=True).message)

TRAP = ([[0, 0], [.5, .5]], [.5, .5])  # trapezoid rule: A-stable but not L-stable
r = check_tableau(*TRAP, claimed=2, l_stable=True)
print("trapezoid claimed L-stable:", r.status, "-", r.message.split(";")[0])
print(f"R(inf) of the trapezoid rule: {r_infinity(*TRAP):.3g} | stage order of SDIRK:", stage_order(SDIRK[0], b=SDIRK[1]))

Output:

RK4: order 4 confirmed by the order conditions (1: 1.1e-16, 2: 0.0e+00, 3: 0.0e+00, 4: 0.0e+00, 5: 1.2e-02)
labelled 3: FAIL - claims order 3 but satisfies the order conditions only up to 2 (largest residual per order: 1: 0.0e+00, 2: 0.0e+00, 3: 1.7e-01, 4: 2.5e-01); linear order 3 is higher than the order 2: a convergence test on y' = lambda*y alone would overstate the order
labelled 2, linear 3: PASS
typo in c: FAIL - row sums of A differ from c by 1.0e-01
SDIRK: order 2 confirmed by the order conditions (1: 0.0e+00, 2: 1.1e-16, 3: 4.0e-02); A-stable as stated; L-stable as stated; stiffly accurate as stated
trapezoid claimed L-stable: FAIL - stated L-stable but it is not (R(inf) = -1.0e+00)
R(inf) of the trapezoid rule: -1 | stage order of SDIRK: 1

What each part checks

Argument	Property	How
claimed	Order: the error is O(h^(p+1)) per step on any smooth ODE	Butcher's order conditions, one per rooted tree with up to p vertices
claimed_linear	Linear order: the order on y' = lambda*y only	The stability function R(z) matches e^z up to z^p
claimed_stage	Stage order: how accurate the internal stages are	Simplifying assumptions C(q) and B(q); a low stage order causes order loss on stiff problems
a_stable	Stable for every decaying linear mode	|R(z)| <= 1 on the imaginary axis (an exact polynomial test) and no poles of R in the left half-plane
l_stable	A-stable, and very stiff modes are damped completely	A-stable and R(infinity) = 0
stiffly_accurate	The step result is the last stage	b equals the last row of A
c given	The nodes are consistent	Each row sum of A equals c

Each property is checked only when you state it, and a mismatch fails in either direction: claiming L-stability for the trapezoid rule fails, and so does claiming that an L-stable method is not.

Why the linear order matters

A convergence test on y' = -y (the "Dahlquist" test equation) measures the linear order, which can be higher than the real order. The table above that claims order 3 is order 2 on general problems but order 3 on y' = -y, so a library whose only test is y' = -y would never notice. A public Runge-Kutta library had exactly this case; physverify found it by the order conditions, and the library now records both orders.

Rounded coefficients

The default tolerance (1e-8) accepts coefficients typed with 9 or 10 digits. When a missing order is only a matter of too few digits, the message says so and suggests a tolerance to retry with.

Lower-level functions

tableau_order, linear_order, stage_order, stability_function, r_infinity, a_stability_margin and is_stiffly_accurate in physverify.butcher return the numbers themselves.

Not yet covered

The extra conditions that link the explicit and implicit tables of an IMEX (additive) pair; see the Roadmap.