RFC: solver ode23s for stiff problems

[acroy]

The implementation is based on the "modified Rosenbrock triple" described in

L.F. Shampine and M.W. Reichelt: "The MATLAB ODE Suite," SIAM Journal on Scientific Computing, Vol. 18, 1997, pp. 1–22

which is also the basis for MATLAB's ode23s. In addition to the usual keywords, the solver accepts a keyword jacobian, which can be used to supply a function giving the Jacobian of the RHS. Otherwise a (very crude) finite-difference method is used (pulled out from oderosenbrock).

A simple example is:

# Example: van der Pol oscillator
import ODE: ode23s

function vdP(t,y,mu)
    dy = zeros(y)

    dy[1] = y[2]
    dy[2] = -y[1] + mu*(1 - y[1]^2)*y[2]

    return dy
end

function jac_vdP(t,y,mu)
    dFdy = zeros(length(y), length(y))

    dFdy[1,1] = 0.
    dFdy[1,2] = 1.

    dFdy[2,1] = -1. - mu*2*y[1]*y[2]
    dFdy[2,2] = mu*(1 - y[1]^2)

    return dFdy
end

tout, yout = ode23s((t,y)->vdP(t,y,1000), [2.,0.], [0., 3000.], jacobian=(t,y)->jac_vdP(t,y,1000), reltol=1e-3, abstol=1e-6)

Yielding for the first component of yout

This PR should also "fix" #1.

[coveralls]

Coverage increased (+0.82%) when pulling 3b20e24 on acroy:ode23s into 913a3a5 on JuliaLang:master.

[tomasaschan]

Nice work!

Are stiff solvers generalizable in the same way we've generalized oderkf to do the work of both ode23 and ode45 for various sets of coefficients?`If so, would it be possible to generalize this code to allow for later addition of other sets?

[acroy]

In principle this should be possible. The Rosenbrock methods (see e.g. generalized rosenbrock formulation and of course Hairer & Wanner) seem to be very similar to RK methods. In Solving Ordinary Differential Equations, II Hairer & Wanner even write that one can basically take DOPRI5 and add a few lines of code for the Jacobian and solving linear systems and et voilà ... (well almost)

In the long run we probably want to have something more general (and not only Rosenbrock based), but for now ode23s should be Ok.

[coveralls]

Changes Unknown when pulling f25ef27 on acroy:ode23s into * on JuliaLang:master*.

[ViralBShah]

Isn't cvode good enough for stiff problems?

[acroy]

Sure. And right now it is probably better/faster/more reliable for almost any problem. But it depends on the Sundials library, while in ODE.jl we are trying to collect pure Julia implementations. This allows it for example to use ode23s et al to be used with user-defined types.

[ViralBShah]

Ok. Thanks - I thought as much, but good to have the confirmation.

[acroy]

Anyone minds if we merge this (after I cleaned up my commits)?

[ViralBShah]

I guess that should be ok, but it would be nice to get a review.

Cc: @vtjnash @stevengj

[ViralBShah]

I'd say go ahead with the cleanup and merging. I do know enough about ode23s to review the code, but we can fix issues as they come up.

[coveralls]

Coverage increased (+0.82%) when pulling 0b6401b on acroy:ode23s into 568c213 on JuliaLang:master.

[acroy]

I agree. And the likelihood of finding an issue is much bigger if the method is on master.