GitHub

OrdinaryDiffEqFIRK

Fully Implicit Runge-Kutta (FIRK) methods for stiff differential equations requiring very high accuracy. These methods solve a fully coupled implicit system at each timestep, providing superior accuracy and stability compared to diagonally implicit methods.

Real Numbers Only

FIRK methods should only be used for problems defined on real numbers, not complex numbers.

Key Properties

FIRK methods provide:

  • Highest-order implicit methods (excluding extrapolation)
  • Superior accuracy for very low tolerance requirements (≤ 1e-9)
  • A-stable and L-stable behavior for stiff problems
  • Higher order per stage than SDIRK methods (order 2s+1 for s stages)
  • Special geometric properties (some methods are symplectic)
  • Excellent for small to medium systems with high accuracy requirements

When to Use FIRK Methods

These methods are recommended for:

  • Very low tolerance problems (1e-9 and below) where accuracy is paramount
  • Small to medium stiff systems (< 200 equations)
  • Problems requiring highest possible accuracy for implicit methods
  • Stiff problems where SDIRK order limitations (max order 5) are insufficient
  • Applications where computational cost is acceptable for maximum accuracy

Mathematical Background

RadauIIA methods are based on Gaussian collocation and achieve order 2s+1 for s stages, making them among the highest-order implicit methods available. They represent the ODE analog of Gaussian quadrature. For more details on recent advances in FIRK methods, see our paper: High-Order Adaptive Time Stepping for the Incompressible Navier-Stokes Equations.

Computational Considerations

Advantages

  • Higher accuracy per stage than diagonal methods
  • Better multithreading for small systems due to larger linear algebra operations
  • No order restrictions like SDIRK methods (which max out at order 5)

Disadvantages

  • Limited to real-valued problems - cannot be used for complex number systems
  • Higher implementation complexity compared to SDIRK methods

Solver Selection Guide

High accuracy requirements

  • AdaptiveRadau: Recommended - adaptive order method that automatically selects optimal order
  • RadauIIA5: 5th-order method, good balance of accuracy and efficiency
  • RadauIIA9: 9th-order method for extremely high accuracy requirements
  • RadauIIA3: 3rd-order method for moderate accuracy needs

System size considerations

  • Systems < 200: FIRK methods are competitive due to better multithreading
  • Systems > 200: Consider SDIRK or BDF methods instead

Performance Guidelines

  • Best for tolerances ≤ 1e-9 where high accuracy justifies the cost
  • Most efficient on small to medium systems where linear algebra cost is manageable
  • Should be tested against parallel implicit extrapolation methods which specialize in similar regimes
  • Compare with high-order SDIRK methods for borderline cases

Installation

To be able to access the solvers in OrdinaryDiffEqFIRK, you must first install them use the Julia package manager:

using Pkg
Pkg.add("OrdinaryDiffEqFIRK")

This will only install the solvers listed at the bottom of this page. If you want to explore other solvers for your problem, you will need to install some of the other libraries listed in the navigation bar on the left.

Example usage

using OrdinaryDiffEqFIRK

function lorenz!(du, u, p, t)
    du[1] = 10.0 * (u[2] - u[1])
    du[2] = u[1] * (28.0 - u[3]) - u[2]
    du[3] = u[1] * u[2] - (8 / 3) * u[3]
end
u0 = [1.0; 0.0; 0.0]
tspan = (0.0, 100.0)
prob = ODEProblem(lorenz!, u0, tspan)
sol = solve(prob, RadauIIA5())

Full list of solvers

OrdinaryDiffEqFIRK.RadauIIA3 — Type
RadauIIA3(; autodiff = AutoForwardDiff(),
            concrete_jac = nothing,
            linsolve = nothing,
            extrapolant = :dense,
            smooth_est = true,
            step_limiter! = trivial_limiter!)

Fully-Implicit Runge-Kutta Method. An A-B-L stable fully implicit Runge-Kutta method with internal tableau complex basis transform for efficiency. Similar to Hairer's SEULEX.

Keyword Arguments

  • autodiff: Uses ADTypes.jl to specify whether to use automatic differentiation via ForwardDiff.jl or finite differencing via FiniteDiff.jl. Defaults to AutoForwardDiff() for automatic differentiation, which by default uses chunksize = 0, and thus uses the internal ForwardDiff.jl algorithm for the choice. To use FiniteDiff.jl, the AutoFiniteDiff() ADType can be used, which has a keyword argument fdtype with default value Val{:forward}(), and alternatives Val{:central}() and Val{:complex}().
  • concrete_jac: Specifies whether a Jacobian should be constructed. Defaults to nothing, which means it will be chosen true/false depending on circumstances of the solver, such as whether a Krylov subspace method is used for linsolve.
  • linsolve: Any LinearSolve.jl compatible linear solver. For example, to use KLU.jl, specify RadauIIA3(linsolve = KLUFactorization()). When nothing is passed, uses DefaultLinearSolver.
  • extrapolant: TBD
  • smooth_est: TBD
  • step_limiter!: function of the form limiter!(u, integrator, p, t)

References

@article{hairer1999stiff, title={Stiff differential equations solved by Radau methods}, author={Hairer, Ernst and Wanner, Gerhard}, journal={Journal of Computational and Applied Mathematics}, volume={111}, number={1-2}, pages={93–111}, year={1999}, publisher={Elsevier}}

source
OrdinaryDiffEqFIRK.RadauIIA5 — Type
RadauIIA5(; autodiff = AutoForwardDiff(),
            concrete_jac = nothing,
            linsolve = nothing,
            extrapolant = :dense,
            smooth_est = true,
            step_limiter! = trivial_limiter!)

Fully-Implicit Runge-Kutta Method. An A-B-L stable fully implicit Runge-Kutta method with internal tableau complex basis transform for efficiency. 5th order method with excellent numerical stability. Good for highly stiff systems, problems requiring high-order implicit integration, systems with complex eigenvalue structures. Best for low tolerance stiff problems (<1e-9).

Keyword Arguments

  • autodiff: Uses ADTypes.jl to specify whether to use automatic differentiation via ForwardDiff.jl or finite differencing via FiniteDiff.jl. Defaults to AutoForwardDiff() for automatic differentiation, which by default uses chunksize = 0, and thus uses the internal ForwardDiff.jl algorithm for the choice. To use FiniteDiff.jl, the AutoFiniteDiff() ADType can be used, which has a keyword argument fdtype with default value Val{:forward}(), and alternatives Val{:central}() and Val{:complex}().
  • concrete_jac: Specifies whether a Jacobian should be constructed. Defaults to nothing, which means it will be chosen true/false depending on circumstances of the solver, such as whether a Krylov subspace method is used for linsolve.
  • linsolve: Any LinearSolve.jl compatible linear solver. For example, to use KLU.jl, specify RadauIIA5(linsolve = KLUFactorization()). When nothing is passed, uses DefaultLinearSolver.
  • extrapolant: TBD
  • smooth_est: TBD
  • step_limiter!: function of the form limiter!(u, integrator, p, t)

References

@article{hairer1999stiff, title={Stiff differential equations solved by Radau methods}, author={Hairer, Ernst and Wanner, Gerhard}, journal={Journal of Computational and Applied Mathematics}, volume={111}, number={1-2}, pages={93–111}, year={1999}, publisher={Elsevier}}

source
OrdinaryDiffEqFIRK.RadauIIA9 — Type
RadauIIA9(; autodiff = AutoForwardDiff(),
            concrete_jac = nothing,
            linsolve = nothing,
            extrapolant = :dense,
            smooth_est = true,
            step_limiter! = trivial_limiter!)

Fully-Implicit Runge-Kutta Method. An A-B-L stable fully implicit Runge-Kutta method with internal tableau complex basis transform for efficiency. Similar to Hairer's SEULEX.

Keyword Arguments

  • autodiff: Uses ADTypes.jl to specify whether to use automatic differentiation via ForwardDiff.jl or finite differencing via FiniteDiff.jl. Defaults to AutoForwardDiff() for automatic differentiation, which by default uses chunksize = 0, and thus uses the internal ForwardDiff.jl algorithm for the choice. To use FiniteDiff.jl, the AutoFiniteDiff() ADType can be used, which has a keyword argument fdtype with default value Val{:forward}(), and alternatives Val{:central}() and Val{:complex}().
  • concrete_jac: Specifies whether a Jacobian should be constructed. Defaults to nothing, which means it will be chosen true/false depending on circumstances of the solver, such as whether a Krylov subspace method is used for linsolve.
  • linsolve: Any LinearSolve.jl compatible linear solver. For example, to use KLU.jl, specify RadauIIA9(linsolve = KLUFactorization()). When nothing is passed, uses DefaultLinearSolver.
  • extrapolant: TBD
  • smooth_est: TBD
  • step_limiter!: function of the form limiter!(u, integrator, p, t)

References

@article{hairer1999stiff, title={Stiff differential equations solved by Radau methods}, author={Hairer, Ernst and Wanner, Gerhard}, journal={Journal of Computational and Applied Mathematics}, volume={111}, number={1-2}, pages={93–111}, year={1999}, publisher={Elsevier}}

source
OrdinaryDiffEqFIRK.AdaptiveRadau — Type
AdaptiveRadau(; autodiff = AutoForwardDiff(),
                concrete_jac = nothing,
                linsolve = nothing,
                extrapolant = :dense,
                smooth_est = true,
                step_limiter! = trivial_limiter!)

Fully-Implicit Runge-Kutta Method. An A-B-L stable fully implicit Runge-Kutta method with internal tableau complex basis transform for efficiency. Fully autonomous derivation of tableau for arbitrary order and order adaptivity.

Keyword Arguments

  • autodiff: Uses ADTypes.jl to specify whether to use automatic differentiation via ForwardDiff.jl or finite differencing via FiniteDiff.jl. Defaults to AutoForwardDiff() for automatic differentiation, which by default uses chunksize = 0, and thus uses the internal ForwardDiff.jl algorithm for the choice. To use FiniteDiff.jl, the AutoFiniteDiff() ADType can be used, which has a keyword argument fdtype with default value Val{:forward}(), and alternatives Val{:central}() and Val{:complex}().
  • concrete_jac: Specifies whether a Jacobian should be constructed. Defaults to nothing, which means it will be chosen true/false depending on circumstances of the solver, such as whether a Krylov subspace method is used for linsolve.
  • linsolve: Any LinearSolve.jl compatible linear solver. For example, to use KLU.jl, specify AdaptiveRadau(linsolve = KLUFactorization()). When nothing is passed, uses DefaultLinearSolver.
  • extrapolant: TBD
  • smooth_est: TBD
  • step_limiter!: function of the form limiter!(u, integrator, p, t)

References

@article{AdaptiveRadauPaper, author={Ekanathan, Shreyas and Smith, Oscar and Rackauckas, Christopher}, booktitle={2025 IEEE High Performance Extreme Computing Conference (HPEC)}, title={A Fully Adaptive Radau Method for the Efficient Solution of Stiff Ordinary Differential Equations at Low Tolerances}, year={2025}, pages={1-9}, doi={10.1109/HPEC67600.2025.11196706}}

source