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 index1 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=1e3
/abstol=1e5
set the relative/absolute local error tolerances 
initstep=1e4
/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 equationF(t,y,a*y+b)=0
via Newton's method, which requires a jacobian of the formdF/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 asjacobian=(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 6th 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 choosey
andy'
such thatF(t,y,y')=0
initially. DASSL tries to guess the initial value ofy'
, 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 aserr=norm(ycy0)
, and accepting the step whenerr<1
. The local error tolerancesreltol
andabstol
are hidden in the definition ofdassl_norm(v, wt)=norm(v./wt)/sqrt(length(v))
, where weightswt
are defined bydassl_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