---
title: 'Radau IIA Method: High-Order Implicit Integration'
url: https://www.emergentmind.com/topics/radau-iia-method
type: topic
---

# Radau IIA Method: High-Order Implicit Integration

The Radau IIA method is a family of implicit Runge–Kutta schemes of collocation type with stiffly accurate Radau quadrature abscissae, distinguished by their order $p=2s-1$ ($s$ stages), A- and L-stability, and systematic advantages for integrating stiff ODEs, DAEs, and optimal control problems at high accuracy thresholds. These methods operate by collocating the solution at the right-hand Radau points on $[0,1]$, including the endpoint, resulting in unique tableau formulas tailored for high-order time integration, error estimation, and robust nonlinear solver design.

## 1. Mathematical Definition and Order Properties

The $s$-stage Radau IIA method for solving $dy/dt=f(t,y)$ with step size $h$ is defined by the fully implicit Runge–Kutta system:
$$
Y_i = y_n + h\sum_{j=1}^s a_{ij}f(t_n + c_j h, Y_j),\quad i = 1,\dots,s
$$
$$
y_{n+1} = y_n + h\sum_{i=1}^s b_i f(t_n + c_i h, Y_i)
$$
with Butcher tableau
\[
\begin{array}{c|ccc}
c_1 & a_{11} & \cdots & a_{1s} \\
\vdots & \vdots & & \vdots \\
c_s & a_{s1} & \cdots & a_{ss} \\ \hline
 & b_1 & \cdots & b_s
\end{array}
\]
where $c_i$ are the roots of the shifted Radau polynomial condition $d^{s-1}/dx^{s-1}[x^{s-1}(x-1)^s]=0$ and $c_s=1$. The collocation theory ensures order $p=2s-1$ due to degree-$s$ quadrature with endpoint inclusion; this high order is proven both algebraically and via the precision properties of the associated quadrature and stage order. Radau IIA methods are stiffly accurate (i.e., $b_i=a_{s,i}$) and guarantee that the value at $c_s=1$ coincides with the new step value [2412.14362], [1604.00355], [2310.12846].

## 2. Stability Analysis

Radau IIA methods possess unconditional A-stability and L-stability. The stability function is
$$
R(z) = 1 + z\, \mathbf{b}^T (I - zA)^{-1}\mathbf{1}
$$
which matches the subdiagonal Padé$(s,s-1)$ approximation to $e^z$. The denominator polynomial $Q_s(z)=\det(I-zA)$ has all zeros in $\Re(z)<0$; thus the method's region of absolute stability contains the entire left half-plane. L-stability is achieved via $\lim_{z\to\infty}R(z)=0$, making Radau IIA strongly favorable for highly stiff problems and damping of fast transient modes [1604.00355], [2412.14362], [1311.0640].

## 3. Tableau Construction and Coefficient Derivation

The tableau coefficients can be computed on the fly for arbitrary order and precision, not limited to double precision or a small set of fixed orders. Given $s=(p+1)/2$, the nodes $c_i$ are obtained from the Radau polynomial, internal coefficients $a_{ij}$ solve
$$
\sum_{j=1}^s a_{ij}c_j^{q-1} = {c_i^q}/{q},\quad q=1,\dots,s
$$
and $b_i=a_{s,i}$ for the weights. Practically, polynomial root-finding and Vandermonde system solution yield $(c, a, b)$, enabling real-time derivation even in arbitrary precision arithmetic (e.g., BigFloat). This tableau flexibility is central to fully adaptive algorithms at low error tolerances and high accuracy [2412.14362].

| Order $p$ | Stages $s$ | Endpoint node $c_s$ | Global order | Stiffly accurate |
|:-------:|:------:|:----------------:|:-----------:|:---------------:|
| 3       | 2      | 1                | 3           | Yes             |
| 5       | 3      | 1                | 5           | Yes             |
| 9       | 5      | 1                | 9           | Yes             |
| 13      | 7      | 1                | 13          | Yes             |

## 4. Embedded Error Estimation and Time-Step Control

Radau IIA implementations typically employ embedded lower-order formulas for adaptive time step and order control. Given the main weights $b_i$ for order $p$, one constructs new weights $\hat{b}_i$ for order $p-1$ by
$$
\sum_{i=1}^s \hat{b}_i c_i^{k-1} = 1/k,\quad k=1,\dots,s
$$
The step error estimate is
$$
e_{n+1} = y_{n+1}^{(p)} - y_{n+1}^{(p-1)}
$$
with a normed error control factor,
$$
E = \|e_{n+1}/(atol + \max(|y_n|, |y_{n+1}|) rtol)\|_\infty
$$
Step size is updated via a predictive controller
$$
h_{new} = h \cdot \min(5, \max(0.2, E^{-1/(p+1)}\cdot f)),\ \ f = ((1+2m)\gamma)/(iter+2m),\ \ \gamma\approx0.8
$$
where $m$ is the Newton iteration budget and $iter$ the actual count. This design provides robust adaptivity for stiff problems, with further order adaptation driven by a low-pass filtered Newton iteration history [2412.14362], [1604.00355].

## 5. Nonlinear and Linear Algebra Procedures

Each step requires solving the coupled nonlinear system for stage values, typically via simplified/sparse Newton methods, sometimes leveraging GMRES or low-rank splitting. The stage Jacobian linear solve in classical implementations scales as $O(sn^3)$, but modern algorithms exploit block-diagonal structures in $A^{-1}$, performing similarity transformations to split into one $n \times n$ system and several $2n \times 2n$ systems (for complex eigenvalue pairs). This enables smaller LU factorizations, improved cache usage, and vectorized SIMD computations. Efficient low-rank splitting (auxiliary abscissae and Crout factorization) further reduces work per step to one $m \times m$ LU and $O(s m^2)$ back-substitutions, yielding substantial speed improvements and rapid Newton convergence [2412.14362], [1302.1037], [1604.00355].

## 6. Applications to Stiff ODEs, DAEs, PDEs, and Optimal Control

Radau IIA methods are among the canonical integrators for extremely stiff ODEs, index-1 and index-2 DAEs, and large-scale semidiscrete PDE systems. In DAEs, the stiffly accurate stage makes endpoint algebraic constraints tractable. In optimal control, Radau IIA collocation ensures commutation of discretization and optimization, delivering full order under additional “A-conditions,” which are automatically satisfied by the underlying tableau [1508.03783], [1311.0640]. In high-index physical DAEs, Radau IIA, embedded within physics-informed neural networks, has yielded state-of-the-art accuracy for both differential and algebraic states, especially at fifth order [2310.12846]. In adaptive PDE solvers, Radau IIA methods allow compressed representations and efficient time integration with adaptive spatial grids [1604.00355]. 

| Application Domain         | Property/Advantage                 | Reference         |
|:--------------------------|:-----------------------------------|:------------------|
| Stiff ODEs                | High order, A-/L-stability         | [2412.14362]      |
| Index-1, Index-2 DAEs     | Endpoint, stiffly accurate         | [2212.02630], [2310.12846] |
| PDEs on adaptive grids     | Sparse block Jacobians, efficient step control | [1604.00355]      |
| Optimal control           | Exponential convergence, commutation of optimize/discretize | [1508.03783], [1311.0640] |

## 7. Modern Implementations and Performance

Recent algorithms provide adaptive-time adaptive-order Radau IIA with tableau generation at arbitrary order and precision (including BigFloat arithmetic), improved tableau caching to avoid recomputation, and advanced linear algebra with block decoupling. Benchmark results show 2x speedups over classic Fortran implementations (e.g., Hairer's RADAU5) and 1.5x acceleration over advanced BDF and Rosenbrock solvers at low error tolerances; in very high precision (≤1e−24), Radau IIA outperforms popular alternatives due to its superconvergent order scaling and mixed-precision capability [2412.14362]. Practical solvers (e.g., AdaptiveRadau in Julia) exploit automatic Jacobian differentiation and structure-aware sparse algebra for large-scale problems. Maple and PINN architectures leverage Radau IIA for direct algebraic constraint enforcement and analytic Jacobian extraction [2212.02630], [2310.12846]. For optimal control and DAE systems, observed convergence rates match theoretical expectations, with high stability and accuracy at orders 3 and 5 [1311.0640], [2310.12846].

## 8. Summary of Critical Properties

Radau IIA methods:
- Possess order $p=2s-1$ with endpoint collocation.
- Are A-stable and L-stable, suitable for extreme stiffness.
- Enable robust adaptive time–order stepping, with embedded error control.
- Benefit from block-diagonal and low-rank splitting for scalable nonlinear and linear algebra.
- Are stiffly accurate, optimal for DAEs and control problems.
- Support arbitrary order and precision in modern implementations.
- Deliver superior performance and accuracy across stiff ODEs, DAEs, PDEs, and optimal control scenarios.

All undeclared metrics, algorithms, and comparative statements are traceable directly to the cited literature [2412.14362], [1604.00355], [1302.1037], [2310.12846], [1508.03783], [1311.0640], [2212.02630].

Source: https://www.emergentmind.com/topics/radau-iia-method