on Q-D recursive Runge-Kutta methods

Overview

He and Huang construct explicit Runge—Kutta methods of arbitrary even order by recursively controlling Q- and D-type residuals (He & Huang, 2026). Their construction should be read against Stepanov’s Q/D formulation for order-ten explicit methods (Stepanov, 2025): Stepanov shows how branch-side and root-side order-condition defects can be made redundant through stage-order layers and repeated-node clusters; He and Huang turn that principle into a uniform algorithm based on Lobatto nodes, coordinate support separation, and two structured linear systems.

An $s$-stage explicit Runge—Kutta method is a triple $(A,b,c)$, where $A=(a_{ij})$ is strictly lower triangular, $b=(b_i)$ is the weight vector, and $c=(c_i)$ is the node vector. Applied to $y^{\prime }=f(t,y)$, it computes

$$Y_i=y_n+h\sum _{j<i}{a}_{ij}f(t_n+c_{j}h,Y_j),i=1,\dots ,s,$$

and

$$y_{n+1}=y_n+h\sum _{i=1}^s{b}_{i}f(t_n+c_{i}h,Y_i).$$

The convention throughout is

$$A1=c,$$

where $1=(1,\dots ,1{)}^T$. This is the standard consistency condition used in the paper.

For each even $p\ge 4$, He and Huang’s construction produces a method with

$$s(p)=\frac{p^2-2p+8}{4}$$

stages. These counts are not asserted to be minimal. The point is the recursive certificate: instead of solving the nonlinear rooted-tree order equations directly, the construction places every residual generated by the order-condition recursion into a subspace that is later annihilated.

🏷️ Definitions and Order Conditions

Runge—Kutta order conditions are naturally indexed by rooted trees. If $\tau$ is a rooted tree, write $|\tau |$ for its number of vertices and $\tau !$ for its tree factorial. Define the stage vector $\Phi (\tau )$ by

$$\Phi (\bullet )=1,$$

and, if $\tau =[{\tau }_1,\dots ,{\tau }_r]$,

$$\Phi (\tau )=\prod _{\alpha =1}^{r}A\Phi ({\tau }_{\alpha }),$$

with componentwise product. The method has order at least $p$ exactly when

$$b^T\Phi (\tau )=\frac{1}{\tau !}\text{for every rooted tree }|\tau |\le p.$$

The first few conditions are

$$b^{T}1=1,b^{T}c=\frac{1}{2},b^T{c}^2=\frac{1}{3},b^{T}Ac=\frac{1}{6}.$$

The first three are quadrature moment conditions. The last one is already an internal-stage condition: it asks whether information propagated through $A$ matches the Taylor expansion.

Butcher’s simplifying assumptions organize these conditions as

$$B(q):b^T{c}^k=\frac{1}{k+1},0\le k\le q-1,$$ $$C(q):A{c}^k=\frac{c^{k+1}}{k+1},0\le k\le q-1,$$

and

$$D(q):(b\odot c^k{)}^{T}A=\frac{b\odot (1-c^{k+1})}{k+1},0\le k\le q-1.$$

Here powers and $\odot$ are componentwise. The condition $B(q)$ is quadrature. The condition $C(q)$ is internal stage order. The condition $D(q)$ is the left-sided companion condition that appears when the root of a tree is separated from its branches.

For explicit methods, high $C(q)$ cannot hold literally. Since the second row of a strictly lower triangular $A$ gives $(Ac{)}_2=0$, the condition $C(2)$ would imply

$$0=(Ac{)}_2=\frac{c_2^2}{2},$$

so $c_2=0$. Similar constraints propagate through early stages and conflict with nontrivial quadrature. The obstruction is not that $C$ and $D$ are irrelevant; it is that explicit methods can only satisfy weakened, residual versions of them.

🏷️ Q/D Residual Calculus

Define the $C$-residuals

$$q_k=A{c}^k-\frac{c^{k+1}}{k+1},k\ge 0,$$

and the $D$-residuals

$$d_k=(b\odot c^k{)}^{T}A-\frac{b\odot (1-c^{k+1})}{k+1},k\ge 0.$$

Thus $q_k=0$ is the $k$th $C$ equation and $d_k=0$ is the $k$th $D$ equation. The central idea is to allow $q_k$ and $d_k$ to be nonzero, but only in directions that are killed by later contractions.

The order-three chain condition illustrates the mechanism. Since

$$Ac=\frac{c^2}{2}+q_1,$$

one has

$$b^{T}Ac=\frac{1}{2}{b}^T{c}^2+b^T{q}_1.$$

If $B(3)$ holds, the first term is $1/6$. The chain condition is therefore true provided $b^T{q}_1=0$. The residual $q_1$ does not need to vanish; it only needs to be invisible to $b$.

At higher order, residuals are repeatedly acted on by $A$, multiplied by powers of $c$, and multiplied componentwise with other residuals. This is why the useful object is not an individual residual but a recursively closed residual space. With ${\Phi }_1=\mathrm{span}\{1,c\}$, He and Huang use

$$Q_0=\{0\},Q_1=\mathrm{span}\{q_0\},$$ $$Q_2=\mathrm{span}\{q_1,A{q}_0,q_0\odot c,q_0\},$$

and, for $k\ge 3$,

$$Q_k=\mathrm{span}\{q_{k-1},A{Q}_{k-1},Q_{k-1}\odot {\Phi }_1\}.$$

The corresponding row-vector spaces are

$$D_0=\{0\},D_1=\mathrm{span}\{d_0\},$$ $$D_2=\mathrm{span}\{d_1,d_{0}A,d_0\odot c,d_0\},$$

and, for $k\ge 3$,

$$D_k=\mathrm{span}\{d_{k-1},D_{k-1}A,D_{k-1}\odot {\Phi }_1\}.$$

These definitions are closure statements. Every operation that appears in the rooted-tree recursion sends an allowed error to a next-level allowed error. This is why the proof can control all trees without enumerating their equations one by one.

🏷️ Stepanov’s Preceding Idea

Stepanov defines Q-type defects for arbitrary rooted trees by

$$Q(t)=A\Phi (t)-\frac{c^{|t|}}{t!},$$

so the usual subquadrature vectors are the height-one cases $Q([{\bullet }^n])=q_n$. He also defines D-type row defects

$$D(t)=(b\odot \Phi (t){)}^{T}A-\frac{b\odot (1-c^{|t|})}{t!},$$

whose height-one cases are $d_n$. In that language, the full order conditions split into quadrature conditions, Q-type conditions, and D-type conditions. Q-type conditions may contain several branch defects, while D-type conditions contain one root-side defect. This branch/root distinction is precisely the distinction that He and Huang encode in the spaces $Q_k$ and $D_k$.

Stepanov’s order-ten construction uses two organizing heuristics. Stage Order Layers make selected stages have enough stage order that many Q-type conditions become redundant. Counterpoised Node Clusters group stages with identical nodes and place mutually orthogonal Q- and D-subspaces on each cluster; this is what lets D-type moment conditions vanish by cancellation inside a repeated-node cluster. This leads to a seven-parameter family of explicit order-ten methods with 15 stages.

He and Huang preserve the cancellation logic but make it systematic. Instead of designing a special order-ten architecture, they prescribe recursive spaces, support regions, and repeated-node clusters for every even order. The price is that their order-ten member uses 22 stages rather than Stepanov’s 15. The gain is a uniform construction with a proof that reduces coefficient generation to two structured linear systems.

🏷️ Sufficient Q/D Conditions

Let $m,n$ satisfy

$$m\ge n-1,m+n+1\ge p.$$

The paper proves that order $p$ follows from

$$B(p):b^T{c}^{k-1}=\frac{1}{k},k=1,\dots ,p,$$ $$QO(m):b\odot Q_m=\{0\},$$ $$DO(n):D_n{c}^{k-1}=\{0\},k=1,\dots ,p-n,$$ $$QD(m,n):Q_m\odot D_n=\{0\},$$

and

$$QR(m):Q_{m_1}\odot Q_{m_2}\subseteq Q_{m_1}(m_2\le m_1\le m).$$

This is a routing theorem. Pure powers of $c$ are handled by $B(p)$. Branch-level residuals are routed into the $Q$ hierarchy and killed by $QO(m)$. Root-level residuals are routed into the $D$ hierarchy and killed by $DO(n)$. Mixed branch/root terms are killed by $QD(m,n)$. Products of Q-residuals remain controlled by $QR(m)$.

The recursion works because it replaces a large nonlinear family of rooted-tree equations by an invariant-subspace statement: every error term produced by the tree recursion stays inside a space that the construction knows how to annihilate.

🏷️ Recursive Tableau Construction

For even $p\ge 4$, set

$$m=\frac{p}{2}-1,n=\frac{p}{2},\ell =\frac{m(m-1)}{2}.$$

Then

$$s=1+\ell +\frac{n(n+1)}{2}=\frac{p^2-2p+8}{4}.$$

The node vector $c$ is built from the $(n+1)$-point Gauss—Lobatto rule on $[0,1]$. Since this rule is exact through degree $2n-1=p-1$, the moment conditions $B(p)$ are built in. When a Lobatto node is repeated across several stages, the corresponding quadrature weight is split across that repeated-node cluster.

The stage indices are divided into coordinate regions. The weights $b$ live on the first region, the $Q_m$ residuals are forced into the middle region, and the $D_n$ residuals are forced into the last region. Hence

$$b\odot Q_m=0,Q_m\odot D_n=0$$

by support separation.

The entries of $A$ are then determined as follows.

1. Choose the Lobatto nodes and initialize $b,c$.
2. Initialize $A=0$.
3. Solve the D-system in the bottom-right part of $A$, enforcing the recursive inclusions for $D_1,\dots ,D_n$.
4. Insert those coefficients into $A$.
5. Solve the Q-system in the left part of $A$, enforcing the recursive inclusions for $Q_1,\dots ,Q_m$.
6. Return the explicit tableau $(A,b,c)$.

The D-system is solved first because the Q-equations use the bottom-right coefficients already fixed by the D-equations. Once the support pattern is chosen, both systems are linear. This is the constructive content of the paper: the nonlinear rooted-tree order problem is converted into two structured linear solves.

The stage counts for the first few even orders are:

order $p$	stages $s(p)$
$4$	$4$
$6$	$8$
$8$	$14$
$10$	$22$
$12$	$32$
$14$	$44$

These counts are upper-bound construction counts, not minimality claims.

🏷️ The Order-Six Mechanism

For $p=6$,

$$m=2,n=3,\ell =1,s=8.$$

The four Lobatto nodes on $[0,1]$ are

$$0,x_{-}=\frac{5-\sqrt{5}}{10},x_{+}=\frac{5+\sqrt{5}}{10},1.$$

The construction arranges

$$c=(0,x_{-},x_{+},x_{-},x_{+},x_{-},x_{+},1)$$

and

$$b=\left(\frac{1}{12},0,\frac{5}{36},\frac{5}{24},\frac{5}{36},\frac{5}{24},\frac{5}{36},\frac{1}{12}\right).$$

The zero weight at stage $2$ is intentional: this stage belongs to the Q-block, so it can carry a branch residual without being seen by $b$.

For the generated order-six tableau,

$$q_0=A1-c=0,$$

but, rounded to six significant digits,

$$q_1\approx (0,-0.038197,0,0,0,0,0,0{)}^T.$$

Thus $C(2)$ fails, but it fails in a harmless coordinate. Since $b_2=0$, this residual is killed by the final weight contraction.

Similarly,

$$d_0=0,$$

while

$$d_1\approx (0,0,-0.005305,0,0,0,0.005305,0).$$

The two nonzero entries sit at repeated copies of the same Lobatto node and cancel in the relevant cluster moments. This is the D-side analogue of the Q-block cancellation.

📊 Numerical Verification

The authoritative verification script is

content/codes/2026 Summer/qd_rk_bigfloat_verification.jl

It generates the tableaux by the recursive algorithm, using Julia BigFloat arithmetic throughout coefficient generation and rooted-tree checking. The saved coefficient file is decimal text, not a binary double array:

content/codes/2026 Summer/qd_rk_tableaux_p4_p14_bigfloat.txt

The script performs the following steps for each even order $p=4,6,\dots ,14$:

1. construct the $(p/2+1)$-point Lobatto rule on $[0,1]$ in high precision;
2. assemble the Q-block nodes, D ghost nodes, and D triangular repeated-node clusters;
3. split each Lobatto weight across the stages carrying that node;
4. solve the D-system in BigFloat arithmetic;
5. solve the Q-system in BigFloat arithmetic using the D coefficients just computed;
6. enumerate all unlabeled rooted trees through order $p$;
7. verify

$$b^T\Phi (\tau )=\frac{1}{\tau !}$$

for every rooted tree $|\tau |\le p$.

The high-precision run is

julia 'content/codes/2026 Summer/qd_rk_bigfloat_verification.jl' --orders 4 6 8 10 12 14 --dps 80 --save 'content/codes/2026 Summer/qd_rk_tableaux_p4_p14_bigfloat.txt'

The completed BigFloat output confirms all even orders $4\le p\le 14$:

order $p$	stages	rooted trees checked	max tree residual	max $B(p)$ residual	row-sum residual
$4$	$4$	$8$	$1.23\cdot 10^{-91}$	$0$	$0$
$6$	$8$	$37$	$3.93\cdot 10^{-90}$	$3.93\cdot 10^{-90}$	$0$
$8$	$14$	$200$	$9.82\cdot 10^{-91}$	$9.82\cdot 10^{-91}$	$3.93\cdot 10^{-90}$
$10$	$22$	$1205$	$2.95\cdot 10^{-90}$	$2.95\cdot 10^{-90}$	$7.85\cdot 10^{-90}$
$12$	$32$	$7813$	$3.93\cdot 10^{-90}$	$3.93\cdot 10^{-90}$	$7.85\cdot 10^{-90}$
$14$	$44$	$53272$	$3.93\cdot 10^{-90}$	$3.93\cdot 10^{-90}$	$3.14\cdot 10^{-89}$

The generated systems also close at high precision:

order $p$	D unknowns	D solve residual	Q unknowns	Q solve residual
$4$	$3$	$0$	$3$	$0$
$6$	$13$	$2.45\cdot 10^{-91}$	$13$	$0$
$8$	$38$	$2.45\cdot 10^{-91}$	$45$	$4.91\cdot 10^{-91}$
$10$	$89$	$3.68\cdot 10^{-91}$	$122$	$1.96\cdot 10^{-90}$
$12$	$180$	$9.82\cdot 10^{-91}$	$276$	$1.96\cdot 10^{-90}$
$14$	$328$	$5.52\cdot 10^{-91}$	$548$	$1.96\cdot 10^{-90}$

As an independent sanity check, the same BigFloat tableaux were applied to

$$y^{\prime }=ty,y(0)=1,0\le t\le 1,$$

and compared with the exact value $e^{1/2}$. This nonautonomous test uses the stage nodes $c_i$ explicitly. The log-log plot shows the expected slopes $4,6,8,10,12,14$ under step refinement.

These ODE runs do not replace the rooted-tree verification; they check the expected convergence behavior on a problem with a known exact solution. The rooted-tree computation is the stronger nonlinear order check.

🏷️ Order-Ten Assembly

For $p=10$,

$$m=4,n=5,\ell =6,s=22.$$

The construction uses a six-point Lobatto rule. The D-system has $89$ unknowns, split by levels as

$$14,13,22,24,16,$$

and the Q-system has $122$ unknowns, split as

$$21,20,36,45.$$

Stepanov’s order-ten construction uses 15 stages, so it is substantially tighter at this order. He and Huang’s order-ten member should therefore not be read as a best-known order-ten method. Its role is different: it is the $p=10$ instance of a uniform even-order recursive family.

🏷️ Interpretation and Limitations

The construction separates three tasks. Lobatto quadrature supplies the moment conditions $B(p)$. Support placement makes Q-residuals invisible to $b$ and separates them from D-residuals. Repeated-node clusters make D-residual moments cancel. The coefficient problem then splits into a D-system followed by a Q-system.

The Q/D conditions are sufficient, not necessary. They do not characterize all explicit Runge—Kutta methods of a given order, and the stage count is an upper-bound construction rather than a minimality theorem. Free parameters can still influence stability regions and error constants; algebraic order alone does not determine practical performance.

Links

• on From ODE45 to others gives the practical MATLAB-side viewpoint on explicit Runge—Kutta solvers and when nonstiff one-step methods are a reasonable default.
• on Wave Equation with PML includes explicit Runge—Kutta—Chebyshev time stepping, where stability rather than only algebraic order controls performance.
• on epsilon extrapolation connects to extrapolation methods, which have the same leading $p^2/4$ stage-growth scale.
• on Crouzeix’s conjecture and the 1+sqrt(2) bound is useful background for thinking about stability through polynomial bounds on matrices.

References

🐻  He, J. & Huang, J. 2026. Low Stage High Order Explicit Runge–Kutta Methods via Q- and D-Conditions: General Theory and Efficient Recursive Construction.

🐻  Stepanov, M. 2025. On Runge–Kutta Methods of Order 10. Electronic Transactions on Numerical Analysis 63, 609–626.