Non-autonomous Linear ODE / Lie Group ODE Solvers

Non-autonomous linear ODE solvers focus on equations in the general form of

and utilize the Lie group structure of the solution to accelerate the numerical methods and capture certain properties in the solution process. One common simplification is for solvers to require state-independent operators, which implies the form:

Another type of solver is needed when the operators are state-dependent, i.e.

Others specifically require linearity, i.e.

where is a constant operator.

Recommendations

It is recommended to always specialize on the properties of the operator as much as possible.

Standard ODE Integrators

The standard ODE integrators will work on Non-autonomous linear ODE problems via an automatic transformation to a first-order ODE. See the ODE solvers page for more details.

Specialized OrdinaryDiffEq.jl Integrators

Unless otherwise specified, the OrdinaryDiffEq algorithms all come with a 3rd order Hermite polynomial interpolation. The algorithms denoted as having a “free” interpolation means that no extra steps are required for the interpolation. For the non-free higher order interpolating functions, the extra steps are computed lazily (i.e. not during the solve).

Note that all of these methods are fixed timestep unless otherwise specified.

Exponential Methods for Linear and Affine Problems

These methods require that is constant.

LinearExponential - Exact solution formula for linear, time-independent problems.

Options:

krylov - symbol. One of
- :off (default) - cache the operator beforehand. Requires Matrix(A) method defined for the operator A.
- :simple - uses simple Krylov approximations with fixed subspace size m.
- :adaptive - uses adaptive Krylov approximations with internal timestepping.
m - integer, default: 30. Controls the size of Krylov subspace if krylov=:simple, and the initial subspace size if krylov=:adaptive.
iop - integer, default: 0. If not zero, determines the length of the incomplete orthogonalization procedure (IOP) ^[1]. Note that if the linear operator/Jacobian is hermitian, then the Lanczos algorithm will always be used and the IOP setting is ignored.

using DifferentialEquations
_A = [2 -1; -3 -5] / 5
A = DiffEqArrayOperator(_A)
prob = ODEProblem(A, [1.0, -1.0], (1.0, 6.0))
sol = solve(prob, LinearExponential())

retcode: Success
Interpolation: 3rd order Hermite
t: 2-element Vector{Float64}:
 1.0
 6.0
u: 2-element Vector{Vector{Float64}}:
 [1.0, -1.0]
 [11.923375744981003, -4.833128734161089]

LinearExponential is exact, and thus it uses dt=tspan[2]-tspan[1] by default. The interpolation used is inexact (3rd order Hermite). Thus, values generated by the interpolation (via sol(t) or saveat) will be inexact with increasing error as the size of the time span grows. To counteract this, directly set dt or use tstops instead of saveat. For more information, see this issue.

State-Independent Solvers

These methods require is only dependent on the independent variable, i.e. .

MagnusMidpoint - Second order Magnus Midpoint method.
MagnusLeapfrog- Second order Magnus Leapfrog method.
MagnusGauss4 - Fourth order Magnus method approximated using a two stage Gauss quadrature.
MagnusGL4- Fourth order Magnus method approximated using Gauss-Legendre quadrature.
MagnusNC6- Sixth order Magnus method approximated using Newton-Cotes quadrature.
MagnusGL6- Sixth order Magnus method approximated using Gauss-Legendre quadrature.
MagnusNC8- Eighth order Magnus method approximated using Newton-Cotes quadrature.
MagnusGL8- Eighth order Magnus method approximated using Gauss-Legendre quadrature.

Example:

function update_func(A, u, p, t)
    A[1, 1] = cos(t)
    A[2, 1] = sin(t)
    A[1, 2] = -sin(t)
    A[2, 2] = cos(t)
end
A = DiffEqArrayOperator(ones(2, 2), update_func = update_func)
prob = ODEProblem(A, ones(2), (1.0, 6.0))
sol = solve(prob, MagnusGL6(), dt = 1 / 10)

retcode: Success
Interpolation: 3rd order Hermite
t: 51-element Vector{Float64}:
 1.0
 1.1
 1.2000000000000002
 1.3000000000000003
 1.4000000000000004
 1.5000000000000004
 1.6000000000000005
 1.7000000000000006
 1.8000000000000007
 1.9000000000000008
 ⋮
 5.199999999999998
 5.299999999999998
 5.399999999999998
 5.499999999999997
 5.599999999999997
 5.699999999999997
 5.799999999999996
 5.899999999999996
 6.0
u: 51-element Vector{Vector{Float64}}:
 [1.0, 1.0]
 [0.9560322602141308, 1.138059339025732]
 [0.8837222756549167, 1.2712953135220313]
 [0.7836511349460551, 1.3924887039087297]
 [0.6585912321071351, 1.4945430086037101]
 [0.5134387044387002, 1.5712483084024886]
 [0.35484956094894243, 1.6179942444961344]
 [0.190617494783206, 1.6323079788211143]
 [0.028885769775000115, 1.6141203517260592]
 [-0.12267917859524063, 1.5657147744869344]
 ⋮
 [0.1649471518131617, 0.19050761906534755]
 [0.19016510949800752, 0.18488701621917145]
 [0.21691662630888853, 0.17939687440026209]
 [0.24554653627671552, 0.17419950752424784]
 [0.27642938801560896, 0.16953149318836871]
 [0.3099648774418323, 0.16572504624198953]
 [0.3465668460745136, 0.16323354028339193]
 [0.3866431012103112, 0.1626603704342466]
 [0.4305628420420221, 0.16478922184291447]

The initial values for are irrelevant in this and similar cases, as the update_func immediately overwrites them. Starting with ones(2,2) is just a convenient way to get a mutable 2x2 matrix.

State-Dependent Solvers

These methods can be used when is dependent on the state variables, i.e. .

CayleyEuler - First order method using Cayley transformations.
LieEuler - First order Lie Euler method.
RKMK2 - Second order Runge—Kutta—Munthe-Kaas method.
RKMK4 - Fourth order Runge—Kutta—Munthe-Kaas method.
LieRK4 - Fourth order Lie Runge-Kutta method.
CG2 - Second order Crouch—Grossman method.
MagnusAdapt4 - Fourth Order Adaptive Magnus method.

Example:

function update_func(A, u, p, t)
    A[1, 1] = 0
    A[2, 1] = sin(u[1])
    A[1, 2] = -1
    A[2, 2] = 0
end
A = DiffEqArrayOperator(ones(2, 2), update_func = update_func)
prob = ODEProblem(A, ones(2), (0, 30.0))
sol = solve(prob, LieRK4(), dt = 1 / 4)

retcode: Success
Interpolation: 3rd order Hermite
t: 121-element Vector{Float64}:
  0.0
  0.25
  0.5
  0.75
  1.0
  1.25
  1.5
  1.75
  2.0
  2.25
  ⋮
 28.0
 28.25
 28.5
 28.75
 29.0
 29.25
 29.5
 29.75
 30.0
u: 121-element Vector{Vector{Float64}}:
 [1.0, 1.0]
 [0.7273968579551914, 1.165794872747704]
 [0.4243104556165744, 1.2459210126860043]
 [0.10953160297897346, 1.2654872693089263]
 [-0.2070137554757648, 1.2681541526375768]
 [-0.5272224527692214, 1.3028146975278558]
 [-0.8648991354908154, 1.4149945285244698]
 [-1.244399456687535, 1.642767157085863]
 [-1.6975604074674677, 2.0038463754183535]
 [-2.2537711317460407, 2.449176910606262]
 ⋮
 [-2.6663888753291105, 2.6954837342889535]
 [-3.3611719087644736, 2.7811209839661157]
 [-4.008187561635523, 2.2967303664591197]
 [-4.470316530373341, 1.3459658518794193]
 [-4.66626299411022, 0.2081263142861722]
 [-4.57289678174422, -0.9494989022112654]
 [-4.200152270252206, -1.99585038589234]
 [-3.6052140155232593, -2.6767378307341736]
 [-2.9109618325260116, -2.781287635056373]

The above example solves a non-stiff Non-Autonomous Linear ODE with a state dependent operator, using the LieRK4 method. Similarly, a stiff Non-Autonomous Linear ODE with state dependent operators can be solved using specialized adaptive algorithms, like MagnusAdapt4.

Example:

function update_func(A, u, p, t)
    A[1, 1] = 0
    A[2, 1] = 1
    A[1, 2] = -2 * (1 - cos(u[2]) - u[2] * sin(u[2]))
    A[2, 2] = 0
end
A = DiffEqArrayOperator(ones(2, 2), update_func = update_func)
prob = ODEProblem(A, ones(2), (30, 150.0))
sol = solve(prob, MagnusAdapt4())

retcode: Success
Interpolation: 3rd order Hermite
t: 592-element Vector{Float64}:
  30.0
  30.051554410975918
  30.14327232699745
  30.247524724052145
  30.375696909344892
  30.510855783890896
  30.65090569337366
  30.788186122035587
  30.91952329587291
  31.04078864675185
   ⋮
 147.96273686915632
 148.1386523776747
 148.34620600715417
 148.59559211110502
 148.8847559565441
 149.1635807149984
 149.44955972217943
 149.73645683847397
 150.0
u: 592-element Vector{Vector{Float64}}:
 [1.0, 1.0]
 [1.0418828143197263, 1.0526121015605732]
 [1.1297712835913454, 1.1520649699813994]
 [1.2527691803606753, 1.276033498183453]
 [1.4404455264543194, 1.4481983217458168]
 [1.6784196223822325, 1.6585827542112963]
 [1.941152341880844, 1.9121745307485782]
 [2.1365107433203687, 2.1934881055190636]
 [2.1250337632314107, 2.476421374656319]
 [1.7927262923353107, 2.7176941957849485]
 ⋮
 [1.3382785710099845, -1.3478235487981378]
 [1.1224319158310865, -1.132482377423453]
 [0.9547031935838253, -0.918358952842659]
 [0.8437842085519446, -0.6957609582428903]
 [0.7910831565788868, -0.4608017890713062]
 [0.7783134037066098, -0.24250843704692368]
 [0.7772133962093504, -0.020168014665712365]
 [0.7777533743732511, 0.20284549420577422]
 [0.7858794742201944, 0.4085646081595665]

Time and State-Dependent Operators

These methods can be used when is dependent on both time and state variables, i.e.

CG3 - Third order Crouch-Grossman method.

1. A description of IOP can be found in this paper.