Stiffly Accurate Radau Quadrature

• Stiffly accurate Radau quadrature is a high-order implicit Runge–Kutta integration scheme that uses right-Radau collocation points to ensure the final stage equals the computed solution.
• It maintains L-stability and prevents order reduction in stiff ODEs and PDEs by aligning the final stage with the overall step update.
• The method supports adaptivity, high-precision arithmetic, and efficient linear algebra to optimize performance over traditional solvers like BDF.

Stiffly accurate Radau quadrature refers to a class of numerical integration schemes for stiff ODEs and PDEs, specifically the implicit Runge–Kutta methods known as Radau IIA. These methods employ collocation at so-called right-Radau quadrature points, producing algebraically stable, high-order integration formulas. The key attribute, stiff accuracy, ensures that the final stage of the method coincides with the computed solution after a step, preventing order reduction and maintaining stability even in the presence of strongly stiff dynamics. Stiffly accurate Radau methods constitute the state-of-the-art for high-fidelity, robust time integration in stiff ODE/PDE contexts, offering adaptivity in both order and precision, and have significant theoretical and practical advantages over alternatives such as BDF and Rosenbrock–Wanner schemes (Ekanathan et al., 2024, Duarte et al., 2014).

1. Radau IIA Collocation Methods and Quadrature Formulas

Radau IIA methods are constructed by collocating the ODE solution at $s$ nodes on the unit interval $[0,1]$, corresponding to the $s-1$ roots in $(0,1)$ of the polynomial

$$\frac{d^{s-1}}{dx^{s-1}} \left[x^{s-1}(x-1)^s\right] = 0,$$

along with $c_s = 1$. These points are the right-Radau nodes, as described by the roots of the derivative of degree $s$ Legendre polynomials mapped onto $[0,1]$. The unique feature is the inclusion of the endpoint $t = 1$ in the collocation nodes (Duarte et al., 2014).

For a general initial value problem $y' = f(t, y)$, $y(t_n) = y_n$ with step $h$, the standard Butcher form is

$$\begin{aligned} k_i &= f\!\left(t_n + c_i h,\;y_n + h\sum_{j=1}^s a_{ij} k_j\right), \quad i=1, \dots, s, \ y_{n+1} &= y_n + h \sum_{i=1}^s b_i k_i. \end{aligned}$$

The $s \times s$ coefficient matrix $A = (a_{ij})$, the weights $b = (b_i)$, and the nodes $c = (c_i)$ together form the Butcher tableau for the method.

By construction, the collocation conditions enforce

$$\sum_{j=1}^s a_{ij} c_j^{q-1} = \frac{c_i^q}{q}, \qquad q = 1, ..., s,$$

for each $i$, and the weights are determined as $b_j = a_{s, j}$. This guarantees order $2s-1$ and stage order $s$ (Ekanathan et al., 2024, Duarte et al., 2014).

2. Stiff Accuracy: Definition and Theoretical Implications

A Runge–Kutta method is stiffly accurate (or B-stiffly accurate) if the weights $b_j$ equal the last row of the $A$ matrix:

$a_{s, j} = b_j, \; \forall j = 1, ..., s.$

This ensures that the final stage,

$$k_s = f\left(t_n + h, y_n + h \sum_{j=1}^s a_{s, j} k_j\right) = f\left(t_n + h, y_{n+1}\right),$$

uses the value $y_{n+1}$, ensuring $u_{n+1} = U_s$, where $U_s$ is the last stage. This property is crucial in the context of stiff problems because it suppresses spurious high-frequency components and avoids the order-reduction phenomenon often observed in the integration of stiff systems. B-stability proofs for linear stiff test equations rely on this identity, and it is critical for enforcing $L$-stability, as the amplification factor $R(z)$ satisfies $R(z)\to 0$ as $z\to\infty$ (Ekanathan et al., 2024, Duarte et al., 2014).

3. Computational Construction of Radau IIA Coefficients

The Radau IIA coefficients can be computed to arbitrary order and arithmetic precision:

• The $s$ nodes $c_i$ are determined by finding the roots of the defining polynomial or by leveraging Legendre polynomial derivatives mapped to $[0,1]$.
• The $s \times s$ Vandermonde-type linear system from the collocation/order conditions is solved for each row to obtain the $a_{ij}$, with weights $b_j = a_{s, j}$ (Ekanathan et al., 2024).
• Numerical implementations may rely on companion matrix eigensolvers for node computation and cache the resulting coefficients~(parameterized by order and arithmetic) to optimize repeated evaluations.
• All calculations (coefficient derivation, time stepping) can be performed in arbitrary precision (e.g., BigFloat), and mixed-precision strategies are supported (Ekanathan et al., 2024).

4. Embedded Tableaux, Local Error Estimation, and Adaptivity

For adaptive step-size and order control, an embedded Radau formula of lower order is derived by enforcing the first $s$ moment conditions,

$$\sum_{j=1}^s \hat{b}_j c_j^{k-1} = \frac{1}{k},\quad k = 1, \ldots, s,$$

for new weights $\hat{b}_j$. This embedded method shares the $A$ and $c$ matrices but is only order $s$. The local error estimate is computed as

$$\mathrm{err} = \left\| h \sum_{j=1}^s (b_j - \hat{b}_j) k_j \right\|.$$

A PI-type or predictive controller adjusts the next step size using this estimate. Order-adaptive strategies track smoothed Newton iteration metrics and adjust order up or down (by 4 stages at a time) depending on mean iteration count, allowing both aggressive order escalation when stiffness is low and stability when iteration cost rises (Ekanathan et al., 2024).

5. Algorithmic Modernizations and Linear Algebra Optimization

Contemporary implementations integrate key modernizations:

• High-precision arithmetic: All stages in tableau construction and time-stepping can use arbitrary (and mixed) floating-point precisions, which is essential for high-fidelity or near-machine-precision integrations.
• Real block-diagonalization of $A^{-1}$: $A^{-1}$ in the Radau IIA method has one real eigenvalue and $(s-1)/2$ pairs of complex-conjugate eigenvalues for odd $s$. By block-diagonalizing $A^{-1}$ into $1 \times 1$ and $2 \times 2$ real blocks, Newton iteration for implicit solves can use specialized (and cheaper) factorizations.
• Step-size and order adaptivity: Combined step-size and order controllers respond to local error and iteration cost metrics, maintaining optimal efficiency even across wide changes in problem stiffness (Ekanathan et al., 2024).

6. Error Bounds, Convergence, and Stability Analysis

Radau IIA methods provide stage-order $q = s$ and global order $p = 2s-1$. The local quadrature error at each stage scales as $O(h^{s+1})$, while the full step error is $O(h^{2s})$. The stiff accuracy property assures that the last-stage update preserves order in the stiff regime, preventing classical order reduction due to fast transients. $A$- and $L$-stability ensure unconditional stability for linear test problems, and $L$-stability further ensures that the numerical solution damps to zero for arbitrarily stiff (large negative real part) eigenvalues. These properties have been rigorously quantified in the deferred-correction–splitting context for coupled ODE/PDE systems (Duarte et al., 2014).

7. Empirical Performance: Evaluations and Benchmarks

Recent benchmarking on canonical stiff ODE test problems (Oregonator, Robertson, Hires, Pollution) demonstrates that fully adaptive Radau IIA implementations achieve:

• Approximately $2\times$ speedup over classic Fortran Radau IIA (orders~5, 9, 13) over tolerances $10^{-5}$ to $10^{-14}$.
• Outperformance versus high-order Rosenbrock–Wanner (e.g. Rodas5P) and BDF integrators at low error tolerances ($10^{-10}$ and lower), including in extended-precision contexts (BigFloat) (Ekanathan et al., 2024).
• The stiff accuracy property is essential in these regimes, as it avoids loss of global order and enables robust error control, permitting large $h$ where appropriate without losing stability.

---

Key references:

• "A Fully Adaptive Radau Method for the Efficient Solution of Stiff Ordinary Differential Equations at Low Tolerances" (Ekanathan et al., 2024)
• "High order schemes based on operator splitting and deferred corrections for stiff time dependent PDEs" (Duarte et al., 2014)