Integration (scipy.integrate)#
Das Unterpaket scipy.integrate bietet mehrere Integrationstechniken, einschließlich eines Integrators für gewöhnliche Differentialgleichungen. Eine Übersicht über das Modul finden Sie über den Befehl 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.
Allgemeine Integration (quad)#
Die Funktion quad dient zur Integration einer Funktion einer Variablen zwischen zwei Punkten. Die Punkte können \(\pm\infty\) (\(\pm\) inf) sein, um unendliche Grenzen anzuzeigen. Angenommen, Sie möchten beispielsweise eine Bessel-Funktion jv(2.5, x) über das Intervall \([0, 4.5].\) integrieren.
Dies könnte mit quad berechnet werden
>>> 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
Das erste Argument von quad ist ein "aufrufbares" Python-Objekt (d.h. eine Funktion, Methode oder Instanz einer Klasse). Beachten Sie in diesem Fall die Verwendung einer Lambda-Funktion als Argument. Die nächsten beiden Argumente sind die Integrationsgrenzen. Der Rückgabewert ist ein Tupel, wobei das erste Element den geschätzten Wert des Integrals und das zweite Element eine Schätzung des absoluten Integrationsfehlers enthält. Beachten Sie, dass der wahre Wert dieses Integrals in diesem Fall
wo
die Fresnel-Sinus-Integral ist. Beachten Sie, dass das numerisch berechnete Integral innerhalb von \(1.04\times10^{-11}\) vom exakten Ergebnis liegt – weit unter der angegebenen Fehlerschätzung.
Wenn die zu integrierende Funktion zusätzliche Parameter aufweist, können diese im Argument args angegeben werden. Angenommen, das folgende Integral soll berechnet werden
Dieses Integral kann mit dem folgenden Code ausgewertet werden
>>> 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)
Unendliche Eingaben sind in quad ebenfalls zulässig, indem \(\pm\) inf als eines der Argumente verwendet wird. Angenommen, Sie möchten beispielsweise einen numerischen Wert für das Exponentialintegral
erhalten (und vergessen, dass dieses Integral als special.expn(n,x) berechnet werden kann). Die Funktionalität der Funktion special.expn kann durch Definieren einer neuen Funktion vec_expint auf Basis der Routine quad nachgebildet werden
>>> 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])
Die integrierte Funktion kann sogar das quad-Argument verwenden (obwohl die Fehlerschranke den Fehler aufgrund von möglichem numerischem Fehler im Integranden aus der Verwendung von quad unterschätzen kann). Das Integral ist in diesem Fall
>>> 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
Das letzte Beispiel zeigt, dass mehrere Integrale durch wiederholte Aufrufe von quad behandelt werden können.
Warnung
Numerische Integrationsalgorithmen sampeln den Integranden an einer endlichen Anzahl von Punkten. Folglich können sie keine genauen Ergebnisse (oder Genauigkeitsschätzungen) für beliebige Integranden und Integrationsgrenzen garantieren. Betrachten Sie beispielsweise das Gaußsche Integral
>>> 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
Da der Integrand außer in der Nähe des Ursprungs fast Null ist, würden wir erwarten, dass große, aber endliche Integrationsgrenzen dasselbe Ergebnis liefern. Jedoch
>>> integrate.quad(gaussian, -10000, 10000)
(1.975190562208035e-203, 0.0)
Dies geschieht, weil die adaptive Quadraturroutine, die in quad implementiert ist, obwohl sie wie vorgesehen funktioniert, den kleinen, wichtigen Teil der Funktion innerhalb eines so großen, endlichen Intervalls nicht bemerkt. Für beste Ergebnisse sollten Integrationsgrenzen verwendet werden, die den wichtigen Teil des Integranden eng umschließen.
>>> integrate.quad(gaussian, -15, 15)
(1.772453850905516, 8.476526631214648e-11)
Integranden mit mehreren wichtigen Regionen können bei Bedarf in Stücke zerlegt werden.
Allgemeine Mehrfachintegration (dblquad, tplquad, nquad)#
Die Mechanik für Doppel- und Dreifachintegration sind in den Funktionen dblquad und tplquad gekapselt. Diese Funktionen nehmen die zu integrierende Funktion und vier bzw. sechs Argumente entgegen. Die Grenzen aller inneren Integrale müssen als Funktionen definiert werden.
Ein Beispiel für die Verwendung der Doppelintegration zur Berechnung mehrerer Werte von \(I_{n}\) ist unten gezeigt
>>> 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)
Als Beispiel für nicht konstante Grenzen betrachten wir das Integral
Dieses Integral kann mit dem folgenden Ausdruck ausgewertet werden (beachten Sie die Verwendung von nicht-konstanten Lambda-Funktionen für die obere Grenze des inneren Integrals)
>>> 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)
Für n-fache Integration bietet SciPy die Funktion nquad. Die Integrationsgrenzen sind ein iterierbares Objekt: entweder eine Liste konstanter Grenzen oder eine Liste von Funktionen für die nicht-konstanten Integrationsgrenzen. Die Reihenfolge der Integration (und damit der Grenzen) erfolgt vom innersten Integral zum äußersten.
Das Integral von oben
kann berechnet werden als
>>> 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)
Beachten Sie, dass die Reihenfolge der Argumente für f mit der Reihenfolge der Integrationsgrenzen übereinstimmen muss; d.h. das innere Integral bezüglich \(t\) liegt im Intervall \([1, \infty]\) und das äußere Integral bezüglich \(x\) liegt im Intervall \([0, \infty]\).
Nicht-konstante Integrationsgrenzen können auf ähnliche Weise behandelt werden; das Beispiel von oben
kann mittels ausgewertet werden
>>> 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)
was dasselbe Ergebnis wie zuvor liefert.
Gaußsche Quadratur#
fixed_quad führt eine Gaußsche Quadratur fester Ordnung über ein festes Intervall durch. Diese Funktion verwendet die Sammlung orthogonaler Polynome, die von scipy.special bereitgestellt wird und die Wurzeln und Quadraturgewichte einer großen Vielfalt orthogonaler Polynome berechnen kann (die Polynome selbst sind als Spezialfunktionen verfügbar, die Instanzen der Polynomklasse zurückgeben – z.B. special.legendre).
Integration mithilfe von Stichproben#
Wenn die Stichproben gleichmäßig verteilt sind und die Anzahl der verfügbaren Stichproben \(2^{k}+1\) für eine ganze Zahl \(k\) beträgt, dann kann die Romberg-Integration (romb) verwendet werden, um hochpräzise Schätzungen des Integrals aus den verfügbaren Stichproben zu erhalten. Die Romberg-Integration verwendet die Trapezregel bei Schrittweiten, die sich mit einer Zweierpotenz verhalten, und führt dann eine Richardson-Extrapolation dieser Schätzungen durch, um das Integral mit höherer Genauigkeit zu approximieren.
Bei beliebig verteilten Stichproben stehen die beiden Funktionen trapezoid und simpson zur Verfügung. Sie verwenden Newton-Coates-Formeln der Ordnung 1 bzw. 2 zur Durchführung der Integration. Die Trapezregel approximiert die Funktion als gerade Linie zwischen benachbarten Punkten, während die Simpson-Regel die Funktion zwischen drei benachbarten Punkten als Parabel approximiert.
Bei einer ungeraden Anzahl von Stichproben, die gleichmäßig verteilt sind, ist die Simpson-Regel exakt, wenn die Funktion ein Polynom 3. Grades oder niedriger ist. Wenn die Stichproben nicht gleichmäßig verteilt sind, ist das Ergebnis nur dann exakt, wenn die Funktion ein Polynom 2. Grades oder niedriger ist.
>>> 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
Dies entspricht genau
während die Integration der zweiten Funktion
>>> y2 = f2(x)
>>> I2 = integrate.simpson(y2, x=x)
>>> print(I2)
61.5
nicht entspricht
da der Grad des Polynoms in f2 größer als zwei ist.
Schnellere Integration mithilfe von Low-Level-Callback-Funktionen#
Ein Benutzer, der reduzierte Integrationszeiten wünscht, kann einen C-Funktionszeiger über scipy.LowLevelCallable an quad, dblquad, tplquad oder nquad übergeben und es wird integriert und ein Ergebnis in Python zurückgegeben. Die Leistungssteigerung ergibt sich aus zwei Faktoren. Die primäre Verbesserung ist die schnellere Funktionsauswertung, die durch die Kompilierung der Funktion selbst bereitgestellt wird. Zusätzlich erhalten wir eine Beschleunigung durch die Entfernung von Funktionsaufrufen zwischen C und Python in quad. Diese Methode kann eine Geschwindigkeitssteigerung von ca. 2x für triviale Funktionen wie Sinus erzielen, aber für komplexere Funktionen eine deutlich spürbarere Verbesserung (10x+) erzielen. Dieses Feature richtet sich daher an Benutzer mit numerisch intensiven Integrationen, die bereit sind, ein wenig C zu schreiben, um die Rechenzeit erheblich zu reduzieren.
Der Ansatz kann beispielsweise über ctypes in wenigen einfachen Schritten verwendet werden
1.) Schreiben Sie eine Integrandenfunktion in C mit der Funktionssignatur double f(int n, double *x, void *user_data), wobei x ein Array mit dem Punkt ist, an dem die Funktion f ausgewertet wird, und user_data für beliebige zusätzliche Daten steht, die Sie bereitstellen möchten.
/* 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.) Kompilieren Sie diese Datei nun zu einer Shared/Dynamic Library (eine schnelle Suche wird dabei helfen, da sie vom Betriebssystem abhängt). Der Benutzer muss alle verwendeten mathematischen Bibliotheken usw. verknüpfen. Unter Linux sieht dies so aus
$ gcc -shared -fPIC -o testlib.so testlib.c
Die Ausgabebibliothek wird als testlib.so referenziert, kann aber eine andere Dateiendung haben. Eine Bibliothek wurde nun erstellt, die mit ctypes in Python geladen werden kann.
3.) Laden Sie die Shared Library mit ctypes in Python und setzen Sie restypes und argtypes - dies ermöglicht SciPy, die Funktion korrekt zu interpretieren
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)
Das letzte void *user_data in der Funktion ist optional und kann weggelassen werden (sowohl in der C-Funktion als auch in den ctypes argtypes), wenn es nicht benötigt wird. Beachten Sie, dass die Koordinaten als Array von Doubles und nicht als separates Argument übergeben werden.
4.) Integrieren Sie nun die Bibliotheksfunktion wie gewohnt, hier mit nquad
>>> integrate.nquad(func, [[0, 10], [-10, 0], [-1, 1]])
(1200.0, 1.1102230246251565e-11)
Das Python-Tupel wird wie erwartet in reduzierter Zeit zurückgegeben. Alle optionalen Parameter können mit dieser Methode verwendet werden, einschließlich der Angabe von Singularitäten, unendlichen Grenzen usw.
Gewöhnliche Differentialgleichungen (solve_ivp)#
Die Integration eines Systems gewöhnlicher Differentialgleichungen (ODEs) mit Anfangsbedingungen ist ein weiteres nützliches Beispiel. Die Funktion solve_ivp ist in SciPy zur Integration einer Differentialgleichung erster Ordnung in Vektorform verfügbar
mit den Anfangsbedingungen \(\mathbf{y}\left(0\right)=\mathbf{y}_{0}\), wobei \(\mathbf{y}\) ein Vektor der Länge \(N\) ist und \(\mathbf{f}\) eine Abbildung von \(\mathbb{R}^{N}\) nach \(\mathbb{R}^{N}.\) Eine gewöhnliche Differentialgleichung höherer Ordnung kann immer auf eine Differentialgleichung dieses Typs reduziert werden, indem Zwischenableitungen in den \(\mathbf{y}\)-Vektor eingeführt werden.
Angenommen, wir möchten beispielsweise die Lösung der folgenden Differentialgleichung zweiter Ordnung finden
mit den Anfangsbedingungen \(w\left(0\right)=\frac{1}{\sqrt[3]{3^{2}}\Gamma\left(\frac{2}{3}\right)}\) und \(\left.\frac{dw}{dz}\right|_{z=0}=-\frac{1}{\sqrt[3]{3}\Gamma\left(\frac{1}{3}\right)}.\) Es ist bekannt, dass die Lösung dieser Differentialgleichung mit diesen Randbedingungen die Airy-Funktion ist
was eine Möglichkeit bietet, den Integrator mithilfe von special.airy zu überprüfen.
Konvertieren Sie diese ODE zunächst in Standardform, indem Sie \(\mathbf{y}=\left[\frac{dw}{dz},w\right]\) und \(t=z\) setzen. Somit wird die Differentialgleichung zu
Mit anderen Worten,
Zur Erinnerung: Wenn \(\mathbf{A}\left(t\right)\) mit \(\int_{0}^{t}\mathbf{A}\left(\tau\right)\, d\tau\) unter Matrizenmultiplikation pendelt, dann hat diese lineare Differentialgleichung eine exakte Lösung mithilfe der Matrixexponentialfunktion
jedoch pendeln \(\mathbf{A}\left(t\right)\) und sein Integral in diesem Fall nicht.
Diese Differentialgleichung kann mit der Funktion solve_ivp gelöst werden. Sie benötigt als Eingabeargumente die Ableitung, *fprime*, den Zeitbereich [t_start, t_end] und den Anfangsbedingungsvektor, *y0*, und gibt ein Objekt zurück, dessen Feld *y* ein Array mit aufeinanderfolgenden Lösungs Werten als Spalten ist. Die Anfangsbedingungen sind daher in der ersten Ausgabespalte gegeben.
>>> 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. ]
Wie zu sehen ist, bestimmt solve_ivp seine Zeitschritte automatisch, wenn nicht anders angegeben. Um die Lösung von solve_ivp mit der *airy*-Funktion zu vergleichen, wird der von solve_ivp erzeugte Zeitvektor an die *airy*-Funktion übergeben.
>>> 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]
Die Lösung von solve_ivp mit seinen Standardparametern zeigt eine große Abweichung zur Airy-Funktion. Um diese Abweichung zu minimieren, können relative und absolute Toleranzen verwendet werden.
>>> 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]
Um benutzerdefinierte Zeitpunkte für die Lösung von solve_ivp anzugeben, bietet solve_ivp zwei Möglichkeiten, die auch komplementär genutzt werden können. Durch Übergabe der Option *t_eval* an den Funktionsaufruf gibt solve_ivp die Lösungen dieser Zeitpunkte von *t_eval* in seiner Ausgabe zurück.
>>> import numpy as np
>>> t = np.linspace(0, 4, 100)
>>> sol3 = solve_ivp(func, t_span, y0, t_eval=t)
Wenn die Jacobi-Matrix der Funktion bekannt ist, kann sie an solve_ivp übergeben werden, um bessere Ergebnisse zu erzielen. Beachten Sie jedoch, dass die Standardintegrationsmethode RK45 keine Jacobi-Matrizen unterstützt und daher eine andere Integrationsmethode gewählt werden muss. Eine der Integrationsmethoden, die eine Jacobi-Matrix unterstützt, ist zum Beispiel die Methode Radau des folgenden Beispiels.
>>> def gradient(t, y):
... return [[0,t], [1,0]]
>>> sol4 = solve_ivp(func, t_span, y0, method='Radau', jac=gradient)
Lösen eines Systems mit einer bandförmigen Jacobi-Matrix#
odeint kann mitgeteilt werden, dass die Jacobi-Matrix *bandförmig* ist. Für ein großes System von Differentialgleichungen, die bekanntermaßen steif sind, kann dies die Leistung erheblich verbessern.
Als Beispiel lösen wir die 1D Gray-Scott partielle Differentialgleichungen mit der Methode der Linien [MOL]. Die Gray-Scott-Gleichungen für die Funktionen \(u(x, t)\) und \(v(x, t)\) auf dem Intervall \(x \in [0, L]\) sind
wobei \(D_u\) und \(D_v\) die Diffusionskoeffizienten der Komponenten \(u\) und \(v\) sind und \(f\) und \(k\) Konstanten sind. (Weitere Informationen über das System finden Sie unter http://groups.csail.mit.edu/mac/projects/amorphous/GrayScott/)
Wir nehmen Neumann-Randbedingungen (d.h. "kein Fluss") an
Um die Methode der Linien anzuwenden, diskretisieren wir die Variable \(x\), indem wir das gleichmäßig verteilte Gitter von \(N\) Punkten \(\left\{x_0, x_1, \ldots, x_{N-1}\right\}\) definieren, mit \(x_0 = 0\) und \(x_{N-1} = L\). Wir definieren \(u_j(t) \equiv u(x_k, t)\) und \(v_j(t) \equiv v(x_k, t)\) und ersetzen die \(x\)-Ableitungen durch endliche Differenzen. Das heißt,
Wir erhalten dann ein System von \(2N\) gewöhnlichen Differentialgleichungen
Der Einfachheit halber wurden die Argumente \((t)\) weggelassen.
Um die Randbedingungen durchzusetzen, führen wir "Ghost"-Punkte \(-1\) und \(x_N\) ein und definieren \(u_{-1}(t) \equiv u_1(t)\), \(u_N(t) \equiv u_{N-2}(t)\); \(v_{-1}(t)\) und \(v_N(t)\) werden analog definiert.
Dann
und
Unser vollständiges System von \(2N\) gewöhnlichen Differentialgleichungen ist (1) für \(k = 1, 2, \ldots, N-2\), zusammen mit (2) und (3).
Wir können nun mit der Implementierung dieses Systems in Code beginnen. Wir müssen \(\{u_k\}\) und \(\{v_k\}\) zu einem einzigen Vektor der Länge \(2N\) zusammenfassen. Die beiden offensichtlichen Wahlmöglichkeiten sind \(\{u_0, u_1, \ldots, u_{N-1}, v_0, v_1, \ldots, v_{N-1}\}\) und \(\{u_0, v_0, u_1, v_1, \ldots, u_{N-1}, v_{N-1}\}\). Mathematisch spielt es keine Rolle, aber die Wahl beeinflusst, wie effizient odeint das System lösen kann. Der Grund liegt darin, wie die Reihenfolge das Muster der Nicht-Null-Elemente der Jacobi-Matrix beeinflusst.
Wenn die Variablen als \(\{u_0, u_1, \ldots, u_{N-1}, v_0, v_1, \ldots, v_{N-1}\}\) geordnet sind, ist das Muster der Nicht-Null-Elemente der Jacobi-Matrix
Das Jacobi-Muster mit verschachtelten Variablen als \(\{u_0, v_0, u_1, v_1, \ldots, u_{N-1}, v_{N-1}\}\) ist
In beiden Fällen gibt es nur fünf nicht-triviale Diagonalen, aber wenn die Variablen verschachtelt sind, ist die Bandbreite viel kleiner. Das heißt, die Hauptdiagonale und die beiden Diagonalen unmittelbar darüber und die beiden unmittelbar darunter sind die Nicht-Null-Diagonalen. Dies ist wichtig, da die Eingaben mu und ml von odeint die oberen und unteren Bandbreiten der Jacobi-Matrix sind. Wenn die Variablen verschachtelt sind, sind mu und ml 2. Wenn die Variablen mit \(\{v_k\}\) nach \(\{u_k\}\) gestapelt sind, sind die oberen und unteren Bandbreiten \(N\).
Mit dieser Entscheidung können wir die Funktion schreiben, die das System von Differentialgleichungen implementiert.
Zuerst definieren wir die Funktionen für die Quell- und Reaktionsterme des Systems
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
Als Nächstes definieren wir die Funktion, die die rechte Seite des Systems von Differentialgleichungen berechnet
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
Wir werden keine Funktion zur Berechnung des Jacobians implementieren, aber wir werden odeint mitteilen, dass die Jacobi-Matrix bandförmig ist. Dies ermöglicht es dem zugrundeliegenden Löser (LSODA), die Berechnung von Werten zu vermeiden, von denen er weiß, dass sie null sind. Für ein großes System verbessert dies die Leistung erheblich, wie die folgende ipython-Sitzung zeigt.
Zuerst definieren wir die erforderlichen Eingaben
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
Zeitliche Berechnung, ohne die bandförmige Struktur der Jacobi-Matrix zu nutzen
In [38]: %timeit sola = odeint(grayscott1d, y0, t, args=(f, k, Du, Dv, dx))
1 loop, best of 3: 25.2 s per loop
Setzen Sie nun ml=2 und mu=2, damit odeint weiß, dass die Jacobi-Matrix bandförmig ist
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
Das ist ziemlich viel schneller!
Stellen wir sicher, dass sie das gleiche Ergebnis berechnet haben
In [41]: np.allclose(sola, solb)
Out[41]: True