[SciPy-dev] ode - generic interface to ODEs solvers

[P.Peterson]

Hi!

I have commited a new module, ode, to the integrate package. It provides
a generic interface to various numerical integrators of ODE systems.
Currently it has support for vode routine. Here is a stiff sample problem
solved with using ode:

    from scipy.integrate import ode
    def f(t,y):
        ydot0 = -0.04*y[0] + 1e4*y[1]*y[2]
        ydot2 = 3e7*y[1]*y[1]
        ydot1 = -ydot0-ydot2
        return [ydot0,ydot1,ydot2]
    def jac(t,y):
        jc = [[-0.04,1e4*y[2]          ,1e4*y[1]],
              [0.04 ,-1e4*y[2]-6e7*y[1],-1e4*y[1]],
              [0.0    ,6e7*y[1]           ,0.0]]
        return jc
    r = ode(f,jac).set_integrator('vode',
                                  rtol=1e-4,
                                  atol=[1e-8,1e-14,1e-6],
                                  method='bdf',
                                  )
    r.set_initial_value([1,0,0])
    print 'At t=%s  y=%s'%(r.t,r.y)
    tout = 0.4
    for i in range(12):
        r.integrate(tout)
        print 'At t=%s  y=%s'%(r.t,r.y)
        tout *= 10

The output is:
At t=0.0  y=[ 1.  0.  0.]
At t=0.4  y=[  9.85172114e-01   3.38639538e-05   1.47940224e-02]
At t=4.0  y=[  9.05518679e-01   2.24047569e-05   9.44589164e-02]
At t=40.0  y=[  7.15827066e-01   9.18553466e-06   2.84163748e-01]
At t=400.0  y=[  4.50518664e-01   3.22290138e-06   5.49478113e-01]
At t=4000.0  y=[  1.83202258e-01   8.94237125e-07   8.16796848e-01]
At t=40000.0  y=[  3.89833604e-02   1.62176759e-07   9.61016477e-01]
At t=400000.0  y=[  4.93828907e-03   1.98499999e-08   9.95061691e-01]
At t=4000000.0  y=[  5.16818452e-04   2.06832993e-09   9.99483179e-01]
At t=40000000.0  y=[  5.20295745e-05   2.08128996e-10   9.99947970e-01]
At t=400000000.0  y=[  5.19687803e-06   2.07876187e-11   9.99994803e-01]
At t=4000000000.0  y=[  5.20558016e-07   2.08223292e-12   9.99999479e-01]
At t=40000000000.0  y=[  4.22779340e-08   1.69111761e-13   9.99999958e-01]



Another example, more realistic and useful in the usage, is in ode.py (see
the function test1()).

You may ask for the reasons why I wrote this module as the integrate
package already provides odeint function for the same task.

The main reason is that numerical intergration of ODE systems is a complex
task that can be very sensitive to the problem in hand. For example,
completely different methods are efficent for stiff and non-stiff
problems, respectively. There are many algorithms available that are very
efficient for certain problems while they can perform very poorly if
applied to wrong problems. Therefore, it is crucial that one can easily
switch algorithms for different problems and also to test which ones are
most efficient. One may even need to change the algorithm during
the integration. ode module tackles all these issues efficiently and at
the same time being rather user-friendly.

Currently the only supported routine vode is a Fortran program. There is
also a C version of vode available (called cvode) that includes many
improvements to vode, including Krylov solver (the former vodpk). So, it
is the next routine in my task list to get supported.

Comments are welcome as always.
Regards,
	Pearu