Solves stiff differential algebraic equations (DAE) using variable stepsize backwards finite difference formula (BDF) in the SciML scientific machine learning organization
This is an implementation of DASSL algorithm for solving algebraic differential equations. To install a stable version run
Pkg.add("DASSL")
This package is compatible with the JuliaDiffEq common solver interface which is documented in the DifferentialEquations.jl documentation. Following the DAE Tutorial, one can use dassl()
as follows:
using DASSL
u0 = [1.0, 0, 0]
du0 = [-0.04, 0.04, 0.0]
tspan = (0.0,100000.0)
function resrob(r,yp,y,p,t)
r[1] = -0.04*y[1] + 1.0e4*y[2]*y[3]
r[2] = -r[1] - 3.0e7*y[2]*y[2] - yp[2]
r[1] -= yp[1]
r[3] = y[1] + y[2] + y[3] - 1.0
end
prob = DAEProblem(resrob,du0,u0,tspan)
sol = solve(prob, dassl())
For more details on using this interface, see the ODE tutorial.
To solve a scalar equation y'(t)+y(t)=0
with initial data y(0)=0.0
up to time t=10.0
run the following code
using DASSL
F(t,y,dy) = dy+y # the equation solved is F(t,y,dy)=0
y0 = 1.0 # the initial value
tspan = [0.0,10.0] # time span over which we integrate
(tn,yn) = dasslSolve(F,y0,tspan) # returns (tn,yn)
You can also change the relative error tolerance rtol
, absolute
error tolerance atol
as well as initial step size h0
as follows
(tn,yn) = dasslSolve(F,y0,tspan)
To test the convergence and execution time for index-1 problem run
convergence.jl
from the test
directory.
Naturally, DASSL.jl also supports multiple equations. For example the pendulum equation
u'-v=0
v'+sin(u)=0
with initial data u(0)=0.0
and v(0)=1.0
can be solved by defining
the following residual function
function F(t,y,dy)
[
dy[1]-y[2], # y[1]=u, y[2]=v
dy[2]+sin(y[1]) # dy[1]=u', dy[2]=v'
]
end
The initial data shoud now be set as a vector
y0 = [0.0,1.0] # y0=[u(0),v(0)]
The solution can be computed by calling
tspan = [0.0,10.0]
(tn,yn) = dasslSolve(F,y0,tspan)
Apart from producing the times tn
and values yn
, dasslSolve also
produces the derivatives dyn
(as the byproduct of BDF
algorithm), e.g.
(tn,yn,dyn) = dasslSolve(F,y0,tspan)
The decision to produce these values is that it is not entirely
trivial to compute y'
from F(t,y,y')=0
when t
and y
are given.
DASSL supports a number of keyword arguments, the names of most of them are compatible with the namse used in ODE package.
reltol=1e-3
/abstol=1e-5
set the relative/absolute local error tolerances
initstep=1e-4
/minstep=0
/maxstep=Inf
set the
initial/minimal/maximal step sizes (when step size drops below
minimum the integration stops)
jacobian
The most expensive step during the integration is solving
the nonlinear equation F(t,y,a*y+b)=0
via Newton's method, which
requires a jacobian of the form dF/dy+a*dF/dy'
. By default, the
solver approximates this Jacobian by a method of finite differences
but you can provide your own method as a function
(t,y,dy,a)->dF/dy+a*dF/dy'
. For the pendulum equation we would
define jacobian as
jacobian=(t,y,dy,a)->[[a,cos(y[1])] [-1,a]]
maxorder=6
Apart from selecting the current step size DASSL method
can also dynamically change the order of BDF method used. BDF is
stable up to 6-th order, which is the defaul upper limit but for
some systems of equations it may make more sense to use lower
orders.
dy0=zero(y)
When solving differential algebraic equations it is
important to start with consistent initial conditions, i.e. to
choose y
and y'
such that F(t,y,y')=0
initially. DASSL tries
to guess the initial value of y'
, but if it fails you can set your
own initial condtions for the derivative.
norm=dassl_norm
/weights=dassl_weights
DASSL computes the error
roughly as err=norm(yc-y0)
, and accepting the step when
err<1
. The local error tolerances reltol
and abstol
are
hidden in the definition of dassl_norm(v, wt)=norm(v./wt)/sqrt(length(v))
, where weights wt
are defined by
dassl_weights(y,reltol,abstol)=reltol*abs(y).+abstol
. You can
supply your own weights and norms when they are more appropriate for
the problem at hand.
factorize_jacobian=true
is a Boolean option which forces the
factorization of Jacobian before storing it. It dramatically
increases performance for large systems, but may decrease the
computation speed for small systems.
DASSL.jl supports an iterative version of solver (implemented via
coroutines, so debugging might be a little tricky) via
dasslIterator
. In the following example the dasslIterator
is used
to stop the integration when the solution y
drops below 0.1
F(t,y,dy)=dy+y
# iterator version of dassl solver
for (t,y,dy) in dasslIterator(F,1.0,0.0)
if y < 0.1
@show (t,y,dy)
break
end
end