OrdinaryDiffEqRosenbrock

Rosenbrock methods for mass matrix differential-algebraic equations (DAEs) and stiff ODEs with singular mass matrices. These methods provide efficient integration for moderately stiff systems with algebraic constraints, offering excellent performance for small to medium-sized DAE problems.

Key Properties

Mass matrix Rosenbrock methods provide:

DAE capability for index-1 differential-algebraic equations
W-method efficiency using approximate Jacobians for computational savings
Mass matrix support for singular and non-diagonal mass matrices
Moderate to high order accuracy (2nd to 6th order available)
Good stability properties with stiffly accurate behavior
Embedded error estimation for adaptive timestepping

When to Use Mass Matrix Rosenbrock Methods

These methods are recommended for:

Index-1 DAE systems with moderate stiffness
Small to medium constrained systems (< 1000 equations)
Semi-explicit DAEs arising from discretized PDEs
Problems requiring good accuracy with moderate computational cost
DAEs with moderate nonlinearity where W-methods are efficient
Electrical circuits and mechanical systems with constraints

Warn

In order to use OrdinaryDiffEqRosenbrock with DAEs that require a non-trivial consistent initialization, a nonlinear solver is required and thus using OrdinaryDiffEqNonlinearSolve is required or you must pass an initializealg with a valid nlsolve choice.

Mathematical Background

Mass matrix DAEs have the form: M du/dt = f(u,t)

Rosenbrock methods linearize around the current solution and solve linear systems of the form: (M/γh - J) k_i = ...

where J is the Jacobian of f and γ is a method parameter.

Solver Selection Guide

Recommended Methods by Tolerance

High tolerances (>1e-2): Rosenbrock23 - efficient low-order method
Medium tolerances (1e-8 to 1e-2): Rodas5P - most efficient choice, or Rodas4P for higher reliability
Low tolerances (<1e-8): Rodas5Pe or higher-order alternatives

Method families

Rodas5P: Recommended - Most efficient 5th-order method for general use
Rodas4P: More reliable 4th-order alternative
Rosenbrock23: Good for high tolerance problems
Rodas5: Standard 5th-order method without embedded pair optimization

Performance Guidelines

When mass matrix Rosenbrock methods excel

Small to medium DAE systems (< 1000 equations)
Moderately stiff problems where full BDF methods are overkill
Problems with efficient Jacobian computation or finite difference approximation
Index-1 DAEs with well-conditioned mass matrices
Semi-explicit index-1 problems from spatial discretizations

System size considerations

Small systems (< 100): Rosenbrock methods often outperform multistep methods
Medium systems (100-1000): Good performance with proper linear algebra
Large systems (> 1000): Consider BDF methods instead

Important DAE Considerations

Initial conditions

Must be consistent with algebraic constraints
Consistent initialization may require nonlinear solver
Index-1 assumption for reliable performance

Mass matrix requirements

Index-1 DAE structure for optimal performance
Non-singular leading submatrix for differential variables
Well-conditioned constraint equations

Alternative Approaches

Consider these alternatives:

Mass matrix BDF methods for larger or highly stiff DAE systems
Implicit Runge-Kutta methods for higher accuracy requirements
Standard Rosenbrock methods for regular ODEs without constraints
IMEX methods if natural explicit/implicit splitting exists

Example Usage

using LinearAlgebra: Diagonal
function rober(du, u, p, t)
    y₁, y₂, y₃ = u
    k₁, k₂, k₃ = p
    du[1] = -k₁ * y₁ + k₃ * y₂ * y₃
    du[2] = k₁ * y₁ - k₃ * y₂ * y₃ - k₂ * y₂^2
    du[3] = y₁ + y₂ + y₃ - 1
    nothing
end
M = Diagonal([1.0, 1.0, 0])  # Singular mass matrix
f = ODEFunction(rober, mass_matrix = M)
prob_mm = ODEProblem(f, [1.0, 0.0, 0.0], (0.0, 1e5), (0.04, 3e7, 1e4))
sol = solve(prob_mm, Rodas5(), reltol = 1e-8, abstol = 1e-8)

Installation

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

using Pkg
Pkg.add("OrdinaryDiffEqRosenbrock")

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 OrdinaryDiffEqRosenbrock

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, Rodas5P())

Full list of solvers

OrdinaryDiffEqRosenbrock.Rosenbrock23 — Type

Rosenbrock23(; chunk_size = Val{0}(),
               standardtag = Val{true}(),
               autodiff = AutoForwardDiff(),
               concrete_jac = nothing,
               diff_type = Val{:forward}(),
               linsolve = nothing,
               precs = DEFAULT_PRECS,
               step_limiter! = OrdinaryDiffEq.trivial_limiter!)

Rosenbrock-Wanner-W(olfbrandt) Method. An Order 2/3 L-Stable Rosenbrock-W method which is good for very stiff equations with oscillations at low tolerances. 2nd order stiff-aware interpolation.

Keyword Arguments

standardtag: Specifies whether to use package-specific tags instead of the ForwardDiff default function-specific tags. For more information, see this blog post. Defaults to Val{true}().
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 Rosenbrock23(linsolve = KLUFactorization()). When nothing is passed, uses DefaultLinearSolver.
precs: Any LinearSolve.jl-compatible preconditioner can be used as a left or right preconditioner. Preconditioners are specified by the Pl,Pr = precs(W,du,u,p,t,newW,Plprev,Prprev,solverdata) function where the arguments are defined as:
- W: the current Jacobian of the nonlinear system. Specified as either $I - \gamma J$ or $I/\gamma - J$ depending on the algorithm. This will commonly be a WOperator type defined by OrdinaryDiffEq.jl. It is a lazy representation of the operator. Users can construct the W-matrix on demand by calling convert(AbstractMatrix,W) to receive an AbstractMatrix matching the jac_prototype.
- du: the current ODE derivative
- u: the current ODE state
- p: the ODE parameters
- t: the current ODE time
- newW: a Bool which specifies whether the W matrix has been updated since the last call to precs. It is recommended that this is checked to only update the preconditioner when newW == true.
- Plprev: the previous Pl.
- Prprev: the previous Pr.
- solverdata: Optional extra data the solvers can give to the precs function. Solver-dependent and subject to change.
The return is a tuple (Pl,Pr) of the LinearSolve.jl-compatible preconditioners. To specify one-sided preconditioning, simply return nothing for the preconditioner which is not used. Additionally, precs must supply the dispatch:
```
Pl, Pr = precs(W, du, u, p, t, ::Nothing, ::Nothing, ::Nothing, solverdata)
```
which is used in the solver setup phase to construct the integrator type with the preconditioners (Pl,Pr). The default is precs=DEFAULT_PRECS where the default preconditioner function is defined as:
```
DEFAULT_PRECS(W, du, u, p, t, newW, Plprev, Prprev, solverdata) = nothing, nothing
```
step_limiter!: function of the form limiter!(u, integrator, p, t)

References

Shampine L.F. and Reichelt M., (1997) The MATLAB ODE Suite, SIAM Journal of

Scientific Computing, 18 (1), pp. 1-22.