---
title: Stiffly Accurate Radau Quadrature
url: https://www.emergentmind.com/topics/stiffly-accurate-radau-quadrature
type: topic
---

# Stiffly Accurate Radau Quadrature

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 [2412.14362], [1407.0195].

## 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 [1407.0195].

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$ [2412.14362], [1407.0195].

## 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$ [2412.14362], [1407.0195].

## 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}$ [2412.14362].
- 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 [2412.14362].

## 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 [2412.14362].

## 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 [2412.14362].

## 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 [1407.0195].

## 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) [2412.14362].
- 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" [2412.14362]  
- "High order schemes based on operator splitting and deferred corrections for stiff time dependent PDEs" [1407.0195]

Source: https://www.emergentmind.com/topics/stiffly-accurate-radau-quadrature