Stable corrections for perturbed diagonally implicit Runge--Kutta methods

1. Problem Statement

The paper addresses the computational challenges associated with solving stiff systems of Ordinary Differential Equations (ODEs) using Diagonally Implicit Runge–Kutta (DIRK) methods.

• The Bottleneck: DIRK methods are preferred for stiff problems due to their large linear stability regions, but they require costly implicit solves at each stage.
• The Motivation: To accelerate computation, researchers have proposed mixed-accuracy frameworks where the implicit stages are solved using lower-precision arithmetic, under-resolved iterative solvers, or simplified linearized models.
• The Challenge: While these low-accuracy approaches reduce computational cost, they introduce perturbation errors ($h(y) = f(y) - f_\epsilon(y)$). Previous work established that these errors could be damped out to maintain accuracy, but explicit corrections proposed to further improve accuracy were found to introduce instabilities when the time-step ($\Delta t$) is large.
• The Goal: The authors aim to rigorously analyze the stability of these perturbed methods and design stabilized correction strategies that maintain high accuracy without sacrificing stability, even for larger time-steps.

2. Methodology

A. Theoretical Framework

The authors analyze a perturbed DIRK method where the implicit stage function $f$ is replaced by a low-accuracy approximation $f_\epsilon$.

• Error Growth Analysis: They derive bounds on the error growth between the exact solution ($z$) and the perturbed solution ($y$). Assuming the function $f$ is contractive (a common property for stiff problems), they prove that the error growth per time-step is bounded by terms involving the perturbation size $\epsilon$, the stiffness $L$, and the time-step $\Delta t$.
• Key Theorem (Theorem 1): They establish that the global error at time $T_f$ grows linearly with time as $O(\epsilon \Delta t L T_f)$. This confirms that the method is stable (errors do not grow exponentially) provided the perturbation is small relative to the stiffness and time-step.
• Limitation of Previous Corrections: They analyze the explicit corrections proposed in prior literature. While these improve accuracy by $O(\Delta t)$ per correction step, they act as fixed-point iterations that converge only if $\alpha \Delta t L \leq 1$. For stiff problems (large $L$) or large $\Delta t$, this condition is violated, leading to instability.

B. Stabilized Correction Strategy

To overcome the instability of explicit corrections, the authors propose a stabilized correction approach.

• The Mechanism: Instead of a pure explicit update, they introduce a stabilization matrix $\Phi$ (or equivalently a matrix $J$) into the correction iteration:
  $$y^{[k+1]} = y^{[k]} + \Phi \left( y_{exp} + \alpha \Delta t f(y^{[k]}) - y^{[k]} \right)$$
  where $\Phi = (I - \mu \Delta t J)^{-1}$.
• Stability Condition: The matrix $J$ is chosen to approximate the Jacobian of the system. The term $(I - \mu \Delta t J)^{-1}$ acts as a damping factor, ensuring the iteration remains stable even when $\Delta t$ is large.
• Strategies for Defining $\Phi$:
  1. Static Jacobian ($\Phi_J$): Computed once using the Jacobian at the initial time ($J_0 = f'(y_0)$).
  2. Static Differential Operator ($\Phi_{EIN}$): Based on the dominant spatial derivative operator (inspired by the Explicit-Implicit-Null method), e.g., $L = -D_x$ for advection or $L = D_{xx}$ for diffusion.
  3. Dynamic Broyden Update ($\Phi_B$): Uses a Broyden-type quasi-Newton update to approximate the inverse Jacobian dynamically without full recomputation, updating $\Phi$ at each time-step or stage.

3. Key Contributions

1. Rigorous Stability Analysis: The paper provides a formal proof (Lemma 1, Lemma 2, Theorem 1) bounding the error growth of perturbed DIRK methods, demonstrating that errors grow linearly rather than exponentially for contractive systems.
2. Novel Stabilized Corrections: The authors introduce a framework for correcting mixed-accuracy DIRK methods that includes a stabilization matrix $\Phi$. This prevents the instabilities associated with pure explicit corrections in stiff regimes.
3. Comparative Evaluation of Stabilization Matrices: They systematically compare three approaches for defining $\Phi$:
   • Frozen Jacobian ($\Phi_J$)
   • Differential Operator ($\Phi_{EIN}$)
   • Dynamic Broyden updates ($\Phi_B$)
4. Comprehensive Numerical Validation: The methods are tested on three distinct nonlinear problems:
   • Inviscid Burgers' Equation: Testing linearization and mixed-precision effects.
   • Shallow Water Equations: A system of equations testing coupled dynamics.
   • Porous Medium Equation: A highly nonlinear diffusion problem.

4. Results

• Accuracy vs. Stability Trade-off:

  • No Corrections: Mixed-accuracy methods degrade to first-order accuracy due to perturbation errors.
  • Explicit Corrections: Improve accuracy to the design order for small $\Delta t$ but cause catastrophic instability for larger $\Delta t$ (especially in higher-order methods like SDIRK4).
  • Stabilized Corrections: Successfully recover the design order of accuracy (2nd, 3rd, or 4th order) while maintaining stability across a much wider range of $\Delta t$.

• Performance of Stabilization Matrices:

  • $\Phi_J$ (Frozen Jacobian): Consistently provided the best balance of stability and accuracy, particularly for the porous medium equation and mixed-precision implementations.
  • $\Phi_{EIN}$ (Differential Operator): Performed well and was comparable to explicit corrections in some cases but generally less effective than $\Phi_J$ for highly nonlinear problems.
  • $\Phi_B$ (Broyden): While theoretically promising for dynamic updates, the "Bad Broyden" approach performed similarly to static methods in the tested scenarios, and the added computational overhead did not always yield significant gains over the static Jacobian approach.

• Mixed Precision Impact:

  • In mixed-precision implementations (e.g., double/single or quad/double), the perturbation error $\epsilon$ depends on the system size ($N_x$). As $N_x$ increases, round-off errors accumulate, degrading accuracy.
  • The stabilized corrections effectively dampened these accumulated errors, allowing high-order methods to function correctly even with significant precision reduction in the implicit solves.

5. Significance

This work is significant for the computational mathematics community for several reasons:

• Enabling Efficient Stiff Solvers: It provides a theoretical and practical pathway to use low-precision arithmetic (e.g., FP16 or FP32) for the most expensive parts of stiff ODE/PDE solvers without sacrificing the stability or accuracy required for scientific simulations.
• Robustness: By solving the instability issue inherent in previous correction schemes, the proposed method allows for larger time-steps, which is critical for long-time integration of stiff systems.
• Hardware Efficiency: The ability to use lower precision or simpler models for implicit solves aligns with modern hardware trends (GPUs, TPUs) that favor lower precision for throughput, potentially offering significant speedups for large-scale simulations in fluid dynamics, climate modeling, and reaction-diffusion systems.
• Generalizability: The framework is applicable to any contractive system and can be adapted to various types of perturbations, including linearization errors, iterative solver truncation, and hardware-induced round-off.

In conclusion, the authors successfully bridge the gap between computational efficiency (via mixed-accuracy) and numerical robustness (via stabilized corrections), offering a viable strategy for next-generation high-performance computing of stiff differential equations.