LSODA.jl

Introduction

LSODA.jl is a Julia package that interfaces to the liblsoda library, developed by Simon Frost (@sdwfrost), thereby providing a way to use the LSODA algorithm from Linda Petzold and Alan Hindmarsh from Julia. Clang.jl has been used to write the library and Sundials.jl was a inspiring source.

Installation

To install this package, run the command add LSODA.

Simplified Functions

To solve an ODE, one can call the simplified solver:

function rhs!(t, x, ydot, data)
	ydot[1]=1.0E4 * x[2] * x[3] - .04E0 * x[1]
	ydot[3]=3.0E7 * x[2] * x[2]
	ydot[2]=-ydot[1] - ydot[3]
  nothing
end

y0 = [1.,0.,0.]
tspan = [0., 0.4]
res =  lsoda(rhs!, y0, tspan, reltol= 1e-4, abstol = Vector([1.e-6,1.e-10,1.e-6]))

To reproduce the test example from liblsoda, on can use:

lsoda_0(rhs!, y0, tspan, reltol= 1e-4, abstol = Vector([1.e-6,1.e-10,1.e-6]))

This should give the following.

at t =   4.0000e-01 y=   9.851712e-01   3.386380e-05   1.479493e-02
at t =   4.0000e+00 y=   9.055333e-01   2.240655e-05   9.444430e-02
at t =   4.0000e+01 y=   7.158403e-01   9.186334e-06   2.841505e-01
at t =   4.0000e+02 y=   4.505250e-01   3.222964e-06   5.494717e-01
at t =   4.0000e+03 y=   1.831976e-01   8.941774e-07   8.168016e-01
at t =   4.0000e+04 y=   3.898729e-02   1.621940e-07   9.610125e-01
at t =   4.0000e+05 y=   4.936362e-03   1.984221e-08   9.950636e-01
at t =   4.0000e+06 y=   5.161832e-04   2.065786e-09   9.994838e-01
at t =   4.0000e+07 y=   5.179811e-05   2.072030e-10   9.999482e-01
at t =   4.0000e+08 y=   5.283524e-06   2.113420e-11   9.999947e-01
at t =   4.0000e+09 y=   4.658945e-07   1.863579e-12   9.999995e-01
at t =   4.0000e+10 y=   1.423392e-08   5.693574e-14   1.000000e+00

JuliaDiffEq Common Interface

The functionality of LSODA.jl can be accessed through the JuliaDiffEq common interface. To do this, you build a problem object for like:

using LSODA, DiffEqBase
function rhs!(du, u, p, t)
    du[1]=1.0E4 * u[2] * u[3] - .04E0 * u[1]
    du[3]=3.0E7 * u[2] * u[2]
    du[2]=-du[1] - du[3]
  nothing
end

y0 = [1.,0.,0.]
tspan = (0., 0.4)
prob = ODEProblem(rhs!,y0,tspan)

This problem is solved by LSODA by using the lsoda() algorithm in the common solve command as follows:

sol = solve(prob,lsoda())

Many keyword arguments can be used to control the solver, its tolerances, and its output formats. For more information, please see the DifferentialEquations.jl documentation.

Citing

If using the algorithm through the DifferentialEquations.jl common interface, please cite:

@article{rackauckas2017differentialequations,
  title={Differentialequations. jl--a performant and feature-rich ecosystem for solving differential equations in julia},
  author={Rackauckas, Christopher and Nie, Qing},
  journal={Journal of Open Research Software},
  volume={5},
  number={1},
  year={2017},
  publisher={Ubiquity Press}
}

For the original algorithm, please cite:

• Alan Hindmarsh, ODEPACK, a Systematized Collection of ODE Solvers, in Scientific Computing, edited by Robert Stepleman, Elsevier, 1983, ISBN13: 978-0444866073, LC: Q172.S35.

• K Radhakrishnan, Alan Hindmarsh, Description and Use of LSODE, the Livermore Solver for Ordinary Differential Equations, Technical report UCRL-ID-113855, Lawrence Livermore National Laboratory, December 1993.

• Linda Petzold, Automatic Selection of Methods for Solving Stiff and Nonstiff Systems of Ordinary Differential Equations, SIAM J. Sci. and Stat. Comput., 4(1), 136–148.