van der Pol equation

[mauro3]

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

[pwl]

I was deducing minimal allowed stepsize from eps(T) where T is the type of the time variable. When eps is changed to eps(t), where t is the current time your example works fine (the epsilon for t=0.0 is ~e-324...).

EDIT: After thinking about it, the van der Pol equation should be time independent, so starting from t=1.0 brings the bug back...

[mauro3]

How does the time step evolve? Your fix above suggests that the time step is >eps(1) at t~1 even though the solution repeats after t~1.6. Maybe it's a problem with the initial steps? Anyway, using eps(t) instead of eps(T) is probably good irrespective.

[pwl]

This fixes the two errors I made when implementing DASSL.

1. I was not using the initial value for the derivative properly (the dy0 keyword argument).
2. When dy0 was not given I was not obtaining it from solving F(t,y0,dy0)=0 for dy0.

These two errors combined were forcing the stepsize to go way below the sensible values thus causing the premature termination of integration due to h<hmin.

In this fix I also rewrote the way error estimates are generated and handled (now it's clearer and simpler) and I slightly changed the order/timestep selection strategy.

Now it is possible to run the solver (with or without Jacobian) with the local error tolerances reduced to 1e-14 and still have h>1e-16 (although for this particular example it leads to larger global error due to roundoff errors).

[mauro3]

thanks