Demonstrate a bug in the handling of a banded Jacobian when using the 'vode' solver method of scipy.integrate.ode.

ode_banded_jac_example.py

import numpy as np
from scipy.integrate import ode
import matplotlib.pyplot as plt


def func(t, y, c):
    return c.dot(y)

def jac(t, y, c):
    return c

def banded_jac(t, y, c):
    # Banded Jacobian of the linear system `func` with ml=1, mu=0.
    jac = np.array([np.diag(c),
                    np.r_[np.diag(c, -1), 0]])
    return jac

def banded_jac_padded(t, y, c):
    # Same as banded_jac, but with an extra row of zeros.
    # This padding shoud not be necessary!
    jac = np.array([np.diag(c),
                    np.r_[np.diag(c, -1), 0],
                    np.zeros(c.shape[0])])
    return jac

c = np.array([[ -0.5,    0,    0],
              [  0.5, -0.5,    0],
              [  0,      0,  -0.5]])

lband = 1
uband = 0

t0 = 0
y0 = np.arange(1, c.shape[0] + 1)

# Using banded_jac with the 'vode' method, the solver raises an exception.
# Change the method to 'lsoda', and it works.
# Or, change the jacobian to banded_jac_padded, and it works with 'vode'.
# The 'vode' method with banded_jac also works if c is transposed (so
# lband=0 and uband=1).
r = ode(func, banded_jac)
#r = ode(func, banded_jac_padded)
r.set_integrator('vode',
                 with_jacobian=True,
                 method='bdf',
                 lband=lband, uband=uband,
                 )
r.set_initial_value(y0, t0)
r.set_f_params(c)
r.set_jac_params(c)

t1 = 10
dt = 0.05
t = [t0]
y = [y0]
while r.successful() and r.t < t1:
    r.integrate(r.t + dt)
    t.append(r.t)
    y.append(r.y)

t = np.array(t)
y = np.array(y)

plt.plot(t, y)
plt.xlim(t0, t1)
plt.grid(True)
plt.show()