scipy.integrate.

cubature#

scipy.integrate.cubature(f, a, b, *, rule='gk21', rtol=1e-08, atol=0, max_subdivisions=10000, args=(), workers=1, points=None)[Quelle]#

Adaptives Kubatur einer mehrdimensionalen vektorwertigen Funktion.

Gegeben eine beliebige Integrationsregel, gibt diese Funktion eine Schätzung des Integrals bis zur angeforderten Toleranz über dem Bereich zurück, der durch die Arrays a und b definiert ist, welche die Ecken eines Hyperwürfels spezifizieren.

Die Konvergenz ist nicht für alle Integrale garantiert.

Parameter:

faufrufbar

Zu integrierende Funktion. f muss die Signatur haben

f(x : ndarray, *args) -> ndarray

f sollte Arrays x der Form

(npoints, ndim)

und Arrays der Form

(npoints, output_dim_1, ..., output_dim_n)

In diesem Fall gibt cubature Arrays der Form zurück

(output_dim_1, ..., output_dim_n)

a, barray_like

Untere und obere Integrationsgrenzen als 1D-Arrays, die die linken und rechten Endpunkte der integrierten Intervalle spezifizieren. Grenzen können unendlich sein.

rulestr, optional

Regel zur Schätzung des Integrals. Wenn ein String übergeben wird, sind die Optionen "gauss-kronrod" (21 Knoten) oder "genz-malik" (Grad 7). Wenn eine Regel wie "gauss-kronrod" für einen n-dimensionalen Integranden spezifiziert wird, wird die entsprechende kartesische Produktregel verwendet. "gk21", "gk15" werden ebenfalls zur Kompatibilität mit quad_vec unterstützt. Siehe Hinweise.

rtol, atolfloat, optional

Relative und absolute Toleranzen. Iterationen werden durchgeführt, bis der Fehler geschätzt wird, kleiner als atol + rtol * abs(est) zu sein. Hier steuert rtol die relative Genauigkeit (Anzahl korrekter Ziffern), während atol die absolute Genauigkeit (Anzahl korrekter Dezimalstellen) steuert. Um die gewünschte rtol zu erreichen, setzen Sie atol kleiner als den kleinsten Wert, der von rtol * abs(y) erwartet werden kann, damit rtol den zulässigen Fehler dominiert. Wenn atol größer als rtol * abs(y) ist, ist die Anzahl der korrekten Ziffern nicht garantiert. Umgekehrt, um die gewünschte atol zu erreichen, setzen Sie rtol so, dass rtol * abs(y) immer kleiner als atol ist. Standardwerte sind 1e-8 für rtol und 0 für atol.

max_subdivisionsint, optional

Obergrenze für die Anzahl der durchzuführenden Unterteilungen. Standard ist 10.000.

argstuple, optional

Zusätzliche positionelle Argumente, die an f übergeben werden, falls vorhanden.

workersint oder Map-ähnliches aufrufbare Objekt, optional

Wenn workers eine Ganzzahl ist, wird ein Teil der Berechnung parallel in so viele Aufgaben unterteilt (mit multiprocessing.pool.Pool). Geben Sie -1 an, um alle für den Prozess verfügbaren Kerne zu nutzen. Alternativ können Sie einen Map-ähnlichen Callable übergeben, wie z. B. multiprocessing.pool.Pool.map, um die Population parallel auszuwerten. Diese Auswertung erfolgt als workers(func, iterable).

pointslist of array_like, optional

Liste von Punkten, an denen die Auswertung von f vermieden werden soll, unter der Bedingung, dass die verwendete Regel f nicht am Rand eines Bereichs auswertet (was bei allen Genz-Malik- und Gauss-Kronrod-Regeln der Fall ist). Dies kann nützlich sein, wenn f eine Singularität an dem angegebenen Punkt hat. Dies sollte eine Liste von Array-ähnlichen Objekten sein, wobei jedes Element die Länge ndim hat. Standard ist eine leere Liste. Siehe Beispiele.

Rückgabe:

resobject

Objekt, das die Ergebnisse der Schätzung enthält. Es hat die folgenden Attribute:

estimatendarray: Schätzung des Werts des Integrals über den gesamten angegebenen Bereich.
errorndarray: Schätzung des Fehlers der Approximation über den gesamten angegebenen Bereich.
statusstr: Ob die Schätzung erfolgreich war. Kann sein: "converged", "not_converged".
subdivisionsint: Anzahl der durchgeführten Unterteilungen.
atol, rtolfloat: Angefragte Toleranzen für die Approximation.
regions: list of object: Liste von Objekten, die die Schätzungen des Integrals über kleinere Bereiche des Definitionsbereichs enthalten.

Jedes Objekt in regions hat die folgenden Attribute:

a, bndarray: Punkte, die die Ecken des Bereichs beschreiben. Wenn das ursprüngliche Integral unendliche Grenzen enthielt oder über einen Bereich integriert wurde, der durch region beschrieben wurde, dann sind a und b in den transformierten Koordinaten.
estimatendarray: Schätzung des Werts des Integrals über diesen Bereich.
errorndarray: Schätzung des Fehlers der Approximation über diesen Bereich.

Hinweise

Der Algorithmus verwendet einen ähnlichen Algorithmus wie quad_vec, der selbst auf der Implementierung der DQAG*-Algorithmen von QUADPACK basiert und eine globale Fehlerkontrolle sowie adaptive Unterteilung implementiert.

Die Quelle der Knoten und Gewichte für die Gauss-Kronrod-Quadratur finden Sie in [1], und der Algorithmus zur Berechnung der Knoten und Gewichte bei Genz-Malik-Kubatur finden Sie in [2].

Die über das Argument rule unterstützten Regeln sind derzeit:

"gauss-kronrod", 21-Knoten-Gauss-Kronrod
"genz-malik", n-Knoten-Genz-Malik

Wenn Gauss-Kronrod für einen n-dimensionalen Integranden verwendet wird, wobei n > 2 ist, wird die entsprechende kartesische Produktregel durch die Bildung des kartesischen Produkts der Knoten im 1D-Fall gefunden. Das bedeutet, dass die Anzahl der Knoten im Gauss-Kronrod-Fall exponentiell mit 21^n skaliert, was bei einer moderaten Anzahl von Dimensionen problematisch sein kann.

Genz-Malik ist in der Regel weniger genau als Gauss-Kronrod, hat aber deutlich weniger Knoten. Daher könnte in dieser Situation die Verwendung von "genz-malik" vorzuziehen sein.

Unendliche Grenzen werden mit einer geeigneten Variablentransformation behandelt. Angenommen a = [a_1, ..., a_n] und b = [b_1, ..., b_n]

Wenn \(a_i = -\infty\) und \(b_i = \infty\) ist, verwendet die i-te Integrationsvariable die Transformation \(x = \frac{1-|t|}{t}\) und \(t \in (-1, 1)\).

Wenn \(a_i \ne \pm\infty\) und \(b_i = \infty\) ist, verwendet die i-te Integrationsvariable die Transformation \(x = a_i + \frac{1-t}{t}\) und \(t \in (0, 1)\).

Wenn \(a_i = -\infty\) und \(b_i \ne \pm\infty\) ist, verwendet die i-te Integrationsvariable die Transformation \(x = b_i - \frac{1-t}{t}\) und \(t \in (0, 1)\).

Referenzen

[1]

R. Piessens, E. de Doncker, Quadpack: A Subroutine Package for Automatic Integration, files: dqk21.f, dqk15.f (1983).

[2]

A.C. Genz, A.A. Malik, Remarks on algorithm 006: An adaptive algorithm for numerical integration over an N-dimensional rectangular region, Journal of Computational and Applied Mathematics, Volume 6, Issue 4, 1980, Pages 295-302, ISSN 0377-0427 DOI:10.1016/0771-050X(80)90039-X

Beispiele

1D-Integral mit Vektor-Ausgabe:

\[\int^1_0 \mathbf f(x) \text dx\]

Wobei f(x) = x^n und n = np.arange(10) ein Vektor ist. Da keine Regel angegeben ist, wird die Standardregel "gk21" verwendet, die Gauss-Kronrod-Integration mit 21 Knoten entspricht.

>>> import numpy as np
>>> from scipy.integrate import cubature
>>> def f(x, n):
...    # Make sure x and n are broadcastable
...    return x[:, np.newaxis]**n[np.newaxis, :]
>>> res = cubature(
...     f,
...     a=[0],
...     b=[1],
...     args=(np.arange(10),),
... )
>>> res.estimate
 array([1.        , 0.5       , 0.33333333, 0.25      , 0.2       ,
        0.16666667, 0.14285714, 0.125     , 0.11111111, 0.1       ])

7D-Integral mit beliebig geformter Array-Ausgabe:

f(x) = cos(2*pi*r + alphas @ x)

für ein bestimmtes r und alphas, und das Integral wird über den Einheits-Hyperwürfel, \([0, 1]^7\), durchgeführt. Da das Integral in einer moderaten Anzahl von Dimensionen liegt, wird "genz-malik" anstelle des Standardwerts "gauss-kronrod" verwendet, um die Erstellung einer Produktregel mit \(21^7 \approx 2 \times 10^9\) Knoten zu vermeiden.

>>> import numpy as np
>>> from scipy.integrate import cubature
>>> def f(x, r, alphas):
...     # f(x) = cos(2*pi*r + alphas @ x)
...     # Need to allow r and alphas to be arbitrary shape
...     npoints, ndim = x.shape[0], x.shape[-1]
...     alphas = alphas[np.newaxis, ...]
...     x = x.reshape(npoints, *([1]*(len(alphas.shape) - 1)), ndim)
...     return np.cos(2*np.pi*r + np.sum(alphas * x, axis=-1))
>>> rng = np.random.default_rng()
>>> r, alphas = rng.random((2, 3)), rng.random((2, 3, 7))
>>> res = cubature(
...     f=f,
...     a=np.array([0, 0, 0, 0, 0, 0, 0]),
...     b=np.array([1, 1, 1, 1, 1, 1, 1]),
...     rtol=1e-5,
...     rule="genz-malik",
...     args=(r, alphas),
... )
>>> res.estimate
 array([[-0.79812452,  0.35246913, -0.52273628],
        [ 0.88392779,  0.59139899,  0.41895111]])

Parallele Berechnung mit workers

>>> from concurrent.futures import ThreadPoolExecutor
>>> with ThreadPoolExecutor() as executor:
...     res = cubature(
...         f=f,
...         a=np.array([0, 0, 0, 0, 0, 0, 0]),
...         b=np.array([1, 1, 1, 1, 1, 1, 1]),
...         rtol=1e-5,
...         rule="genz-malik",
...         args=(r, alphas),
...         workers=executor.map,
...      )
>>> res.estimate
 array([[-0.79812452,  0.35246913, -0.52273628],
        [ 0.88392779,  0.59139899,  0.41895111]])

2D-Integral mit unendlichen Grenzen:

\[\int^{ \infty }_{ -\infty } \int^{ \infty }_{ -\infty } e^{-x^2-y^2} \text dy \text dx\]

>>> def gaussian(x):
...     return np.exp(-np.sum(x**2, axis=-1))
>>> res = cubature(gaussian, [-np.inf, -np.inf], [np.inf, np.inf])
>>> res.estimate
 3.1415926

1D-Integral mit Singularitäten, die mit points vermieden werden

\[\int^{ 1 }_{ -1 } \frac{\sin(x)}{x} \text dx\]

Es ist notwendig, den Parameter points zu verwenden, um die Auswertung von f am Ursprung zu vermeiden.

>>> def sinc(x):
...     return np.sin(x)/x
>>> res = cubature(sinc, [-1], [1], points=[[0]])
>>> res.estimate
 1.8921661