ScalarRadau
This module implements the 5th order, Radau IIA method for a scalar ordinary differential equation (ODE), in Julia. The algorithm is famously effective for stiff ODEs. Implementation mostly follows the description in chapter IV.8 in Solving Ordinary Differential Equations II, by Ernst Hairer and Gerhard Wanner, with a couple small changes that were found to be beneficial for scalar equations.
- Step size is adaptive and the initial step size is chosen automatically.
- Functions implemented here expect to use
Float64
numbers. - Dense output for continuous solutions is implemented using cubic Hermite interpolation.
- Approximate Jacobian evaluation is performed with a finite difference.
- Because the equation is scalar and the method has three stages, the Jacobian is always a 3 x 3 matrix. Static arrays are used for efficient Newton iterations.
The implementation here is useful for any scenario where a stiff, scalar ODE must be solved repeatedly under different conditions. For example, you might need to solve the same stiff ODE with a range of different initial conditions or with many sets of system parameters. The solver functions specialize directly on the ODE provided, making them fast. This is slightly different than DifferentialEquations.jl, which uses a two-step system of defining an ODE problem with one function then solving it with another function, but if you need to solve a stiff system of ODEs instead of a scalar equation, look here.
For a nice overview of Radau methods and their utility, check out: Stiff differential equations solved by Radau methods.
Quick Start
Install using Julia's package manager
julia> ]add ScalarRadau
To solve an ODE, first define it as a function, then pass it to the radau
function.
using ScalarRadau
F(x, y, param) = -y
x, y = radau(F, 1, 0, 2, 25)
The snippet above solves the equation dy/dx = -y
, starting at y=1
, between x=0
and x=2
, and returns 25 evenly spaced points in the solution interval.
In-place Solution
For full control over output points, the in-place function is
radau!(yout, xout, F, y₀, x₀, xₙ, param=nothing; rtol=1e-6, atol=1e-6, facmax=100.0, facmin=0.01, κ=1e-3, ϵ=0.25, maxnwt=7, maxstp=1000000)
Mandatory function arguments are
yout
- vector where output points will be writtenxout
- sorted vector ofx
values where output points should be sampledF
- scalar ODE in the formF(x, y, param)
y₀
- initial value fory
x₀
- starting point forx
xₙ
- end point of the integration
The optional param
argument is nothing
by default, but it may be any type and is meant for scenarios where extra information must be accessible to the ODE function. It is passed to your ODE function whenever it's evaluated.
The coordinates of the output points (xout
) should be between x₀
and xₙ
and they should be in ascending order. They are not checked for integrity before integrating, though. The only check performed is xₙ > x₀
, or that the integration isn't going backward.
Keyword arguments are
rtol
- relative error toleranceatol
- absolute error tolerancefacmax
- maximum fraction that the step size may increase, compared to the previous stepfacmin
- minimum fraction that the step size may decrease, compared to the previous stepκ
(kappa) - stopping tolerance for Newton iterationsϵ
(epsilon) - fraction of current step size used for finite difference Jacobian approximationmaxnwt
- maximum number of Newton iterations before step size reductionmaxstp
- maximum number of steps for the solver stops and throws an error
Two other functions are available to make different output options convenient. Both of them use the in-place version internally.
Evenly Spaced Output
For evenly spaced output points (as in the example above) the function definition is
radau(F, y₀, x₀, xₙ, nout, param=nothing; kwargs...)
In this case, you must specify the number of output points with the nout
argument. Keyword arguments and default values are the same as above. Solution vectors for x
and y
are returned.
End-point Output
To compute only the y
value at the end of the integration interval (xₙ
), the function is
radau(F, y₀, x₀, xₙ, param=nothing; kwargs...)
Again, keyword arguments and default values are identical to the in-place function.