Интегрирование (scipy.integrate)#
The scipy.integrate подпакет предоставляет несколько методов интегрирования, включая интегратор обыкновенных дифференциальных уравнений. Обзор модуля предоставляется командой help:
>>> 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.
Общее интегрирование (quad)#
Функция quad предоставляется для интегрирования функции одной
переменной между двумя точками. Точки могут быть \(\pm\infty\)
(\(\pm\) inf) для обозначения бесконечных пределов. Например, предположим, вы хотите проинтегрировать функцию Бесселя jv(2.5, x) вдоль интервала \([0, 4.5].\)
Это может быть вычислено с использованием 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
Первый аргумент quad — это вызываемый объект Python (например, функция, метод или экземпляр класса). Обратите внимание на использование лямбда-функции в данном случае в качестве аргумента. Следующие два аргумента — это пределы интегрирования. Возвращаемое значение — это кортеж, где первый элемент содержит оценённое значение интеграла, а второй — оценку абсолютной ошибки интегрирования. Обратите внимание, что в данном случае истинное значение этого интеграла
где
является интегралом Френеля от синуса. Заметим, что численно вычисленный интеграл находится в пределах \(1.04\times10^{-11}\) точного результата — значительно ниже заявленной оценки ошибки.
Если интегрируемая функция принимает дополнительные параметры, они могут быть предоставлены в args аргумент. Предположим, что требуется вычислить следующий интеграл:
Этот интеграл можно вычислить с помощью следующего кода:
>>> 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)
Бесконечные входные данные также допустимы в quad используя \(\pm\)
inf в качестве одного из аргументов. Например, предположим, что числовое значение для экспоненциального интеграла:
желательно (и тот факт, что этот интеграл может быть вычислен как
special.expn(n,x) забывается). Функциональность функции
special.expn может быть воспроизведено путем определения новой функции
vec_expint на основе процедуры 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])
Функция, которая интегрируется, может даже использовать аргумент quad (хотя оценка ошибки может недооценивать ошибку из-за возможной численной ошибки в подынтегральном выражении от использования quad ). Интеграл в этом случае равен
>>> 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
Этот последний пример показывает, что множественное интегрирование может быть выполнено с использованием
повторных вызовов quad.
Предупреждение
Алгоритмы численного интегрирования берут выборки подынтегральной функции в конечном числе точек. Следовательно, они не могут гарантировать точные результаты (или оценки точности) для произвольных подынтегральных функций и пределов интегрирования. Рассмотрим, например, интеграл Гаусса:
>>> 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
Поскольку подынтегральная функция почти нулевая, за исключением области вблизи начала координат, можно ожидать, что большие, но конечные пределы интегрирования дадут тот же результат. Однако:
>>> integrate.quad(gaussian, -10000, 10000)
(1.975190562208035e-203, 0.0)
Это происходит потому, что адаптивная квадратурная процедура, реализованная в quad,
хотя и работает как задумано, не замечает небольшую, важную часть функции
в таком большом конечном интервале. Для лучших результатов рассмотрите использование пределов
интегрирования, которые тесно окружают важную часть подынтегрального выражения.
>>> integrate.quad(gaussian, -15, 15)
(1.772453850905516, 8.476526631214648e-11)
Подынтегральные выражения с несколькими важными областями могут быть разбиты на части по мере необходимости.
Общее многократное интегрирование (dblquad, tplquad, nquad)#
Механика двойного и тройного интегрирования была обёрнута в функции dblquad и tplquad. Эти функции принимают функцию для интегрирования и четыре или шесть аргументов соответственно. Пределы всех внутренних интегралов должны быть определены как функции.
Пример использования двойного интегрирования для вычисления нескольких значений \(I_{n}\) показано ниже:
>>> 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)
В качестве примера для непостоянных пределов рассмотрим интеграл
Этот интеграл может быть вычислен с использованием выражения ниже (Обратите внимание на использование непостоянных лямбда-функций для верхнего предела внутреннего интеграла):
>>> 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)
Для n-кратного интегрирования scipy предоставляет функцию nquad.
Границы интегрирования являются итерируемым объектом: либо список постоянных границ,
либо список функций для непостоянных границ интегрирования. Порядок
интегрирования (и, следовательно, границ) — от самого внутреннего интеграла к
самому внешнему.
Интеграл сверху
может быть вычислено как
>>> 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)
Обратите внимание, что порядок аргументов для f должен соответствовать порядку границ интегрирования; т.е., внутренний интеграл по \(t\) находится на интервале \([1, \infty]\) и внешний интеграл по \(x\) находится в интервале \([0, \infty]\).
Непостоянные границы интегрирования могут быть обработаны аналогичным образом; пример из вышеприведённого
может быть вычислено с помощью
>>> 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)
что даёт тот же результат, что и ранее.
Квадратура Гаусса#
fixed_quad выполняет квадратуру Гаусса фиксированного порядка на фиксированном интервале. Эта функция использует набор ортогональных полиномов, предоставленных scipy.special,
который может вычислять корни и веса квадратуры для большого разнообразия ортогональных
полиномов (сами полиномы доступны как специальные функции, возвращающие
экземпляры класса полиномов — например,
special.legendre).
Интегрирование с использованием выборок#
Если выборки равноотстоящие и количество доступных выборок \(2^{k}+1\) для некоторого целого числа \(k\), затем Ромберг romb
интегрирование может использоваться для получения высокоточных оценок
интеграла с использованием доступных выборок. Интегрирование Ромберга использует
правило трапеций с шагами, связанными степенью двойки, а затем
выполняет экстраполяцию Ричардсона по этим оценкам для приближения
интеграла с более высокой степенью точности.
В случае произвольно расположенных выборок можно использовать две функции trapezoid
и simpson доступны. Они используют формулы Ньютона-Котеса порядка 1 и 2 соответственно для выполнения интегрирования. Трапециевидное правило аппроксимирует функцию прямой линией между соседними точками, в то время как правило Симпсона аппроксимирует функцию между тремя соседними точками как параболу.
Для нечётного количества равноотстоящих выборок правило Симпсона точно, если функция является полиномом порядка 3 или меньше. Если выборки не равноотстоящие, результат точен только если функция является полиномом порядка 2 или меньше.
>>> 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=x)
>>> print(I1)
21.0
Это точно соответствует
в то время как интегрирование второй функции
>>> y2 = f2(x)
>>> I2 = integrate.simpson(y2, x=x)
>>> print(I2)
61.5
не соответствует
потому что порядок полинома в f2 больше двух.
Более быстрая интеграция с использованием низкоуровневых функций обратного вызова#
Пользователь, желающий сократить время интегрирования, может передать указатель на C-функцию
через scipy.LowLevelCallable to quad, dblquad,
tplquad или nquad и он будет проинтегрирован и вернет результат в Python. Увеличение производительности здесь возникает из двух факторов. Основное улучшение — более быстрая оценка функции, которая обеспечивается компиляцией самой функции. Дополнительно у нас есть ускорение, обеспечиваемое удалением вызовов функций между C и Python в
quad. Этот метод может обеспечить ускорение примерно в 2 раза для тривиальных функций, таких как синус, но может дать гораздо более заметное улучшение (10+ раз) для более сложных функций. Эта функция ориентирована на пользователей с численно интенсивными интеграциями, готовых написать немного кода на C для значительного сокращения времени вычислений.
Подход может быть использован, например, через ctypes в несколько простых шагов:
1.) Написать подынтегральную функцию на C с сигнатурой функции
double f(int n, double *x, void *user_data), где x — это массив, содержащий точку, в которой вычисляется функция f, и user_data
для произвольных дополнительных данных, которые вы хотите предоставить.
/* 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.) Теперь скомпилируйте этот файл в общую/динамическую библиотеку (быстрый поиск поможет с этим, так как это зависит от ОС). Пользователь должен связать любые математические библиотеки и т.д., которые используются. В Linux это выглядит так:
$ gcc -shared -fPIC -o testlib.so testlib.c
Выходная библиотека будет называться testlib.so, но он может иметь другое расширение файла. Теперь создана библиотека, которую можно загрузить в Python с помощью ctypes.
3.) Загрузить общую библиотеку в Python с помощью ctypes и установить restypes и
argtypes - это позволяет SciPy правильно интерпретировать функцию:
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)
Последний void *user_data в функции является необязательным и может быть опущено (как в C-функции, так и в ctypes argtypes), если не требуется. Обратите внимание, что координаты передаются в виде массива double, а не отдельного аргумента.
4.) Теперь интегрируйте библиотечную функцию как обычно, здесь используя nquad:
>>> integrate.nquad(func, [[0, 10], [-10, 0], [-1, 1]])
(1200.0, 1.1102230246251565e-11)
Кортеж Python возвращается, как и ожидалось, за меньшее время. Все дополнительные параметры могут быть использованы с этим методом, включая указание особенностей, бесконечных границ и т.д.
Обыкновенные дифференциальные уравнения (solve_ivp)#
Интегрирование системы обыкновенных дифференциальных уравнений (ОДУ) при заданных
начальных условиях — ещё один полезный пример. Функция
solve_ivp доступен в SciPy для интегрирования векторного дифференциального уравнения первого порядка:
заданные начальные условия \(\mathbf{y}\left(0\right)=\mathbf{y}_{0}\), где \(\mathbf{y}\) имеет длину \(N\) вектор и \(\mathbf{f}\) является отображением из \(\mathbb{R}^{N}\) to \(\mathbb{R}^{N}.\) Обыкновенное дифференциальное уравнение высшего порядка всегда может быть сведено к дифференциальному уравнению этого типа путём введения промежуточных производных в \(\mathbf{y}\) вектор.
Например, предположим, что требуется найти решение следующего дифференциального уравнения второго порядка:
с начальными условиями \(w\left(0\right)=\frac{1}{\sqrt[3]{3^{2}}\Gamma\left(\frac{2}{3}\right)}\) и \(\left.\frac{dw}{dz}\right|_{z=0}=-\frac{1}{\sqrt[3]{3}\Gamma\left(\frac{1}{3}\right)}.\) Известно, что решение этого дифференциального уравнения с этими граничными условиями является функцией Эйри
что дает возможность проверить интегратор с помощью special.airy.
Сначала преобразуйте это ОДУ в стандартную форму, установив \(\mathbf{y}=\left[\frac{dw}{dz},w\right]\) и \(t=z\). Таким образом, дифференциальное уравнение становится
Другими словами,
Как интересное напоминание, если \(\mathbf{A}\left(t\right)\) коммутирует с \(\int_{0}^{t}\mathbf{A}\left(\tau\right)\, d\tau\) при матричном умножении, тогда это линейное дифференциальное уравнение имеет точное решение с использованием матричной экспоненты:
Однако в этом случае, \(\mathbf{A}\left(t\right)\) и его интеграл не коммутируют.
Это дифференциальное уравнение может быть решено с помощью функции solve_ivp.
Требует производную, fprime, временной интервал [t_start, t_end]
и вектор начальных условий, y0, в качестве входных аргументов и возвращает
объект, чьи y поле представляет собой массив с последовательными значениями решения в виде
столбцов. Начальные условия, таким образом, заданы в первом выходном столбце.
>>> 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. ]
Как можно видеть solve_ivp определяет свои временные шаги автоматически, если не
указано иное. Чтобы сравнить решение solve_ivp с airy
функция временного вектора, созданного solve_ivp передается в airy функция.
>>> 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]
Решение solve_ivp со стандартными параметрами показывает большое отклонение от функции Эйри. Чтобы минимизировать это отклонение, можно использовать относительные и абсолютные допуски.
>>> 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]
Для указания пользовательских временных точек для решения solve_ivp, solve_ivp
предлагает две возможности, которые также могут использоваться взаимодополняюще. Передавая t_eval
опция для вызова функции solve_ivp возвращает решения в эти моменты времени t_eval в его выводе.
>>> import numpy as np
>>> t = np.linspace(0, 4, 100)
>>> sol3 = solve_ivp(func, t_span, y0, t_eval=t)
Если матрица Якоби функции известна, её можно передать в solve_ivp
для достижения лучших результатов. Однако учтите, что метод интегрирования по умолчанию
RK45 не поддерживает матрицы Якоби, поэтому необходимо выбрать другой метод
интегрирования. Один из методов интегрирования, поддерживающих матрицу Якоби, — это, например, Radau метод следующего примера.
>>> def gradient(t, y):
... return [[0,t], [1,0]]
>>> sol4 = solve_ivp(func, t_span, y0, method='Radau', jac=gradient)
Решение системы с ленточной матрицей Якоби#
odeint можно указать, что якобиан является ленточная. Для большой системы дифференциальных уравнений, которые известны как жёсткие, это может значительно улучшить производительность.
В качестве примера мы решим одномерные уравнения Грея-Скотта в частных производных методом линий [MOL]. Уравнения Грея-Скотта для функций \(u(x, t)\) и \(v(x, t)\) на интервале \(x \in [0, L]\) являются
где \(D_u\) и \(D_v\) are the diffusion coefficients of the components \(u\) и \(v\), соответственно, и \(f\) и \(k\) являются константами. (Для получения дополнительной информации о системе, см. http://groups.csail.mit.edu/mac/projects/amorphous/GrayScott/)
Мы предположим граничные условия Неймана (т.е., "без потока"):
Чтобы применить метод линий, мы дискретизируем \(x\) переменной, определив равномерно распределённую сетку \(N\) точки \(\left\{x_0, x_1, \ldots, x_{N-1}\right\}\), с \(x_0 = 0\) и \(x_{N-1} = L\). Мы определяем \(u_j(t) \equiv u(x_k, t)\) и \(v_j(t) \equiv v(x_k, t)\), и заменить \(x\) производные с конечными разностями. То есть,
Тогда у нас есть система \(2N\) обыкновенные дифференциальные уравнения:
Для удобства, \((t)\) аргументы были удалены.
Чтобы обеспечить граничные условия, мы вводим «призрачные» точки \(x_{-1}\) и \(x_N\), и определить \(u_{-1}(t) \equiv u_1(t)\), \(u_N(t) \equiv u_{N-2}(t)\); \(v_{-1}(t)\) и \(v_N(t)\) определяются аналогично.
Затем
и
Наша полная система \(2N\) обыкновенных дифференциальных уравнений есть (1) для \(k = 1, 2, \ldots, N-2\), вместе с (2) и (3).
Теперь мы можем начать реализацию этой системы в коде. Мы должны объединить
\(\{u_k\}\) и \(\{v_k\}\) в единый вектор длины \(2N\).
Два очевидных варианта —
\(\{u_0, u_1, \ldots, u_{N-1}, v_0, v_1, \ldots, v_{N-1}\}\)
и
\(\{u_0, v_0, u_1, v_1, \ldots, u_{N-1}, v_{N-1}\}\).
Математически это не имеет значения, но выбор влияет на то, насколько эффективно odeint может решить систему. Причина в том, как
порядок влияет на структуру ненулевых элементов матрицы Якоби.
Когда переменные упорядочены как \(\{u_0, u_1, \ldots, u_{N-1}, v_0, v_1, \ldots, v_{N-1}\}\), шаблон ненулевых элементов матрицы Якоби
Шаблон Якобиана с перемешанными переменными как \(\{u_0, v_0, u_1, v_1, \ldots, u_{N-1}, v_{N-1}\}\) является
В обоих случаях есть только пять нетривиальных диагоналей, но
когда переменные перемежаются, ширина полосы значительно
меньше.
То есть главная диагональ и две диагонали непосредственно
над ней и две непосредственно под главной диагональю
являются ненулевыми диагоналями.
Это важно, потому что входные данные mu и ml
of odeint это верхняя и нижняя полосы пропускания матрицы Якоби. Когда переменные перемежаются,
mu и ml являются 2. Когда переменные складываются
с \(\{v_k\}\) следующий \(\{u_k\}\), верхняя
и нижняя полосы пропускания равны \(N\).
После принятия этого решения мы можем написать функцию, реализующую систему дифференциальных уравнений.
Сначала определим функции для источника и членов реакции системы:
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
Далее мы определяем функцию, которая вычисляет правую часть системы дифференциальных уравнений:
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
Мы не будем реализовывать функцию для вычисления якобиана, но мы укажем
odeint что матрица Якоби является ленточной. Это позволяет
основному решателю (LSODA) избегать вычисления значений, которые, как известно, равны нулю. Для большой
системы это значительно улучшает производительность, как показано в
следующей сессии ipython.
Сначала определим необходимые входные данные:
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
Измерить время вычисления без использования преимущества ленточной структуры матрицы Якоби:
In [38]: %timeit sola = odeint(grayscott1d, y0, t, args=(f, k, Du, Dv, dx))
1 loop, best of 3: 25.2 s per loop
Теперь установим ml=2 и mu=2, поэтому odeint знает, что матрица Якоби
является ленточной:
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
Это значительно быстрее!
Убедимся, что они вычислили одинаковый результат:
In [41]: np.allclose(sola, solb)
Out[41]: True