sciml / dassl.jl

Solves stiff differential algebraic equations (DAE) using variable stepsize backwards finite difference formula (BDF) in the SciML scientific machine learning organization

DASSL.jl

This is an implementation of DASSL algorithm for solving algebraic differential equations. To install a stable version run

Pkg.add("DASSL")

Common Interface Example

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.

Examples

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)

Output

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.

Keyword arguments

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.

Iterator version

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

dassl.jl's Issues

Move the interpolating polynomials to a separate package

Remove the definitions of InterpolateAt, InterpolateDerivativeAt etc. from DASSL.jl and move them to the InterPoly.jl package (or other suitable Julia package).

Improve efficiency

File src/DASSL.jl line 686 has the function:

# generate a function that computes approximate jacobian using forward
# finite differences
function numerical_jacobian(F,reltol,abstol,weights)
    function numjac(t, y, dy, a)
        ep      = eps(one(t))   # this is the machine epsilon
        h       = 1/a           # h ~ 1/a
        wt      = weights(y,reltol,abstol)
        # delta for approximation of jacobian.  I removed the
        # sign(h_next*dy0) from the definition of delta because it was
        # causing trouble when dy0==0 (which happens for ord==1)
        edelta  = diagm(0=>max.(abs.(y),abs.(h*dy),wt)*sqrt(ep))

        b=dy-a*y
        f(y1) = F(t,y1,a*y1+b)

        n   = length(y)
        jac = Array{eltype(y)}(undef,n,n)
        for i=1:n
            jac[:,i]=(f(y+edelta[:,i])-f(y))/edelta[i,i]
        end

        jac
    end
end

This (central) function should be improved:

1. In the for loop f(y) is called in every iteration, although y does not change. This means that there are n-1 unnecessary calls of the right hand side function. The call of f(y) should be moved out of the for loop and the result stored in a vector.

2. The statement jac = Array{eltype(y)}(undef,n,n) generates the Jacobian matrix whenever this function is called. So, when this function is called 100 times, then memory for 100 Jacobians are allocated. This should be improved: At the start of the simulation, memory for one Jacobian should be allocated that is then updated within this function.

convergence tests fails

I tried to run the convergence test:

...v0.3/DASSL/test((HEAD detached at a17ab3f))  >> julia convergence.jl 
Warning: could not import Base.Text into Tk
INFO: Stepsize too small (h=3.637978807091713e-16 at t=0.0.
INFO: Stepsize too small (h=3.637978807091713e-16 at t=0.0.
ERROR: dimensions must match
 in - at array.jl:723
 in dasslTestConvergence at /home/mauro/.julia/v0.3/DASSL/test/convergence.jl:35
 in dasslDrawConvergence at /home/mauro/.julia/v0.3/DASSL/test/convergence.jl:59
 in include at ./boot.jl:245
 in include_from_node1 at loading.jl:128
 in process_options at ./client.jl:285
 in _start at ./client.jl:354
while loading /home/mauro/.julia/v0.3/DASSL/test/convergence.jl, in expression starting on line 114

[PackageEvaluator.jl] Your package DASSL may have a testing issue.

This issue is being filed by a script, but if you reply, I will see it.

PackageEvaluator.jl is a script that runs nightly. It attempts to load all Julia packages and run their test (if available) on both the stable version of Julia (0.2) and the nightly build of the unstable version (0.3).

The results of this script are used to generate a package listing enhanced with testing results.

The status of this package, DASSL, on...

• Julia 0.2 is 'Tests fail, but package loads.'
• Julia 0.3 is 'Tests fail, but package loads.'

'No tests, but package loads.' can be due to their being no tests (you should write some if you can!) but can also be due to PackageEvaluator not being able to find your tests. Consider adding a test/runtests.jl file.

'Package doesn't load.' is the worst-case scenario. Sometimes this arises because your package doesn't have BinDeps support, or needs something that can't be installed with BinDeps. If this is the case for your package, please file an issue and an exception can be made so your package will not be tested.

This automatically filed issue is a one-off message. Starting soon, issues will only be filed when the testing status of your package changes in a negative direction (gets worse). If you'd like to opt-out of these status-change messages, reply to this message.

allow any callable object, not just ::Function arguments

On the mailing list, someone was trying to call DASSL from Python via pyjulia. This should work, but the Python function gets passed as a (callable) ::PyObject argument, not a ::Function.

It would be better to leave the arguments untyped rather than specify ::Function, in order that any callable object can be passed. Note that there is no performance advantage to declaring argument types; the main reason to declare arguments is to avoid method ambiguities.

Example is failed

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.04y[1] + 1.0e4y[2]y[3]
r[2] = -r[1] - 3.0e7y[2]*y[2] - yp[2]
r[1] -= yp[1]
r[3] = y[1] + y[2] + y[3] - 1.0
end

prob = DAEProblem(resrob,u0,du0,tspan)
sol = solve(prob, dassl())

ERROR: LinearAlgebra.LAPACKException(1)
Stacktrace:
[1] check_channel_state at .\channels.jl:120 [inlined]
[2] take_unbuffered(::Channel{Any}) at .\channels.jl:318
[3] take! at .\channels.jl:306 [inlined]
[4] iterate(::Channel{Any}, ::Nothing) at .\channels.jl:386
[5] iterate at .\channels.jl:385 [inlined]
[6] #dasslSolve#16(::Array{Float64,1}, ::Base.Iterators.Pairs{Symbol,Real,NTuple{7,Symbol},NamedTuple{(:abstol, :reltol, :maxstep, :minstep, :initstep, :maxorder, :factorize_jacobian),Tuple{Float64,Float64,Float64,Float64,Float64,Int64,Bool}}}, ::Function, ::Function, ::Array{Float64,1}, ::Array{Float64,1}) at C:\Users\hzgzh.julia\packages\DASSL\ZVHVR\src\DASSL.jl:163
[7] (::getfield(DASSL, Symbol("#kw##dasslSolve")))(::NamedTuple{(:abstol, :reltol, :maxstep, :minstep, :initstep, :maxorder, :factorize_jacobian),Tuple{Float64,Float64,Float64,Float64,Float64,Int64,Bool}}, ::typeof(dasslSolve), ::Function, ::Array{Float64,1}, ::Array{Float64,1}) at .\none:0
[8] #solve#2(::Bool, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Nothing, ::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, ::typeof(solve), ::DAEProblem{Array{Float64,1},Array{Float64,1},Tuple{Float64,Float64},true,Nothing,DAEFunction{true,typeof(resrob),Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,Nothing}, ::dassl) at C:\Users\hzgzh.julia\packages\DASSL\ZVHVR\src\common.jl:36
[9] solve(::DAEProblem{Array{Float64,1},Array{Float64,1},Tuple{Float64,Float64},true,Nothing,DAEFunction{true,typeof(resrob),Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,Nothing}, ::dassl)
at C:\Users\hzgzh.julia\packages\DASSL\ZVHVR\src\common.jl:15
[10] top-level scope at none:0

Travis Clang?

Is there a reason why you have the Travis tests using the Clang compiler instead of the normal Travis setup? You don't have any binary dependencies so it seems very unnecessary. Also, release is kept as v0.4 now, should set the environment separately. I could put this to a standard Travis script if you want, but I just want to double check that you aren't using the compiler for something special.

Fix MAXORDER=3

If MAXORDER constant is changed to 3 the BoundsError() shows up.

van der Pol equation

Solving the van der Pol (VDP) equation with stringent error tolerances (reltol=abstol<0.2e-8) leads to a "step size too small" error. However, the tests on http://www.dm.uniba.it/~testset/problems/vdpol.php#stats show that DASSL works for VDP up to 10 significant digits, with DASSL.jl I only get 7 significant digits.

Test case:

# VDPOL with DASSL:  works
using DASSL
using ODE

# ode23s works at least as low as rtol=atol=1e-12 producing a relative error of 2e-9
# DASSL errors with  rtol=atol=1e-0 with too small step size
atol = .2e-8
rtol = .2e-8

dof=2
mu = 1e3
eps = 1/mu^2 # rescaled parameter
function fn(t,y)
    # the ode function dydt = f(t,y)
    out = zeros(Float64,dof)
    out[1] = y[2]
    prod = 
    out[2] = ( (1.0 - y[1]^2)*y[2] - y[1]) / eps
    return out
end
function jac(t,y)
    # the Jacobian of f
    #
    # Set to nothing if it is not known
    dfdy = zeros(Float64,dof,dof)
    dfdy[1,1] = 0.0
    dfdy[1,2] = 1.0
    dfdy[2,1] = (-2.0 * y[1]*y[2]-1.0)/eps
    dfdy[2,2] = (1.0 - y[1]^2)/eps
    return dfdy
end

# implicit eqs:
fni(t,y,dy)  = fn(t,y)-dy
Fy(t,y,dy) = jac(t,y)
Fdy(t,y,dy) = -eye(dof)

y0 = [2., 0]
tspan = [0., 2] # integration interval
(tn,yn_nojac,dyn)=DASSL.dasslSolve(fni, y0, tspan, reltol=atol, abstol=rtol)
(tn,yn,dyn)=DASSL.dasslSolve(fni, y0, tspan, Fy = Fy, Fdy = Fdy, reltol=atol, abstol=rtol)
tout,yout = ode23s(fn, y0, tspan, jacobian=jac, reltol=atol, abstol=rtol)

# reference:
refsol = [0.1706167732170483e1, -0.8928097010247975e0] # reference solution at tspan[2]

relerror(sol, refsol) = norm(abs((sol-refsol)./refsol),Inf)
calc_error_scd(sol, refsol) = -log10(norm(relerror(sol,refsol), Inf))

@show relerror(yn_nojac[end],refsol)
@show relerror(yn[end],refsol)
@show relerror(yout[end],refsol)

@show calc_error_scd(yn_nojac[end],refsol)
@show calc_error_scd(yn[end],refsol)
@show calc_error_scd(yout[end],refsol)

produces:

>> julia vdpol_dassl.jl
relerror(yn_nojac[end],refsol) => 3.714141830332293e-8
relerror(yn[end],refsol) => 4.626878704902238e-8
relerror(yout[end],refsol) => 3.559129266522373e-7
calc_error_scd(yn_nojac[end],refsol) => 7.430141516063499
calc_error_scd(yn[end],refsol) => 7.334711885518931
calc_error_scd(yout[end],refsol) => 6.4486562382625605

Common API

Since the common API has pretty much settled, we should work towards getting it mapped onto DASSL. It shouldn't be difficult at all to get a pretty basic setup together.

Move to JuliaDiffEq?

I would like to move this to JuliaDiffEq so it could be updated to v0.5 testing/compatibility and be added as a backend to AlgebraicDiffEq.jl. What do you think about that?

Pass in exact Jacobian

It would be nice if a function could be passed in for the Jacobian instead of hard-coded finite differences.

Enable AppVeyor Testing

I think AppVeyor testing should be enabled on all of the JuliaDiffEq repositories in addition to Travis. Windows is a large set of users (myself included) and AppVeyor many times will catch Windows-specific issues (usually dealing with the way integers can be different, but also sometimes with how roundoff is handled or obscure facts that are hard to catch by yourself). This code is simple enough that I don't expect it to be different than the Travis results, but it's a very simple no-cost double-check and so I think it should be enabled.

Register 2.4.0

@JuliaRegistrator register()

TagBot trigger issue

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!