1Integration (:mod:`scipy.integrate`) 2==================================== 3 4.. sectionauthor:: Travis E. Oliphant 5 6.. currentmodule:: scipy.integrate 7 8The :mod:`scipy.integrate` sub-package provides several integration 9techniques including an ordinary differential equation integrator. An 10overview of the module is provided by the help command: 11 12.. literalinclude:: examples/4-1 13 14 15General integration (:func:`quad`) 16---------------------------------- 17 18The function :obj:`quad` is provided to integrate a function of one 19variable between two points. The points can be :math:`\pm\infty` 20(:math:`\pm` ``inf``) to indicate infinite limits. For example, 21suppose you wish to integrate a bessel function ``jv(2.5, x)`` along 22the interval :math:`[0, 4.5].` 23 24.. math:: 25 26 I=\int_{0}^{4.5}J_{2.5}\left(x\right)\, dx. 27 28 29This could be computed using :obj:`quad`: 30 31 >>> import scipy.integrate as integrate 32 >>> import scipy.special as special 33 >>> result = integrate.quad(lambda x: special.jv(2.5,x), 0, 4.5) 34 >>> result 35 (1.1178179380783249, 7.8663172481899801e-09) 36 37 >>> from numpy import sqrt, sin, cos, pi 38 >>> I = sqrt(2/pi)*(18.0/27*sqrt(2)*cos(4.5) - 4.0/27*sqrt(2)*sin(4.5) + 39 ... sqrt(2*pi) * special.fresnel(3/sqrt(pi))[0]) 40 >>> I 41 1.117817938088701 42 43 >>> print(abs(result[0]-I)) 44 1.03761443881e-11 45 46The first argument to quad is a "callable" Python object (i.e., a 47function, method, or class instance). Notice the use of a lambda- 48function in this case as the argument. The next two arguments are the 49limits of integration. The return value is a tuple, with the first 50element holding the estimated value of the integral and the second 51element holding an upper bound on the error. Notice, that in this 52case, the true value of this integral is 53 54.. math:: 55 56 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), 57 58where 59 60.. math:: 61 62 \textrm{Si}\left(x\right)=\int_{0}^{x}\sin\left(\frac{\pi}{2}t^{2}\right)\, dt. 63 64is the Fresnel sine integral. Note that the numerically-computed integral is 65within :math:`1.04\times10^{-11}` of the exact result --- well below the 66reported error bound. 67 68 69If the function to integrate takes additional parameters, they can be provided 70in the `args` argument. Suppose that the following integral shall be calculated: 71 72.. math:: 73 74 I(a,b)=\int_{0}^{1} ax^2+b \, dx. 75 76 77This integral can be evaluated by using the following code: 78 79>>> from scipy.integrate import quad 80>>> def integrand(x, a, b): 81... return a*x**2 + b 82... 83>>> a = 2 84>>> b = 1 85>>> I = quad(integrand, 0, 1, args=(a,b)) 86>>> I 87(1.6666666666666667, 1.8503717077085944e-14) 88 89 90Infinite inputs are also allowed in :obj:`quad` by using :math:`\pm` 91``inf`` as one of the arguments. For example, suppose that a numerical 92value for the exponential integral: 93 94.. math:: 95 96 E_{n}\left(x\right)=\int_{1}^{\infty}\frac{e^{-xt}}{t^{n}}\, dt. 97 98is desired (and the fact that this integral can be computed as 99``special.expn(n,x)`` is forgotten). The functionality of the function 100:obj:`special.expn <scipy.special.expn>` can be replicated by defining a new function 101``vec_expint`` based on the routine :obj:`quad`: 102 103 >>> from scipy.integrate import quad 104 >>> def integrand(t, n, x): 105 ... return np.exp(-x*t) / t**n 106 ... 107 108 >>> def expint(n, x): 109 ... return quad(integrand, 1, np.inf, args=(n, x))[0] 110 ... 111 112 >>> vec_expint = np.vectorize(expint) 113 114 >>> vec_expint(3, np.arange(1.0, 4.0, 0.5)) 115 array([ 0.1097, 0.0567, 0.0301, 0.0163, 0.0089, 0.0049]) 116 >>> import scipy.special as special 117 >>> special.expn(3, np.arange(1.0,4.0,0.5)) 118 array([ 0.1097, 0.0567, 0.0301, 0.0163, 0.0089, 0.0049]) 119 120The function which is integrated can even use the quad argument (though the 121error bound may underestimate the error due to possible numerical error in the 122integrand from the use of :obj:`quad` ). The integral in this case is 123 124.. math:: 125 126 I_{n}=\int_{0}^{\infty}\int_{1}^{\infty}\frac{e^{-xt}}{t^{n}}\, dt\, dx=\frac{1}{n}. 127 128>>> result = quad(lambda x: expint(3, x), 0, np.inf) 129>>> print(result) 130(0.33333333324560266, 2.8548934485373678e-09) 131 132>>> I3 = 1.0/3.0 133>>> print(I3) 1340.333333333333 135 136>>> print(I3 - result[0]) 1378.77306560731e-11 138 139This last example shows that multiple integration can be handled using 140repeated calls to :func:`quad`. 141 142 143General multiple integration (:func:`dblquad`, :func:`tplquad`, :func:`nquad`) 144------------------------------------------------------------------------------ 145 146The mechanics for double and triple integration have been wrapped up into the 147functions :obj:`dblquad` and :obj:`tplquad`. These functions take the function 148to integrate and four, or six arguments, respectively. The limits of all 149inner integrals need to be defined as functions. 150 151An example of using double integration to compute several values of 152:math:`I_{n}` is shown below: 153 154 >>> from scipy.integrate import quad, dblquad 155 >>> def I(n): 156 ... return dblquad(lambda t, x: np.exp(-x*t)/t**n, 0, np.inf, lambda x: 1, lambda x: np.inf) 157 ... 158 159 >>> print(I(4)) 160 (0.2500000000043577, 1.29830334693681e-08) 161 >>> print(I(3)) 162 (0.33333333325010883, 1.3888461883425516e-08) 163 >>> print(I(2)) 164 (0.4999999999985751, 1.3894083651858995e-08) 165 166 167As example for non-constant limits consider the integral 168 169.. math:: 170 171 I=\int_{y=0}^{1/2}\int_{x=0}^{1-2y} x y \, dx\, dy=\frac{1}{96}. 172 173 174This integral can be evaluated using the expression below (Note the use of the 175non-constant lambda functions for the upper limit of the inner integral): 176 177>>> from scipy.integrate import dblquad 178>>> area = dblquad(lambda x, y: x*y, 0, 0.5, lambda x: 0, lambda x: 1-2*x) 179>>> area 180(0.010416666666666668, 1.1564823173178715e-16) 181 182 183For n-fold integration, scipy provides the function :obj:`nquad`. The 184integration bounds are an iterable object: either a list of constant bounds, 185or a list of functions for the non-constant integration bounds. The order of 186integration (and therefore the bounds) is from the innermost integral to the 187outermost one. 188 189The integral from above 190 191.. math:: 192 193 I_{n}=\int_{0}^{\infty}\int_{1}^{\infty}\frac{e^{-xt}}{t^{n}}\, dt\, dx=\frac{1}{n} 194 195can be calculated as 196 197>>> from scipy import integrate 198>>> N = 5 199>>> def f(t, x): 200... return np.exp(-x*t) / t**N 201... 202>>> integrate.nquad(f, [[1, np.inf],[0, np.inf]]) 203(0.20000000000002294, 1.2239614263187945e-08) 204 205Note that the order of arguments for `f` must match the order of the 206integration bounds; i.e., the inner integral with respect to :math:`t` is on 207the interval :math:`[1, \infty]` and the outer integral with respect to 208:math:`x` is on the interval :math:`[0, \infty]`. 209 210Non-constant integration bounds can be treated in a similar manner; the 211example from above 212 213.. math:: 214 215 I=\int_{y=0}^{1/2}\int_{x=0}^{1-2y} x y \, dx\, dy=\frac{1}{96}. 216 217can be evaluated by means of 218 219>>> from scipy import integrate 220>>> def f(x, y): 221... return x*y 222... 223>>> def bounds_y(): 224... return [0, 0.5] 225... 226>>> def bounds_x(y): 227... return [0, 1-2*y] 228... 229>>> integrate.nquad(f, [bounds_x, bounds_y]) 230(0.010416666666666668, 4.101620128472366e-16) 231 232which is the same result as before. 233 234Gaussian quadrature 235------------------- 236 237A few functions are also provided in order to perform simple Gaussian 238quadrature over a fixed interval. The first is :obj:`fixed_quad`, which 239performs fixed-order Gaussian quadrature. The second function is 240:obj:`quadrature`, which performs Gaussian quadrature of multiple 241orders until the difference in the integral estimate is beneath some 242tolerance supplied by the user. These functions both use the module 243``scipy.special.orthogonal``, which can calculate the roots and quadrature 244weights of a large variety of orthogonal polynomials (the polynomials 245themselves are available as special functions returning instances of 246the polynomial class --- e.g., :obj:`special.legendre <scipy.special.legendre>`). 247 248 249Romberg Integration 250------------------- 251 252Romberg's method [WPR]_ is another method for numerically evaluating an 253integral. See the help function for :func:`romberg` for further details. 254 255 256Integrating using Samples 257------------------------- 258 259If the samples are equally-spaced and the number of samples available 260is :math:`2^{k}+1` for some integer :math:`k`, then Romberg :obj:`romb` 261integration can be used to obtain high-precision estimates of the 262integral using the available samples. Romberg integration uses the 263trapezoid rule at step-sizes related by a power of two and then 264performs Richardson extrapolation on these estimates to approximate 265the integral with a higher degree of accuracy. 266 267In case of arbitrary spaced samples, the two functions :obj:`trapezoid` 268and :obj:`simpson` are available. They are using Newton-Coates formulas 269of order 1 and 2 respectively to perform integration. The trapezoidal rule 270approximates the function as a straight line between adjacent points, while 271Simpson's rule approximates the function between three adjacent points as a 272parabola. 273 274For an odd number of samples that are equally spaced Simpson's rule is exact 275if the function is a polynomial of order 3 or less. If the samples are not 276equally spaced, then the result is exact only if the function is a polynomial 277of order 2 or less. 278 279>>> import numpy as np 280>>> def f1(x): 281... return x**2 282... 283>>> def f2(x): 284... return x**3 285... 286>>> x = np.array([1,3,4]) 287>>> y1 = f1(x) 288>>> from scipy import integrate 289>>> I1 = integrate.simpson(y1, x) 290>>> print(I1) 29121.0 292 293 294This corresponds exactly to 295 296.. math:: 297 298 \int_{1}^{4} x^2 \, dx = 21, 299 300whereas integrating the second function 301 302>>> y2 = f2(x) 303>>> I2 = integrate.simpson(y2, x) 304>>> print(I2) 30561.5 306 307does not correspond to 308 309.. math:: 310 311 \int_{1}^{4} x^3 \, dx = 63.75 312 313because the order of the polynomial in f2 is larger than two. 314 315.. _quad-callbacks: 316 317Faster integration using low-level callback functions 318----------------------------------------------------- 319 320A user desiring reduced integration times may pass a C function 321pointer through `scipy.LowLevelCallable` to `quad`, `dblquad`, 322`tplquad` or `nquad` and it will be integrated and return a result in 323Python. The performance increase here arises from two factors. The 324primary improvement is faster function evaluation, which is provided 325by compilation of the function itself. Additionally we have a speedup 326provided by the removal of function calls between C and Python in 327:obj:`quad`. This method may provide a speed improvements of ~2x for 328trivial functions such as sine but can produce a much more noticeable 329improvements (10x+) for more complex functions. This feature then, is 330geared towards a user with numerically intensive integrations willing 331to write a little C to reduce computation time significantly. 332 333The approach can be used, for example, via `ctypes` in a few simple steps: 334 3351.) Write an integrand function in C with the function signature 336``double f(int n, double *x, void *user_data)``, where ``x`` is an 337array containing the point the function f is evaluated at, and ``user_data`` 338to arbitrary additional data you want to provide. 339 340.. code-block:: c 341 342 /* testlib.c */ 343 double f(int n, double *x, void *user_data) { 344 double c = *(double *)user_data; 345 return c + x[0] - x[1] * x[2]; /* corresponds to c + x - y * z */ 346 } 347 3482.) Now compile this file to a shared/dynamic library (a quick search will help 349with this as it is OS-dependent). The user must link any math libraries, 350etc., used. On linux this looks like:: 351 352 $ gcc -shared -fPIC -o testlib.so testlib.c 353 354The output library will be referred to as ``testlib.so``, but it may have a 355different file extension. A library has now been created that can be loaded 356into Python with `ctypes`. 357 3583.) Load shared library into Python using `ctypes` and set ``restypes`` and 359``argtypes`` - this allows SciPy to interpret the function correctly: 360 361.. code:: python 362 363 import os, ctypes 364 from scipy import integrate, LowLevelCallable 365 366 lib = ctypes.CDLL(os.path.abspath('testlib.so')) 367 lib.f.restype = ctypes.c_double 368 lib.f.argtypes = (ctypes.c_int, ctypes.POINTER(ctypes.c_double), ctypes.c_void_p) 369 370 c = ctypes.c_double(1.0) 371 user_data = ctypes.cast(ctypes.pointer(c), ctypes.c_void_p) 372 373 func = LowLevelCallable(lib.f, user_data) 374 375The last ``void *user_data`` in the function is optional and can be omitted 376(both in the C function and ctypes argtypes) if not needed. Note that the 377coordinates are passed in as an array of doubles rather than a separate argument. 378 3794.) Now integrate the library function as normally, here using `nquad`: 380 381>>> integrate.nquad(func, [[0, 10], [-10, 0], [-1, 1]]) 382(1200.0, 1.1102230246251565e-11) 383 384The Python tuple is returned as expected in a reduced amount of time. All 385optional parameters can be used with this method including specifying 386singularities, infinite bounds, etc. 387 388Ordinary differential equations (:func:`solve_ivp`) 389--------------------------------------------------- 390 391Integrating a set of ordinary differential equations (ODEs) given 392initial conditions is another useful example. The function 393:obj:`solve_ivp` is available in SciPy for integrating a first-order 394vector differential equation: 395 396.. math:: 397 398 \frac{d\mathbf{y}}{dt}=\mathbf{f}\left(\mathbf{y},t\right), 399 400given initial conditions :math:`\mathbf{y}\left(0\right)=y_{0}`, where 401:math:`\mathbf{y}` is a length :math:`N` vector and :math:`\mathbf{f}` 402is a mapping from :math:`\mathcal{R}^{N}` to :math:`\mathcal{R}^{N}.` 403A higher-order ordinary differential equation can always be reduced to 404a differential equation of this type by introducing intermediate 405derivatives into the :math:`\mathbf{y}` vector. 406 407For example, suppose it is desired to find the solution to the 408following second-order differential equation: 409 410.. math:: 411 412 \frac{d^{2}w}{dz^{2}}-zw(z)=0 413 414with initial conditions :math:`w\left(0\right)=\frac{1}{\sqrt[3]{3^{2}}\Gamma\left(\frac{2}{3}\right)}` and :math:`\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 415boundary conditions is the Airy function 416 417.. math:: 418 419 w=\textrm{Ai}\left(z\right), 420 421which gives a means to check the integrator using :func:`special.airy <scipy.special.airy>`. 422 423First, convert this ODE into standard form by setting 424:math:`\mathbf{y}=\left[\frac{dw}{dz},w\right]` and :math:`t=z`. Thus, 425the differential equation becomes 426 427.. math:: 428 429 \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}. 430 431In other words, 432 433.. math:: 434 435 \mathbf{f}\left(\mathbf{y},t\right)=\mathbf{A}\left(t\right)\mathbf{y}. 436 437As an interesting reminder, if :math:`\mathbf{A}\left(t\right)` 438commutes with :math:`\int_{0}^{t}\mathbf{A}\left(\tau\right)\, d\tau` 439under matrix multiplication, then this linear differential equation 440has an exact solution using the matrix exponential: 441 442.. math:: 443 444 \mathbf{y}\left(t\right)=\exp\left(\int_{0}^{t}\mathbf{A}\left(\tau\right)d\tau\right)\mathbf{y}\left(0\right), 445 446However, in this case, :math:`\mathbf{A}\left(t\right)` and its integral do not commute. 447 448This differential equation can be solved using the function :obj:`solve_ivp`. 449It requires the derivative, *fprime*, the time span `[t_start, t_end]` 450and the initial conditions vector, *y0*, as input arguments and returns 451an object whose *y* field is an array with consecutive solution values as 452columns. The initial conditions are therefore given in the first output column. 453 454>>> from scipy.integrate import solve_ivp 455>>> from scipy.special import gamma, airy 456>>> y1_0 = +1 / 3**(2/3) / gamma(2/3) 457>>> y0_0 = -1 / 3**(1/3) / gamma(1/3) 458>>> y0 = [y0_0, y1_0] 459>>> def func(t, y): 460... return [t*y[1],y[0]] 461... 462>>> t_span = [0, 4] 463>>> sol1 = solve_ivp(func, t_span, y0) 464>>> print("sol1.t: {}".format(sol1.t)) 465sol1.t: [0. 0.10097672 1.04643602 1.91060117 2.49872472 3.08684827 466 3.62692846 4. ] 467 468As it can be seen `solve_ivp` determines its time steps automatically if not 469specified otherwise. To compare the solution of `solve_ivp` with the `airy` 470function the time vector created by `solve_ivp` is passed to the `airy` function. 471 472>>> print("sol1.y[1]: {}".format(sol1.y[1])) 473sol1.y[1]: [0.35502805 0.328952 0.12801343 0.04008508 0.01601291 0.00623879 474 0.00356316 0.00405982] 475>>> print("airy(sol.t)[0]: {}".format(airy(sol1.t)[0])) 476airy(sol.t)[0]: [0.35502805 0.328952 0.12804768 0.03995804 0.01575943 0.00562799 477 0.00201689 0.00095156] 478 479The solution of `solve_ivp` with its standard parameters shows a big deviation 480to the airy function. To minimize this deviation, relative and absolute 481tolerances can be used. 482 483>>> rtol, atol = (1e-8, 1e-8) 484>>> sol2 = solve_ivp(func, t_span, y0, rtol=rtol, atol=atol) 485>>> print("sol2.y[1][::6]: {}".format(sol2.y[1][0::6])) 486sol2.y[1][::6]: [0.35502805 0.19145234 0.06368989 0.0205917 0.00554734 0.00106409] 487>>> print("airy(sol2.t)[0][::6]: {}".format(airy(sol2.t)[0][::6])) 488airy(sol2.t)[0][::6]: [0.35502805 0.19145234 0.06368989 0.0205917 0.00554733 0.00106406] 489 490To specify user defined time points for the solution of `solve_ivp`, `solve_ivp` 491offers two possibilities that can also be used complementarily. By passing the `t_eval` 492option to the function call `solve_ivp` returns the solutions of these time points 493of `t_eval` in its output. 494 495>>> import numpy as np 496>>> t = np.linspace(0, 4, 100) 497>>> sol3 = solve_ivp(func, t_span, y0, t_eval=t) 498 499If the jacobian matrix of function is known, it can be passed to the `solve_ivp` 500to achieve better results. Please be aware however that the default integration method 501`RK45` does not support jacobian matrices and thereby another integration method has 502to be chosen. One of the integration methods that support a jacobian matrix is the for 503example the `Radau` method of following example. 504 505>>> def gradient(t, y): 506... return [[0,t], [1,0]] 507>>> sol4 = solve_ivp(func, t_span, y0, method='Radau', jac=gradient) 508 509Solving a system with a banded Jacobian matrix 510~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 511 512`odeint` can be told that the Jacobian is *banded*. For a large 513system of differential equations that are known to be stiff, this 514can improve performance significantly. 515 516As an example, we'll solve the 1-D Gray-Scott partial 517differential equations using the method of lines [MOL]_. The Gray-Scott equations 518for the functions :math:`u(x, t)` and :math:`v(x, t)` on the interval 519:math:`x \in [0, L]` are 520 521.. math:: 522 523 \begin{split} 524 \frac{\partial u}{\partial t} = D_u \frac{\partial^2 u}{\partial x^2} - uv^2 + f(1-u) \\ 525 \frac{\partial v}{\partial t} = D_v \frac{\partial^2 v}{\partial x^2} + uv^2 - (f + k)v \\ 526 \end{split} 527 528where :math:`D_u` and :math:`D_v` are the diffusion coefficients of the 529components :math:`u` and :math:`v`, respectively, and :math:`f` and :math:`k` 530are constants. (For more information about the system, see 531http://groups.csail.mit.edu/mac/projects/amorphous/GrayScott/) 532 533We'll assume Neumann (i.e., "no flux") boundary conditions: 534 535.. math:: 536 537 \frac{\partial u}{\partial x}(0,t) = 0, \quad 538 \frac{\partial v}{\partial x}(0,t) = 0, \quad 539 \frac{\partial u}{\partial x}(L,t) = 0, \quad 540 \frac{\partial v}{\partial x}(L,t) = 0 541 542To apply the method of lines, we discretize the :math:`x` variable by defining 543the uniformly spaced grid of :math:`N` points :math:`\left\{x_0, x_1, \ldots, x_{N-1}\right\}`, with 544:math:`x_0 = 0` and :math:`x_{N-1} = L`. 545We define :math:`u_j(t) \equiv u(x_k, t)` and :math:`v_j(t) \equiv v(x_k, t)`, and 546replace the :math:`x` derivatives with finite differences. That is, 547 548.. math:: 549 550 \frac{\partial^2 u}{\partial x^2}(x_j, t) \rightarrow 551 \frac{u_{j-1}(t) - 2 u_{j}(t) + u_{j+1}(t)}{(\Delta x)^2} 552 553We then have a system of :math:`2N` ordinary differential equations: 554 555.. math:: 556 :label: interior 557 558 \begin{split} 559 \frac{du_j}{dt} = \frac{D_u}{(\Delta x)^2} \left(u_{j-1} - 2 u_{j} + u_{j+1}\right) 560 -u_jv_j^2 + f(1 - u_j) \\ 561 \frac{dv_j}{dt} = \frac{D_v}{(\Delta x)^2} \left(v_{j-1} - 2 v_{j} + v_{j+1}\right) 562 + u_jv_j^2 - (f + k)v_j 563 \end{split} 564 565For convenience, the :math:`(t)` arguments have been dropped. 566 567To enforce the boundary conditions, we introduce "ghost" points 568:math:`x_{-1}` and :math:`x_N`, and define :math:`u_{-1}(t) \equiv u_1(t)`, 569:math:`u_N(t) \equiv u_{N-2}(t)`; :math:`v_{-1}(t)` and :math:`v_N(t)` 570are defined analogously. 571 572Then 573 574.. math:: 575 :label: boundary0 576 577 \begin{split} 578 \frac{du_0}{dt} = \frac{D_u}{(\Delta x)^2} \left(2u_{1} - 2 u_{0}\right) 579 -u_0v_0^2 + f(1 - u_0) \\ 580 \frac{dv_0}{dt} = \frac{D_v}{(\Delta x)^2} \left(2v_{1} - 2 v_{0}\right) 581 + u_0v_0^2 - (f + k)v_0 582 \end{split} 583 584and 585 586.. math:: 587 :label: boundaryL 588 589 \begin{split} 590 \frac{du_{N-1}}{dt} = \frac{D_u}{(\Delta x)^2} \left(2u_{N-2} - 2 u_{N-1}\right) 591 -u_{N-1}v_{N-1}^2 + f(1 - u_{N-1}) \\ 592 \frac{dv_{N-1}}{dt} = \frac{D_v}{(\Delta x)^2} \left(2v_{N-2} - 2 v_{N-1}\right) 593 + u_{N-1}v_{N-1}^2 - (f + k)v_{N-1} 594 \end{split} 595 596Our complete system of :math:`2N` ordinary differential equations is :eq:`interior` 597for :math:`k = 1, 2, \ldots, N-2`, along with :eq:`boundary0` and :eq:`boundaryL`. 598 599We can now starting implementing this system in code. We must combine 600:math:`\{u_k\}` and :math:`\{v_k\}` into a single vector of length :math:`2N`. 601The two obvious choices are 602:math:`\{u_0, u_1, \ldots, u_{N-1}, v_0, v_1, \ldots, v_{N-1}\}` 603and 604:math:`\{u_0, v_0, u_1, v_1, \ldots, u_{N-1}, v_{N-1}\}`. 605Mathematically, it does not matter, but the choice affects how 606efficiently `odeint` can solve the system. The reason is in how 607the order affects the pattern of the nonzero elements of the Jacobian matrix. 608 609 610When the variables are ordered 611as :math:`\{u_0, u_1, \ldots, u_{N-1}, v_0, v_1, \ldots, v_{N-1}\}`, 612the pattern of nonzero elements of the Jacobian matrix is 613 614.. math:: 615 616 \begin{smallmatrix} 617 * & * & 0 & 0 & 0 & 0 & 0 & * & 0 & 0 & 0 & 0 & 0 & 0 \\ 618 * & * & * & 0 & 0 & 0 & 0 & 0 & * & 0 & 0 & 0 & 0 & 0 \\ 619 0 & * & * & * & 0 & 0 & 0 & 0 & 0 & * & 0 & 0 & 0 & 0 \\ 620 0 & 0 & * & * & * & 0 & 0 & 0 & 0 & 0 & * & 0 & 0 & 0 \\ 621 0 & 0 & 0 & * & * & * & 0 & 0 & 0 & 0 & 0 & * & 0 & 0 \\ 622 0 & 0 & 0 & 0 & * & * & * & 0 & 0 & 0 & 0 & 0 & * & 0 \\ 623 0 & 0 & 0 & 0 & 0 & * & * & 0 & 0 & 0 & 0 & 0 & 0 & * \\ 624 * & 0 & 0 & 0 & 0 & 0 & 0 & * & * & 0 & 0 & 0 & 0 & 0 \\ 625 0 & * & 0 & 0 & 0 & 0 & 0 & * & * & * & 0 & 0 & 0 & 0 \\ 626 0 & 0 & * & 0 & 0 & 0 & 0 & 0 & * & * & * & 0 & 0 & 0 \\ 627 0 & 0 & 0 & * & 0 & 0 & 0 & 0 & 0 & * & * & * & 0 & 0 \\ 628 0 & 0 & 0 & 0 & * & 0 & 0 & 0 & 0 & 0 & * & * & * & 0 \\ 629 0 & 0 & 0 & 0 & 0 & * & 0 & 0 & 0 & 0 & 0 & * & * & * \\ 630 0 & 0 & 0 & 0 & 0 & 0 & * & 0 & 0 & 0 & 0 & ) & * & * \\ 631 \end{smallmatrix} 632 633The Jacobian pattern with variables interleaved 634as :math:`\{u_0, v_0, u_1, v_1, \ldots, u_{N-1}, v_{N-1}\}` is 635 636.. math:: 637 \begin{smallmatrix} 638 * & * & * & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 639 * & * & 0 & * & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 640 * & 0 & * & * & * & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 641 0 & * & * & * & 0 & * & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 642 0 & 0 & * & 0 & * & * & * & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 643 0 & 0 & 0 & * & * & * & 0 & * & 0 & 0 & 0 & 0 & 0 & 0 \\ 644 0 & 0 & 0 & 0 & * & 0 & * & * & * & 0 & 0 & 0 & 0 & 0 \\ 645 0 & 0 & 0 & 0 & 0 & * & * & * & 0 & * & 0 & 0 & 0 & 0 \\ 646 0 & 0 & 0 & 0 & 0 & 0 & * & 0 & * & * & * & 0 & 0 & 0 \\ 647 0 & 0 & 0 & 0 & 0 & 0 & 0 & * & * & * & 0 & * & 0 & 0 \\ 648 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & * & 0 & * & * & * & 0 \\ 649 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & * & * & * & 0 & * \\ 650 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & * & 0 & * & * \\ 651 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & * & * & * \\ 652 \end{smallmatrix} 653 654In both cases, there are just five nontrivial diagonals, but 655when the variables are interleaved, the bandwidth is much 656smaller. 657That is, the main diagonal and the two diagonals immediately 658above and the two immediately below the main diagonal 659are the nonzero diagonals. 660This is important, because the inputs ``mu`` and ``ml`` 661of `odeint` are the upper and lower bandwidths of the 662Jacobian matrix. When the variables are interleaved, 663``mu`` and ``ml`` are 2. When the variables are stacked 664with :math:`\{v_k\}` following :math:`\{u_k\}`, the upper 665and lower bandwidths are :math:`N`. 666 667With that decision made, we can write the function that 668implements the system of differential equations. 669 670First, we define the functions for the source and reaction 671terms of the system:: 672 673 def G(u, v, f, k): 674 return f * (1 - u) - u*v**2 675 676 def H(u, v, f, k): 677 return -(f + k) * v + u*v**2 678 679Next, we define the function that computes the right-hand side 680of the system of differential equations:: 681 682 def grayscott1d(y, t, f, k, Du, Dv, dx): 683 """ 684 Differential equations for the 1-D Gray-Scott equations. 685 686 The ODEs are derived using the method of lines. 687 """ 688 # The vectors u and v are interleaved in y. We define 689 # views of u and v by slicing y. 690 u = y[::2] 691 v = y[1::2] 692 693 # dydt is the return value of this function. 694 dydt = np.empty_like(y) 695 696 # Just like u and v are views of the interleaved vectors 697 # in y, dudt and dvdt are views of the interleaved output 698 # vectors in dydt. 699 dudt = dydt[::2] 700 dvdt = dydt[1::2] 701 702 # Compute du/dt and dv/dt. The end points and the interior points 703 # are handled separately. 704 dudt[0] = G(u[0], v[0], f, k) + Du * (-2.0*u[0] + 2.0*u[1]) / dx**2 705 dudt[1:-1] = G(u[1:-1], v[1:-1], f, k) + Du * np.diff(u,2) / dx**2 706 dudt[-1] = G(u[-1], v[-1], f, k) + Du * (- 2.0*u[-1] + 2.0*u[-2]) / dx**2 707 dvdt[0] = H(u[0], v[0], f, k) + Dv * (-2.0*v[0] + 2.0*v[1]) / dx**2 708 dvdt[1:-1] = H(u[1:-1], v[1:-1], f, k) + Dv * np.diff(v,2) / dx**2 709 dvdt[-1] = H(u[-1], v[-1], f, k) + Dv * (-2.0*v[-1] + 2.0*v[-2]) / dx**2 710 711 return dydt 712 713We won't implement a function to compute the Jacobian, but we will tell 714`odeint` that the Jacobian matrix is banded. This allows the underlying 715solver (LSODA) to avoid computing values that it knows are zero. For a large 716system, this improves the performance significantly, as demonstrated in the 717following ipython session. 718 719First, we define the required inputs:: 720 721 In [30]: rng = np.random.default_rng() 722 723 In [31]: y0 = rng.standard_normal(5000) 724 725 In [32]: t = np.linspace(0, 50, 11) 726 727 In [33]: f = 0.024 728 729 In [34]: k = 0.055 730 731 In [35]: Du = 0.01 732 733 In [36]: Dv = 0.005 734 735 In [37]: dx = 0.025 736 737Time the computation without taking advantage of the banded structure 738of the Jacobian matrix:: 739 740 In [38]: %timeit sola = odeint(grayscott1d, y0, t, args=(f, k, Du, Dv, dx)) 741 1 loop, best of 3: 25.2 s per loop 742 743Now set ``ml=2`` and ``mu=2``, so `odeint` knows that the Jacobian matrix 744is banded:: 745 746 In [39]: %timeit solb = odeint(grayscott1d, y0, t, args=(f, k, Du, Dv, dx), ml=2, mu=2) 747 10 loops, best of 3: 191 ms per loop 748 749That is quite a bit faster! 750 751Let's ensure that they have computed the same result:: 752 753 In [41]: np.allclose(sola, solb) 754 Out[41]: True 755 756References 757~~~~~~~~~~ 758 759.. [WPR] https://en.wikipedia.org/wiki/Romberg's_method 760 761.. [MOL] https://en.wikipedia.org/wiki/Method_of_lines 762