Stiff ODE-solver

Question

I need an ODE-solver for a stiff problem similar to MATLAB ode15s.

For my problem I need to check how many steps (calculations) is needed for different initial values and compare this to my own ODE-solver.

I tried using

solver = scipy.integrate.ode(f)
solver.set_integrator('vode', method='bdf', order=15, nsteps=3000)
solver.set_initial_value(u0, t0)

And then integrating with:

i = 0
while solver.successful() and solver.t<tf:
    solver.integrate(tf, step=True)
    i += 1
print(i)

Where tf is the end of my time interval.

The function used is defined as:

def func(self, t, u):
    u1 = u[1]
    u2 = mu * (1-numpy.dot(u[0], u[0]))*u[1] - u[0]
    return numpy.array([u1, u2])

Which with the initial value u0 = [ 2, 0] is a stiff problem.

This means that the number of steps should not depend on my constant mu.

But it does.

I think the odeint-method can solve this as a stiff problem - but then I have to send in the whole t-vector and therefore need to set the amount of steps that is done and this ruins the point of my assignment.

Is there anyway to use odeint with adaptive stepsize between two t0 and tf?

Or can you see anything I miss in the use of the vode-integrator?

Answer 1

I'm seeing something similar; with the 'vode' solver, changing methods between 'adams' and 'bdf' doesn't change the number of steps by very much. (By the way, there is no point in using order=15; the maximum order of the 'bdf' method of the 'vode' solver is 5 (and the maximum order of the 'adams' solver is 12). If you leave the argument out, it should use the maximum by default.)

odeint is a wrapper of LSODA. ode also provides a wrapper of LSODA: change 'vode' to 'lsoda'. Unfortunately the 'lsoda' solver ignores the step=True argument of the integrate method.

The 'lsoda' solver does much better than 'vode' with method='bdf'. You can get an upper bound on the number of steps that were used by initializing tvals = [], and in func, do tvals.append(t). When the solver completes, set tvals = np.unique(tvals). The length of tvals tells you the number of time values at which your function was evaluated. This is not exactly what you want, but it does show a huge difference between using the 'lsoda' solver and the 'vode' solver with method 'bdf'. The number of steps used by the 'lsoda' solver is on the same order as you quoted for matlab in your comment. (I used mu=10000, tf = 10.)

Update: It turns out that, at least for a stiff problem, it make a huge difference for the 'vode' solver if you provide a function to compute the Jacobian matrix.

The script below runs the 'vode' solver with both methods, and it runs the 'lsoda' solver. In each case, it runs the solver with and without the Jacobian function. Here's the output it generates:

vode   adams    jac=None  len(tvals) = 517992
vode   adams    jac=jac   len(tvals) = 195
vode   bdf      jac=None  len(tvals) = 516284
vode   bdf      jac=jac   len(tvals) = 55
lsoda           jac=None  len(tvals) = 49
lsoda           jac=jac   len(tvals) = 49

The script:

from __future__ import print_function

import numpy as np
from scipy.integrate import ode


def func(t, u, mu):
    tvals.append(t)
    u1 = u[1]
    u2 = mu*(1 - u[0]*u[0])*u[1] - u[0]
    return np.array([u1, u2])


def jac(t, u, mu):
    j = np.empty((2, 2))
    j[0, 0] = 0.0
    j[0, 1] = 1.0
    j[1, 0] = -mu*2*u[0]*u[1] - 1
    j[1, 1] = mu*(1 - u[0]*u[0])
    return j


mu = 10000.0
u0 = [2, 0]
t0 = 0.0
tf = 10

for name, kwargs in [('vode', dict(method='adams')),
                     ('vode', dict(method='bdf')),
                     ('lsoda', {})]:
    for j in [None, jac]:
        solver = ode(func, jac=j)
        solver.set_integrator(name, atol=1e-8, rtol=1e-6, **kwargs)
        solver.set_f_params(mu)
        solver.set_jac_params(mu)
        solver.set_initial_value(u0, t0)

        tvals = []
        i = 0
        while solver.successful() and solver.t < tf:
            solver.integrate(tf, step=True)
            i += 1

        print("%-6s %-8s jac=%-5s " %
              (name, kwargs.get('method', ''), j.func_name if j else None),
              end='')

        tvals = np.unique(tvals)
        print("len(tvals) =", len(tvals))