scipy.integrate.

tanhsinh#

scipy.integrate.tanhsinh(f, a, b, *, args=(), log=False, maxlevel=None, minlevel=2, atol=None, rtol=None, preserve_shape=False, callback=None)[Quelle]#

Evaluiert ein konvergentes Integral numerisch mittels Tanh-Sinh-Quadratur.

In der Praxis erreicht die Tanh-Sinh-Quadratur für viele Integranden eine quadratische Konvergenz: die Anzahl der genauen *Ziffern* skaliert annähernd linear mit der Anzahl der Funktionsauswertungen [1].

Entweder oder beide Integrationsgrenzen können unendlich sein, und Singularitäten an den Endpunkten sind akzeptabel. Divergente Integrale und Integranden mit nicht-endlichen Ableitungen oder Singularitäten innerhalb eines Intervalls sind nicht abgedeckt, aber letztere können durch Aufruf von tanhsinh für jedes Teilintervall einzeln ausgewertet werden.

Parameter:
faufrufbar

Die zu integrierende Funktion. Die Signatur muss lauten

f(xi: ndarray, *argsi) -> ndarray

wobei jedes Element von xi eine endliche reelle Zahl ist und argsi ein Tupel ist, das eine beliebige Anzahl von Arrays enthalten kann, die mit xi broadcastbar sind. f muss eine elementweise Funktion sein: siehe Dokumentation des Parameters preserve_shape für Details. Sie darf das Array xi oder die Arrays in argsi nicht verändern. Wenn f einen Wert mit komplexem dtype zurückgibt, wenn er an einem der Endpunkte ausgewertet wird, haben nachfolgende Argumente x einen komplexen dtype (aber einen imaginären Teil von Null).

a, bfloat array_like

Reelle untere und obere Integrationsgrenzen. Müssen miteinander und mit Arrays in args broadcastbar sein. Elemente können unendlich sein.

argstuple of array_like, optional

Zusätzliche positionelle Array-Argumente, die an f übergeben werden sollen. Arrays müssen miteinander und mit den Arrays von a und b broadcastbar sein. Wenn die aufzurufende Funktion, für die die Wurzel gesucht wird, Argumente benötigt, die nicht mit x broadcastbar sind, wickeln Sie diese aufzurufende Funktion mit f so ein, dass f nur x und broadcastbare *args akzeptiert.

logbool, default: False

Wenn True gesetzt, bedeutet dies, dass f den Logarithmus des Integranden zurückgibt und dass atol und rtol als Logarithmen der absoluten und relativen Fehler ausgedrückt werden. In diesem Fall enthält das Ergebnisobjekt den Logarithmus des Integrals und des Fehlers. Dies ist nützlich für Integranden, bei denen numerisches Underflow oder Overflow zu Ungenauigkeiten führen würden. Wenn log=True, muss der Integrand (die Exponentialfunktion von f) reell sein, er kann jedoch negativ sein, in welchem Fall der Logarithmus des Integranden eine komplexe Zahl mit einem imaginären Teil ist, der ein ungerades Vielfaches von π ist.

maxlevelint, default: 10

Die maximale Verfeinerungsstufe des Algorithmus.

Auf der nullten Stufe wird f einmal aufgerufen, was 16 Funktionsauswertungen durchführt. Auf jeder nachfolgenden Stufe wird f einmal mehr aufgerufen, wodurch sich die Anzahl der durchgeführten Funktionsauswertungen ungefähr verdoppelt. Daher wird für viele Integranden jede aufeinanderfolgende Stufe die Anzahl der genauen Ziffern im Ergebnis verdoppeln (bis zu den Grenzen der Gleitkommapräzision).

Der Algorithmus wird nach Abschluss der Stufe maxlevel oder nach Erfüllung einer anderen Abbruchbedingung beendet, je nachdem, was zuerst eintritt.

minlevelint, default: 2

Die Stufe, auf der die Iteration beginnt (Standard: 2). Dies ändert nicht die Gesamtzahl der Funktionsauswertungen oder die Abszissen, an denen die Funktion ausgewertet wird. Es ändert nur die *Anzahl der Aufrufe* von f. Wenn minlevel=k ist, wird der Integrand an allen Abszissen von den Stufen 0 bis k in einem einzigen Aufruf ausgewertet. Beachten Sie, dass, wenn minlevel größer als maxlevel ist, das bereitgestellte minlevel ignoriert wird und minlevel gleich maxlevel gesetzt wird.

atol, rtolfloat, optional

Absolute Abbruchtoleranz (Standard: 0) und relative Abbruchtoleranz (Standard: eps**0.75, wobei eps die Präzision des Ergebnis-Datentyps ist), jeweils. Muss nicht-negativ und endlich sein, wenn log False ist, und muss als Logarithmus einer nicht-negativen und endlichen Zahl ausgedrückt werden, wenn log True ist. Die Iteration wird gestoppt, wenn res.error < atol oder res.error < res.integral * rtol.

preserve_shapebool, default: False

Im Folgenden bezieht sich „Argumente von f“ auf das Array xi und alle Arrays innerhalb von argsi. Sei shape die broadcastete Form von a, b und allen Elementen von args (was konzeptionell von xi und argsi, die in f übergeben werden, getrennt ist).

  • Wenn preserve_shape=False (Standard), muss f Argumente beliebiger broadcastbarer Formen akzeptieren.

  • Wenn preserve_shape=True, muss f Argumente der Form shape *oder* shape + (n,) akzeptieren, wobei (n,) die Anzahl der Abszissen ist, an denen die Funktion ausgewertet wird.

In beiden Fällen muss das von f zurückgegebene Array für jedes Skalar-Element xi[j] innerhalb von xi das Skalar f(xi[j]) an derselben Stelle enthalten. Folglich ist die Form des Ausgabewerts immer die Form des Eingabewerts xi.

Siehe Beispiele.

callbackcallable, optional

Eine optionale, vom Benutzer bereitgestellte Funktion, die vor der ersten Iteration und nach jeder Iteration aufgerufen wird. Aufruf als callback(res), wobei res ein _RichResult ist, ähnlich dem von tanhsinh zurückgegebenen (aber die Werte aller Variablen des aktuellen Iterationsschritts enthaltend). Wenn callback eine StopIteration auslöst, wird der Algorithmus sofort beendet und tanhsinh gibt ein Ergebnisobjekt zurück. callback darf res oder dessen Attribute nicht verändern.

Rückgabe:
res_RichResult

Ein Objekt, das einem Instanz von scipy.optimize.OptimizeResult ähnelt, mit den folgenden Attributen. (Die Beschreibungen sind so formuliert, als ob die Werte Skalare wären; wenn f jedoch ein Array zurückgibt, sind die Ausgaben Arrays gleicher Form.)

successBool-Array

True, wenn der Algorithmus erfolgreich beendet wurde (Status 0). False andernfalls.

statusInteger-Array

Eine Ganzzahl, die den Exit-Status des Algorithmus darstellt.

0 : Der Algorithmus konvergierte zu den angegebenen Toleranzen. -1 : (ungenutzt) -2 : Die maximale Anzahl von Iterationen wurde erreicht. -3 : Ein nicht-endlicher Wert wurde angetroffen. -4 : Die Iteration wurde durch callback beendet. 1 : Der Algorithmus läuft normal (nur in callback).

integralfloat array

Eine Schätzung des Integrals.

errorFloat-Array

Eine Schätzung des Fehlers. Nur verfügbar, wenn Stufe zwei oder höher abgeschlossen ist; andernfalls NaN.

maxlevelint array

Die maximale verwendete Verfeinerungsstufe.

nfevInteger-Array

Die Anzahl der Punkte, an denen f ausgewertet wurde.

Siehe auch

quad

Hinweise

Implementiert den Algorithmus, wie in [1] beschrieben, mit geringfügigen Anpassungen für die Gleitkommaarithmetik, einschließlich einiger, die von [2] und [3] beschrieben wurden. Das Tanh-Sinh-Schema wurde ursprünglich in [4] eingeführt.

Zwei Fehlerschätzungsmechanismen werden in [1] Abschnitt 5 beschrieben: einer versucht, eine quadratische Konvergenz zu erkennen und auszunutzen; der andere vergleicht einfach die Integralschätzungen aufeinanderfolgender Stufen. Obwohl keiner theoretisch rigoros oder konservativ ist, funktionieren beide in der Praxis gut. Unsere Fehlerschätzung verwendet das Minimum dieser beiden Mechanismen mit einer unteren Grenze von eps * res.integral.

Aufgrund von Gleitkommafehlern in den Abszissen kann die Funktion während der Iterationen an den Endpunkten des Intervalls ausgewertet werden, aber die von der Funktion an den Endpunkten zurückgegebenen Werte werden ignoriert.

Referenzen

[1] (1,2,3)

Bailey, David H., Karthik Jeyabalan und Xiaoye S. Li. „A comparison of three high-precision quadrature schemes.“ Experimental Mathematics 14.3 (2005): 317-329.

[2]

Vanherck, Joren, Bart Sorée und Wim Magnus. „Tanh-sinh quadrature for single and multiple integration using floating-point arithmetic.“ arXiv preprint arXiv:2007.15057 (2020).

[3]

van Engelen, Robert A. „Improving the Double Exponential Quadrature Tanh-Sinh, Sinh-Sinh and Exp-Sinh Formulas.“ https://www.genivia.com/files/qthsh.pdf

[4]

Takahasi, Hidetosi und Masatake Mori. „Double exponential formulas for numerical integration.“ Publications of the Research Institute for Mathematical Sciences 9.3 (1974): 721-741.

Beispiele

Evaluiert das Gaußsche Integral

>>> import numpy as np
>>> from scipy.integrate import tanhsinh
>>> def f(x):
...     return np.exp(-x**2)
>>> res = tanhsinh(f, -np.inf, np.inf)
>>> res.integral  # true value is np.sqrt(np.pi), 1.7724538509055159
1.7724538509055159
>>> res.error  # actual error is 0
4.0007963937534104e-16

Der Wert der Gaußschen Funktion (Glockenkurve) ist für Argumente, die ausreichend weit von Null entfernt sind, annähernd Null, sodass der Wert des Integrals über ein endliches Intervall annähernd derselbe ist.

>>> tanhsinh(f, -20, 20).integral
1.772453850905518

Bei ungünstigen Integrationsgrenzen kann das Integrationsschema jedoch die wichtige Region nicht finden.

>>> tanhsinh(f, -np.inf, 1000).integral
4.500490856616431

In solchen Fällen oder wenn es Singularitäten innerhalb des Intervalls gibt, brechen Sie das Integral in Teile auf, wobei die Endpunkte an den wichtigen Punkten liegen.

>>> tanhsinh(f, -np.inf, 0).integral + tanhsinh(f, 0, 1000).integral
1.772453850905404

Für Integrationen mit sehr großen oder sehr kleinen Beträgen verwenden Sie die Log-Integration. (Zu Illustrationszwecken zeigt das folgende Beispiel einen Fall, in dem sowohl reguläre als auch Log-Integration funktionieren, aber bei extremeren Integrationsgrenzen würde die Log-Integration das Underflow vermeiden, das bei der normalen Auswertung des Integrals auftritt.)

>>> res = tanhsinh(f, 20, 30, rtol=1e-10)
>>> res.integral, res.error
(4.7819613911309014e-176, 4.670364401645202e-187)
>>> def log_f(x):
...     return -x**2
>>> res = tanhsinh(log_f, 20, 30, log=True, rtol=np.log(1e-10))
>>> np.exp(res.integral), np.exp(res.error)
(4.7819613911306924e-176, 4.670364401645093e-187)

Die Integrationsgrenzen und Elemente von args können broadcastbare Arrays sein, und die Integration wird elementweise durchgeführt.

>>> from scipy import stats
>>> dist = stats.gausshyper(13.8, 3.12, 2.51, 5.18)
>>> a, b = dist.support()
>>> x = np.linspace(a, b, 100)
>>> res = tanhsinh(dist.pdf, a, x)
>>> ref = dist.cdf(x)
>>> np.allclose(res.integral, ref)
True

Standardmäßig ist preserve_shape False, und daher kann die aufrufbare Funktion f mit Arrays beliebiger broadcastbarer Formen aufgerufen werden. Zum Beispiel

>>> shapes = []
>>> def f(x, c):
...    shape = np.broadcast_shapes(x.shape, c.shape)
...    shapes.append(shape)
...    return np.sin(c*x)
>>>
>>> c = [1, 10, 30, 100]
>>> res = tanhsinh(f, 0, 1, args=(c,), minlevel=1)
>>> shapes
[(4,), (4, 34), (4, 32), (3, 64), (2, 128), (1, 256)]

Um zu verstehen, woher diese Formen stammen – und um besser zu verstehen, wie tanhsinh genaue Ergebnisse liefert – beachten Sie, dass höhere Werte von c mit höherfrequenten Sinusoiden korrespondieren. Die höherfrequenten Sinusoiden machen den Integranden komplizierter, sodass mehr Funktionsauswertungen erforderlich sind, um die Zielgenauigkeit zu erreichen.

>>> res.nfev
array([ 67, 131, 259, 515], dtype=int32)

Die anfängliche shape, (4,), entspricht der Auswertung des Integranden an einer einzelnen Abszisse und allen vier Frequenzen; dies wird zur Eingabevalidierung und zur Bestimmung der Größe und des Datentyps der Arrays, die Ergebnisse speichern, verwendet. Die nächste Form entspricht der Auswertung des Integranden an einem anfänglichen Gitter von Abszissen und allen vier Frequenzen. Nachfolgende Aufrufe der Funktion verdoppeln die Gesamtzahl der Abszissen, an denen die Funktion ausgewertet wurde. In späteren Funktionsauswertungen wird der Integrand jedoch an weniger Frequenzen ausgewertet, da das entsprechende Integral bereits die erforderliche Toleranz erreicht hat. Dies spart Funktionsauswertungen zur Leistungsverbesserung, erfordert jedoch, dass die Funktion Argumente beliebiger Formen akzeptiert.

„Vektorwertige“ Integranden, wie sie für die Verwendung mit scipy.integrate.quad_vec geschrieben werden, erfüllen diese Anforderung wahrscheinlich nicht. Betrachten Sie zum Beispiel

>>> def f(x):
...    return [x, np.sin(10*x), np.cos(30*x), x*np.sin(100*x)**2]

Dieser Integrand ist nicht mit tanhsinh in der vorliegenden Form kompatibel; beispielsweise ist die Form der Ausgabe nicht dieselbe wie die Form von x. Eine solche Funktion *könnte* mit der Einführung zusätzlicher Parameter in eine kompatible Form umgewandelt werden, aber dies wäre umständlich. In solchen Fällen wäre eine einfachere Lösung die Verwendung von preserve_shape.

>>> shapes = []
>>> def f(x):
...     shapes.append(x.shape)
...     x0, x1, x2, x3 = x
...     return [x0, np.sin(10*x1), np.cos(30*x2), x3*np.sin(100*x3)]
>>>
>>> a = np.zeros(4)
>>> res = tanhsinh(f, a, 1, preserve_shape=True)
>>> shapes
[(4,), (4, 66), (4, 64), (4, 128), (4, 256)]

Hier ist die broadcastete Form von a und b (4,). Mit preserve_shape=True kann die Funktion mit dem Argument x der Form (4,) oder (4, n) aufgerufen werden, und dies beobachten wir auch.

Auf dieser Seite