# Integration (`scipy.integrate`)

The `scipy.integrate` sub-package provides several integration techniques including an ordinary differential equation integrator. An overview of the module is provided by the help command:

    >>> help(integrate)
     Methods for Integrating Functions given function object.

       quad          -- General purpose integration.
       dblquad       -- General purpose double integration.
       tplquad       -- General purpose triple integration.
       fixed_quad    -- Integrate func(x) using Gaussian quadrature of order n.
       quadrature    -- Integrate with given tolerance using Gaussian quadrature.
       romberg       -- Integrate func using Romberg integration.

     Methods for Integrating Functions given fixed samples.

       trapezoid            -- Use trapezoidal rule to compute integral.
       cumulative_trapezoid -- Use trapezoidal rule to cumulatively compute integral.
       simpson              -- Use Simpson's rule to compute integral from samples.
       romb                 -- Use Romberg Integration to compute integral from
                            -- (2**k + 1) evenly-spaced samples.

       See the special module's orthogonal polynomials (special) for Gaussian
          quadrature roots and weights for other weighting factors and regions.

     Interface to numerical integrators of ODE systems.

       odeint        -- General integration of ordinary differential equations.
       ode           -- Integrate ODE using VODE and ZVODE routines.

## General integration (`quad`)

The function `quad` is provided to integrate a function of one variable between two points. The points can be \(\pm\infty\) (\(\pm\) `inf`) to indicate infinite limits. For example, suppose you wish to integrate a bessel function `jv(2.5, x)` along the interval \([0, 4.5].\)

\[I=\int_{0}^{4.5}J_{2.5}\left(x\right)\, dx.\]

This could be computed using `quad`:

    >>> import scipy.integrate as integrate
    >>> import scipy.special as special
    >>> result = integrate.quad(lambda x: special.jv(2.5,x), 0, 4.5)
    >>> result
    (1.1178179380783249, 7.8663172481899801e-09)

    >>> from numpy import sqrt, sin, cos, pi
    >>> I = sqrt(2/pi)*(18.0/27*sqrt(2)*cos(4.5) - 4.0/27*sqrt(2)*sin(4.5) +
    ...                 sqrt(2*pi) * special.fresnel(3/sqrt(pi))[0])
    >>> I
    1.117817938088701

    >>> print(abs(result[0]-I))
    1.03761443881e-11

The first argument to quad is a “callable” Python object (i.e., a function, method, or class instance). Notice the use of a lambda- function in this case as the argument. The next two arguments are the limits of integration. The return value is a tuple, with the first element holding the estimated value of the integral and the second element holding an upper bound on the error. Notice, that in this case, the true value of this integral is

\[I=\sqrt{\frac{2}{\pi}}\left(\frac{18}{27}\sqrt{2}\cos\left(4.5\right)-\frac{4}{27}\sqrt{2}\sin\left(4.5\right)+\sqrt{2\pi}\textrm{Si}\left(\frac{3}{\sqrt{\pi}}\right)\right),\]

where

\[\textrm{Si}\left(x\right)=\int_{0}^{x}\sin\left(\frac{\pi}{2}t^{2}\right)\, dt.\]

is the Fresnel sine integral. Note that the numerically-computed integral is within \(1.04\times10^{-11}\) of the exact result — well below the reported error bound.

If the function to integrate takes additional parameters, they can be provided in the args argument. Suppose that the following integral shall be calculated:

\[I(a,b)=\int_{0}^{1} ax^2+b \, dx.\]

This integral can be evaluated by using the following code:

    >>> from scipy.integrate import quad
    >>> def integrand(x, a, b):
    ...     return a*x**2 + b
    ...
    >>> a = 2
    >>> b = 1
    >>> I = quad(integrand, 0, 1, args=(a,b))
    >>> I
    (1.6666666666666667, 1.8503717077085944e-14)

Infinite inputs are also allowed in `quad` by using \(\pm\) `inf` as one of the arguments. For example, suppose that a numerical value for the exponential integral:

\[E_{n}\left(x\right)=\int_{1}^{\infty}\frac{e^{-xt}}{t^{n}}\, dt.\]

is desired (and the fact that this integral can be computed as `special.expn(n,x)` is forgotten). The functionality of the function `special.expn` can be replicated by defining a new function `vec_expint` based on the routine `quad`:

    >>> from scipy.integrate import quad
    >>> import numpy as np
    >>> def integrand(t, n, x):
    ...     return np.exp(-x*t) / t**n
    ...

    >>> def expint(n, x):
    ...     return quad(integrand, 1, np.inf, args=(n, x))[0]
    ...

    >>> vec_expint = np.vectorize(expint)

    >>> vec_expint(3, np.arange(1.0, 4.0, 0.5))
    array([ 0.1097,  0.0567,  0.0301,  0.0163,  0.0089,  0.0049])
    >>> import scipy.special as special
    >>> special.expn(3, np.arange(1.0,4.0,0.5))
    array([ 0.1097,  0.0567,  0.0301,  0.0163,  0.0089,  0.0049])

The function which is integrated can even use the quad argument (though the error bound may underestimate the error due to possible numerical error in the integrand from the use of `quad` ). The integral in this case is

\[I_{n}=\int_{0}^{\infty}\int_{1}^{\infty}\frac{e^{-xt}}{t^{n}}\, dt\, dx=\frac{1}{n}.\]

    >>> result = quad(lambda x: expint(3, x), 0, np.inf)
    >>> print(result)
    (0.33333333324560266, 2.8548934485373678e-09)

    >>> I3 = 1.0/3.0
    >>> print(I3)
    0.333333333333

    >>> print(I3 - result[0])
    8.77306560731e-11

This last example shows that multiple integration can be handled using repeated calls to `quad`.

Warning

Numerical integration algorithms sample the integrand at a finite number of points. Consequently, they cannot guarantee accurate results (or accuracy estimates) for arbitrary integrands and limits of integration. Consider the Gaussian integral, for example:

    >>> def gaussian(x):
    ...     return np.exp(-x**2)
    >>> res = integrate.quad(gaussian, -np.inf, np.inf)
    >>> res
    (1.7724538509055159, 1.4202636756659625e-08)
    >>> np.allclose(res[0], np.sqrt(np.pi))  # compare against theoretical result
    True

Since the integrand is nearly zero except near the origin, we would expect large but finite limits of integration to yield the same result. However:

    >>> integrate.quad(gaussian, -10000, 10000)
    (1.975190562208035e-203, 0.0)

This happens because the adaptive quadrature routine implemented in `quad`, while working as designed, does not notice the small, important part of the function within such a large, finite interval. For best results, consider using integration limits that tightly surround the important part of the integrand.

    >>> integrate.quad(gaussian, -15, 15)
    (1.772453850905516, 8.476526631214648e-11)

Integrands with several important regions can be broken into pieces as necessary.

## General multiple integration (`dblquad`, `tplquad`, `nquad`)

The mechanics for double and triple integration have been wrapped up into the functions `dblquad` and `tplquad`. These functions take the function to integrate and four, or six arguments, respectively. The limits of all inner integrals need to be defined as functions.

An example of using double integration to compute several values of \(I_{n}\) is shown below:

    >>> from scipy.integrate import quad, dblquad
    >>> def I(n):
    ...     return dblquad(lambda t, x: np.exp(-x*t)/t**n, 0, np.inf, lambda x: 1, lambda x: np.inf)
    ...

    >>> print(I(4))
    (0.2500000000043577, 1.29830334693681e-08)
    >>> print(I(3))
    (0.33333333325010883, 1.3888461883425516e-08)
    >>> print(I(2))
    (0.4999999999985751, 1.3894083651858995e-08)

As example for non-constant limits consider the integral

\[I=\int_{y=0}^{1/2}\int_{x=0}^{1-2y} x y \, dx\, dy=\frac{1}{96}.\]

This integral can be evaluated using the expression below (Note the use of the non-constant lambda functions for the upper limit of the inner integral):

    >>> from scipy.integrate import dblquad
    >>> area = dblquad(lambda x, y: x*y, 0, 0.5, lambda x: 0, lambda x: 1-2*x)
    >>> area
    (0.010416666666666668, 1.1564823173178715e-16)

For n-fold integration, scipy provides the function `nquad`. The integration bounds are an iterable object: either a list of constant bounds, or a list of functions for the non-constant integration bounds. The order of integration (and therefore the bounds) is from the innermost integral to the outermost one.

The integral from above

\[I_{n}=\int_{0}^{\infty}\int_{1}^{\infty}\frac{e^{-xt}}{t^{n}}\, dt\, dx=\frac{1}{n}\]

can be calculated as

    >>> from scipy import integrate
    >>> N = 5
    >>> def f(t, x):
    ...    return np.exp(-x*t) / t**N
    ...
    >>> integrate.nquad(f, [[1, np.inf],[0, np.inf]])
    (0.20000000000002294, 1.2239614263187945e-08)

Note that the order of arguments for f must match the order of the integration bounds; i.e., the inner integral with respect to \(t\) is on the interval \([1, \infty]\) and the outer integral with respect to \(x\) is on the interval \([0, \infty]\).

Non-constant integration bounds can be treated in a similar manner; the example from above

\[I=\int_{y=0}^{1/2}\int_{x=0}^{1-2y} x y \, dx\, dy=\frac{1}{96}.\]

can be evaluated by means of

    >>> from scipy import integrate
    >>> def f(x, y):
    ...     return x*y
    ...
    >>> def bounds_y():
    ...     return [0, 0.5]
    ...
    >>> def bounds_x(y):
    ...     return [0, 1-2*y]
    ...
    >>> integrate.nquad(f, [bounds_x, bounds_y])
    (0.010416666666666668, 4.101620128472366e-16)

which is the same result as before.

## Gaussian quadrature

A few functions are also provided in order to perform simple Gaussian quadrature over a fixed interval. The first is `fixed_quad`, which performs fixed-order Gaussian quadrature. The second function is `quadrature`, which performs Gaussian quadrature of multiple orders until the difference in the integral estimate is beneath some tolerance supplied by the user. These functions both use the module `scipy.special.orthogonal`, which can calculate the roots and quadrature weights of a large variety of orthogonal polynomials (the polynomials themselves are available as special functions returning instances of the polynomial class — e.g., `special.legendre`).

## Romberg Integration

Romberg’s method [WPR] is another method for numerically evaluating an integral. See the help function for `romberg` for further details.

## Integrating using Samples

If the samples are equally-spaced and the number of samples available is \(2^{k}+1\) for some integer \(k\), then Romberg `romb` integration can be used to obtain high-precision estimates of the integral using the available samples. Romberg integration uses the trapezoid rule at step-sizes related by a power of two and then performs Richardson extrapolation on these estimates to approximate the integral with a higher degree of accuracy.

In case of arbitrary spaced samples, the two functions `trapezoid` and `simpson` are available. They are using Newton-Coates formulas of order 1 and 2 respectively to perform integration. The trapezoidal rule approximates the function as a straight line between adjacent points, while Simpson’s rule approximates the function between three adjacent points as a parabola.

For an odd number of samples that are equally spaced Simpson’s rule is exact if the function is a polynomial of order 3 or less. If the samples are not equally spaced, then the result is exact only if the function is a polynomial of order 2 or less.

    >>> import numpy as np
    >>> def f1(x):
    ...    return x**2
    ...
    >>> def f2(x):
    ...    return x**3
    ...
    >>> x = np.array([1,3,4])
    >>> y1 = f1(x)
    >>> from scipy import integrate
    >>> I1 = integrate.simpson(y1, x)
    >>> print(I1)
    21.0

This corresponds exactly to

\[\int_{1}^{4} x^2 \, dx = 21,\]

whereas integrating the second function

    >>> y2 = f2(x)
    >>> I2 = integrate.simpson(y2, x)
    >>> print(I2)
    61.5

does not correspond to

\[\int_{1}^{4} x^3 \, dx = 63.75\]

because the order of the polynomial in f2 is larger than two.

## Faster integration using low-level callback functions

A user desiring reduced integration times may pass a C function pointer through `scipy.LowLevelCallable` to `quad`, `dblquad`, `tplquad` or `nquad` and it will be integrated and return a result in Python. The performance increase here arises from two factors. The primary improvement is faster function evaluation, which is provided by compilation of the function itself. Additionally we have a speedup provided by the removal of function calls between C and Python in `quad`. This method may provide a speed improvements of ~2x for trivial functions such as sine but can produce a much more noticeable improvements (10x+) for more complex functions. This feature then, is geared towards a user with numerically intensive integrations willing to write a little C to reduce computation time significantly.

The approach can be used, for example, via `ctypes` in a few simple steps:

1.) Write an integrand function in C with the function signature `double f(int n, double *x, void *user_data)`, where `x` is an array containing the point the function f is evaluated at, and `user_data` to arbitrary additional data you want to provide.

    /* testlib.c */
    double f(int n, double *x, void *user_data) {
        double c = *(double *)user_data;
        return c + x[0] - x[1] * x[2]; /* corresponds to c + x - y * z */
    }

2.) Now compile this file to a shared/dynamic library (a quick search will help with this as it is OS-dependent). The user must link any math libraries, etc., used. On linux this looks like:

    $ gcc -shared -fPIC -o testlib.so testlib.c

The output library will be referred to as `testlib.so`, but it may have a different file extension. A library has now been created that can be loaded into Python with `ctypes`.

3.) Load shared library into Python using `ctypes` and set `restypes` and `argtypes` - this allows SciPy to interpret the function correctly:

    import os, ctypes
    from scipy import integrate, LowLevelCallable

    lib = ctypes.CDLL(os.path.abspath('testlib.so'))
    lib.f.restype = ctypes.c_double
    lib.f.argtypes = (ctypes.c_int, ctypes.POINTER(ctypes.c_double), ctypes.c_void_p)

    c = ctypes.c_double(1.0)
    user_data = ctypes.cast(ctypes.pointer(c), ctypes.c_void_p)

    func = LowLevelCallable(lib.f, user_data)

The last `void *user_data` in the function is optional and can be omitted (both in the C function and ctypes argtypes) if not needed. Note that the coordinates are passed in as an array of doubles rather than a separate argument.

4.) Now integrate the library function as normally, here using `nquad`:

    >>> integrate.nquad(func, [[0, 10], [-10, 0], [-1, 1]])
    (1200.0, 1.1102230246251565e-11)

The Python tuple is returned as expected in a reduced amount of time. All optional parameters can be used with this method including specifying singularities, infinite bounds, etc.

## Ordinary differential equations (`solve_ivp`)

Integrating a set of ordinary differential equations (ODEs) given initial conditions is another useful example. The function `solve_ivp` is available in SciPy for integrating a first-order vector differential equation:

\[\frac{d\mathbf{y}}{dt}=\mathbf{f}\left(\mathbf{y},t\right),\]

given initial conditions \(\mathbf{y}\left(0\right)=y_{0}\), where \(\mathbf{y}\) is a length \(N\) vector and \(\mathbf{f}\) is a mapping from \(\mathcal{R}^{N}\) to \(\mathcal{R}^{N}.\) A higher-order ordinary differential equation can always be reduced to a differential equation of this type by introducing intermediate derivatives into the \(\mathbf{y}\) vector.

For example, suppose it is desired to find the solution to the following second-order differential equation:

\[\frac{d^{2}w}{dz^{2}}-zw(z)=0\]

with initial conditions \(w\left(0\right)=\frac{1}{\sqrt[3]{3^{2}}\Gamma\left(\frac{2}{3}\right)}\) and \(\left.\frac{dw}{dz}\right|_{z=0}=-\frac{1}{\sqrt[3]{3}\Gamma\left(\frac{1}{3}\right)}.\) It is known that the solution to this differential equation with these boundary conditions is the Airy function

\[w=\textrm{Ai}\left(z\right),\]

which gives a means to check the integrator using `special.airy`.

First, convert this ODE into standard form by setting \(\mathbf{y}=\left[\frac{dw}{dz},w\right]\) and \(t=z\). Thus, the differential equation becomes

\[\begin{split}\frac{d\mathbf{y}}{dt}=\left[\begin{array}{c} ty_{1}\\ y_{0}\end{array}\right]=\left[\begin{array}{cc} 0 & t\\ 1 & 0\end{array}\right]\left[\begin{array}{c} y_{0}\\ y_{1}\end{array}\right]=\left[\begin{array}{cc} 0 & t\\ 1 & 0\end{array}\right]\mathbf{y}.\end{split}\]

In other words,

\[\mathbf{f}\left(\mathbf{y},t\right)=\mathbf{A}\left(t\right)\mathbf{y}.\]

As an interesting reminder, if \(\mathbf{A}\left(t\right)\) commutes with \(\int_{0}^{t}\mathbf{A}\left(\tau\right)\, d\tau\) under matrix multiplication, then this linear differential equation has an exact solution using the matrix exponential:

\[\mathbf{y}\left(t\right)=\exp\left(\int_{0}^{t}\mathbf{A}\left(\tau\right)d\tau\right)\mathbf{y}\left(0\right),\]

However, in this case, \(\mathbf{A}\left(t\right)\) and its integral do not commute.

This differential equation can be solved using the function `solve_ivp`. It requires the derivative, fprime, the time span [t_start, t_end] and the initial conditions vector, y0, as input arguments and returns an object whose y field is an array with consecutive solution values as columns. The initial conditions are therefore given in the first output column.

    >>> from scipy.integrate import solve_ivp
    >>> from scipy.special import gamma, airy
    >>> y1_0 = +1 / 3**(2/3) / gamma(2/3)
    >>> y0_0 = -1 / 3**(1/3) / gamma(1/3)
    >>> y0 = [y0_0, y1_0]
    >>> def func(t, y):
    ...     return [t*y[1],y[0]]
    ...
    >>> t_span = [0, 4]
    >>> sol1 = solve_ivp(func, t_span, y0)
    >>> print("sol1.t: {}".format(sol1.t))
    sol1.t:    [0.         0.10097672 1.04643602 1.91060117 2.49872472 3.08684827
     3.62692846 4.        ]

As it can be seen `solve_ivp` determines its time steps automatically if not specified otherwise. To compare the solution of `solve_ivp` with the airy function the time vector created by `solve_ivp` is passed to the airy function.

    >>> print("sol1.y[1]: {}".format(sol1.y[1]))
    sol1.y[1]: [0.35502805 0.328952   0.12801343 0.04008508 0.01601291 0.00623879
     0.00356316 0.00405982]
    >>> print("airy(sol.t)[0]:  {}".format(airy(sol1.t)[0]))
    airy(sol.t)[0]: [0.35502805 0.328952   0.12804768 0.03995804 0.01575943 0.00562799
     0.00201689 0.00095156]

The solution of `solve_ivp` with its standard parameters shows a big deviation to the airy function. To minimize this deviation, relative and absolute tolerances can be used.

    >>> rtol, atol = (1e-8, 1e-8)
    >>> sol2 = solve_ivp(func, t_span, y0, rtol=rtol, atol=atol)
    >>> print("sol2.y[1][::6]: {}".format(sol2.y[1][0::6]))
    sol2.y[1][::6]: [0.35502805 0.19145234 0.06368989 0.0205917  0.00554734 0.00106409]
    >>> print("airy(sol2.t)[0][::6]: {}".format(airy(sol2.t)[0][::6]))
    airy(sol2.t)[0][::6]: [0.35502805 0.19145234 0.06368989 0.0205917  0.00554733 0.00106406]

To specify user defined time points for the solution of `solve_ivp`, `solve_ivp` offers two possibilities that can also be used complementarily. By passing the t_eval option to the function call `solve_ivp` returns the solutions of these time points of t_eval in its output.

    >>> import numpy as np
    >>> t = np.linspace(0, 4, 100)
    >>> sol3 = solve_ivp(func, t_span, y0, t_eval=t)

If the jacobian matrix of function is known, it can be passed to the `solve_ivp` to achieve better results. Please be aware however that the default integration method `RK45` does not support jacobian matrices and thereby another integration method has to be chosen. One of the integration methods that support a jacobian matrix is the for example the `Radau` method of following example.

    >>> def gradient(t, y):
    ...     return [[0,t], [1,0]]
    >>> sol4 = solve_ivp(func, t_span, y0, method='Radau', jac=gradient)

### Solving a system with a banded Jacobian matrix

`odeint` can be told that the Jacobian is banded. For a large system of differential equations that are known to be stiff, this can improve performance significantly.

As an example, we’ll solve the 1-D Gray-Scott partial differential equations using the method of lines [MOL]. The Gray-Scott equations for the functions \(u(x, t)\) and \(v(x, t)\) on the interval \(x \in [0, L]\) are

\[\begin{split}\begin{split} \frac{\partial u}{\partial t} = D_u \frac{\partial^2 u}{\partial x^2} - uv^2 + f(1-u) \\ \frac{\partial v}{\partial t} = D_v \frac{\partial^2 v}{\partial x^2} + uv^2 - (f + k)v \\ \end{split}\end{split}\]

where \(D_u\) and \(D_v\) are the diffusion coefficients of the components \(u\) and \(v\), respectively, and \(f\) and \(k\) are constants. (For more information about the system, see http://groups.csail.mit.edu/mac/projects/amorphous/GrayScott/)

We’ll assume Neumann (i.e., “no flux”) boundary conditions:

\[\frac{\partial u}{\partial x}(0,t) = 0, \quad \frac{\partial v}{\partial x}(0,t) = 0, \quad \frac{\partial u}{\partial x}(L,t) = 0, \quad \frac{\partial v}{\partial x}(L,t) = 0\]

To apply the method of lines, we discretize the \(x\) variable by defining the uniformly spaced grid of \(N\) points \(\left\{x_0, x_1, \ldots, x_{N-1}\right\}\), with \(x_0 = 0\) and \(x_{N-1} = L\). We define \(u_j(t) \equiv u(x_k, t)\) and \(v_j(t) \equiv v(x_k, t)\), and replace the \(x\) derivatives with finite differences. That is,

\[\frac{\partial^2 u}{\partial x^2}(x_j, t) \rightarrow \frac{u_{j-1}(t) - 2 u_{j}(t) + u_{j+1}(t)}{(\Delta x)^2}\]

We then have a system of \(2N\) ordinary differential equations:

(1)\[\begin{split} \begin{split} \frac{du_j}{dt} = \frac{D_u}{(\Delta x)^2} \left(u_{j-1} - 2 u_{j} + u_{j+1}\right) -u_jv_j^2 + f(1 - u_j) \\ \frac{dv_j}{dt} = \frac{D_v}{(\Delta x)^2} \left(v_{j-1} - 2 v_{j} + v_{j+1}\right) + u_jv_j^2 - (f + k)v_j \end{split}\end{split}\]

For convenience, the \((t)\) arguments have been dropped.

To enforce the boundary conditions, we introduce “ghost” points \(x_{-1}\) and \(x_N\), and define \(u_{-1}(t) \equiv u_1(t)\), \(u_N(t) \equiv u_{N-2}(t)\); \(v_{-1}(t)\) and \(v_N(t)\) are defined analogously.

Then

(2)\[\begin{split} \begin{split} \frac{du_0}{dt} = \frac{D_u}{(\Delta x)^2} \left(2u_{1} - 2 u_{0}\right) -u_0v_0^2 + f(1 - u_0) \\ \frac{dv_0}{dt} = \frac{D_v}{(\Delta x)^2} \left(2v_{1} - 2 v_{0}\right) + u_0v_0^2 - (f + k)v_0 \end{split}\end{split}\]

and

(3)\[\begin{split} \begin{split} \frac{du_{N-1}}{dt} = \frac{D_u}{(\Delta x)^2} \left(2u_{N-2} - 2 u_{N-1}\right) -u_{N-1}v_{N-1}^2 + f(1 - u_{N-1}) \\ \frac{dv_{N-1}}{dt} = \frac{D_v}{(\Delta x)^2} \left(2v_{N-2} - 2 v_{N-1}\right) + u_{N-1}v_{N-1}^2 - (f + k)v_{N-1} \end{split}\end{split}\]

Our complete system of \(2N\) ordinary differential equations is (1) for \(k = 1, 2, \ldots, N-2\), along with (2) and (3).

We can now starting implementing this system in code. We must combine \(\{u_k\}\) and \(\{v_k\}\) into a single vector of length \(2N\). The two obvious choices are \(\{u_0, u_1, \ldots, u_{N-1}, v_0, v_1, \ldots, v_{N-1}\}\) and \(\{u_0, v_0, u_1, v_1, \ldots, u_{N-1}, v_{N-1}\}\). Mathematically, it does not matter, but the choice affects how efficiently `odeint` can solve the system. The reason is in how the order affects the pattern of the nonzero elements of the Jacobian matrix.

When the variables are ordered as \(\{u_0, u_1, \ldots, u_{N-1}, v_0, v_1, \ldots, v_{N-1}\}\), the pattern of nonzero elements of the Jacobian matrix is

\[\begin{split}\begin{smallmatrix} * & * & 0 & 0 & 0 & 0 & 0 & * & 0 & 0 & 0 & 0 & 0 & 0 \\ * & * & * & 0 & 0 & 0 & 0 & 0 & * & 0 & 0 & 0 & 0 & 0 \\ 0 & * & * & * & 0 & 0 & 0 & 0 & 0 & * & 0 & 0 & 0 & 0 \\ 0 & 0 & * & * & * & 0 & 0 & 0 & 0 & 0 & * & 0 & 0 & 0 \\ 0 & 0 & 0 & * & * & * & 0 & 0 & 0 & 0 & 0 & * & 0 & 0 \\ 0 & 0 & 0 & 0 & * & * & * & 0 & 0 & 0 & 0 & 0 & * & 0 \\ 0 & 0 & 0 & 0 & 0 & * & * & 0 & 0 & 0 & 0 & 0 & 0 & * \\ * & 0 & 0 & 0 & 0 & 0 & 0 & * & * & 0 & 0 & 0 & 0 & 0 \\ 0 & * & 0 & 0 & 0 & 0 & 0 & * & * & * & 0 & 0 & 0 & 0 \\ 0 & 0 & * & 0 & 0 & 0 & 0 & 0 & * & * & * & 0 & 0 & 0 \\ 0 & 0 & 0 & * & 0 & 0 & 0 & 0 & 0 & * & * & * & 0 & 0 \\ 0 & 0 & 0 & 0 & * & 0 & 0 & 0 & 0 & 0 & * & * & * & 0 \\ 0 & 0 & 0 & 0 & 0 & * & 0 & 0 & 0 & 0 & 0 & * & * & * \\ 0 & 0 & 0 & 0 & 0 & 0 & * & 0 & 0 & 0 & 0 & ) & * & * \\ \end{smallmatrix}\end{split}\]

The Jacobian pattern with variables interleaved as \(\{u_0, v_0, u_1, v_1, \ldots, u_{N-1}, v_{N-1}\}\) is

\[\begin{split}\begin{smallmatrix} * & * & * & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ * & * & 0 & * & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ * & 0 & * & * & * & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & * & * & * & 0 & * & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & * & 0 & * & * & * & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & * & * & * & 0 & * & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & * & 0 & * & * & * & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & * & * & * & 0 & * & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & * & 0 & * & * & * & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & * & * & * & 0 & * & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & * & 0 & * & * & * & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & * & * & * & 0 & * \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & * & 0 & * & * \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & * & * & * \\ \end{smallmatrix}\end{split}\]

In both cases, there are just five nontrivial diagonals, but when the variables are interleaved, the bandwidth is much smaller. That is, the main diagonal and the two diagonals immediately above and the two immediately below the main diagonal are the nonzero diagonals. This is important, because the inputs `mu` and `ml` of `odeint` are the upper and lower bandwidths of the Jacobian matrix. When the variables are interleaved, `mu` and `ml` are 2. When the variables are stacked with \(\{v_k\}\) following \(\{u_k\}\), the upper and lower bandwidths are \(N\).

With that decision made, we can write the function that implements the system of differential equations.

First, we define the functions for the source and reaction terms of the system:

    def G(u, v, f, k):
        return f * (1 - u) - u*v**2

    def H(u, v, f, k):
        return -(f + k) * v + u*v**2

Next, we define the function that computes the right-hand side of the system of differential equations:

    def grayscott1d(y, t, f, k, Du, Dv, dx):
        """
        Differential equations for the 1-D Gray-Scott equations.

        The ODEs are derived using the method of lines.
        """
        # The vectors u and v are interleaved in y.  We define
        # views of u and v by slicing y.
        u = y[::2]
        v = y[1::2]

        # dydt is the return value of this function.
        dydt = np.empty_like(y)

        # Just like u and v are views of the interleaved vectors
        # in y, dudt and dvdt are views of the interleaved output
        # vectors in dydt.
        dudt = dydt[::2]
        dvdt = dydt[1::2]

        # Compute du/dt and dv/dt.  The end points and the interior points
        # are handled separately.
        dudt[0]    = G(u[0],    v[0],    f, k) + Du * (-2.0*u[0] + 2.0*u[1]) / dx**2
        dudt[1:-1] = G(u[1:-1], v[1:-1], f, k) + Du * np.diff(u,2) / dx**2
        dudt[-1]   = G(u[-1],   v[-1],   f, k) + Du * (- 2.0*u[-1] + 2.0*u[-2]) / dx**2
        dvdt[0]    = H(u[0],    v[0],    f, k) + Dv * (-2.0*v[0] + 2.0*v[1]) / dx**2
        dvdt[1:-1] = H(u[1:-1], v[1:-1], f, k) + Dv * np.diff(v,2) / dx**2
        dvdt[-1]   = H(u[-1],   v[-1],   f, k) + Dv * (-2.0*v[-1] + 2.0*v[-2]) / dx**2

        return dydt

We won’t implement a function to compute the Jacobian, but we will tell `odeint` that the Jacobian matrix is banded. This allows the underlying solver (LSODA) to avoid computing values that it knows are zero. For a large system, this improves the performance significantly, as demonstrated in the following ipython session.

First, we define the required inputs:

    In [30]: rng = np.random.default_rng()

    In [31]: y0 = rng.standard_normal(5000)

    In [32]: t = np.linspace(0, 50, 11)

    In [33]: f = 0.024

    In [34]: k = 0.055

    In [35]: Du = 0.01

    In [36]: Dv = 0.005

    In [37]: dx = 0.025

Time the computation without taking advantage of the banded structure of the Jacobian matrix:

    In [38]: %timeit sola = odeint(grayscott1d, y0, t, args=(f, k, Du, Dv, dx))
    1 loop, best of 3: 25.2 s per loop

Now set `ml=2` and `mu=2`, so `odeint` knows that the Jacobian matrix is banded:

    In [39]: %timeit solb = odeint(grayscott1d, y0, t, args=(f, k, Du, Dv, dx), ml=2, mu=2)
    10 loops, best of 3: 191 ms per loop

That is quite a bit faster!

Let’s ensure that they have computed the same result:

    In [41]: np.allclose(sola, solb)
    Out[41]: True

### References

[WPR]

https://en.wikipedia.org/wiki/Romberg’s_method

[MOL]

https://en.wikipedia.org/wiki/Method_of_lines