11

Wenn man numerische Ableitungen berechnen möchte, ist die von Bengt Fornberg hier vorgestellte (und hier beschriebene ) Methode sehr praktisch (sowohl präzise als auch einfach zu implementieren). Da das ursprüngliche Papier aus dem Jahr 1988 stammt, möchte ich wissen, ob es heute eine bessere Alternative gibt (als (oder fast als) einfach und präziser).

finite-difference precision Vincent
quelle

1

Schwer zu sagen, ohne zu wissen, was Sie unterscheiden möchten. Haben Sie über eine automatische Unterscheidung nachgedacht ?

Biswajit Banerjee

@BiswajitBanerjee: Für endliche Differenzkoeffizienten gilt die automatische Differenzierung nicht.

Geoff Oxberry

@GeoffOxberry: Ich bezog mich auf das eigentliche physikalische Problem, dh den wissenschaftlichen Teil der "Computerwissenschaft".

Biswajit Banerjee

@ Vincent: Versuchen Sie, eine Funktion oder eine Tabelle zu unterscheiden? Wenn Tabellendaten, sind die Tabellendaten verrauscht? Versuchen Sie, eine PDE zu diskretisieren?

user14717

9

Überblick

Gute Frage. Es gibt einen Artikel mit dem Titel "Verbesserung der Genauigkeit der Matrixdifferenzierungsmethode für beliebige Kollokationspunkte" von R. Baltensperger. Meiner Meinung nach ist das keine große Sache, aber es hat einen Punkt (der bereits vor dem Erscheinen im Jahr 2000 bekannt war): Es betont die Bedeutung einer genauen Darstellung der Tatsache, dass die Ableitung der konstanten Funktion $f(x)=1$ sollte Null sein (dies gilt genau im mathematischen Sinne, aber nicht unbedingt in der numerischen Darstellung).

Es ist leicht zu erkennen, dass dies erfordert, dass die Zeilensummen der n-ten Ableitungsmatrizen $D^{(n)}$ Null sind. Es ist üblich, diese Einschränkung durch Anpassen des diagonalen Eintrags zu erzwingen, dh durch Setzen von

\begin{matrix} (1) & D_{j j}^{(n)} := - \sum_{\begin{matrix} i = 1 \\ i \neq j \end{matrix}}^{N} D_{i j} . \end{matrix}

$\tag 1 D^{(n)}_{jj} := -\sum_{\substack{i=1\\i \neq j }}^N D_{ij}\,.$ Es ist klar, dass diese Funktion aufgrund von Rundungsfehlern bei Gleitkommaberechnungen nicht genau für die Arbeit an einem Computer gilt. Überraschender ist, dass diese Fehler bei Verwendung der analytischen Formeln für die Ableitungsmatrix (die für viele klassische Kollokationspunkte verfügbar sind, z. B. Gauss-Lobatto) noch schwerwiegender sind.

In dem Papier (und den darin enthaltenen Referenzen) heißt es nun, dass der Fehler der Ableitung in der Größenordnung der Abweichung der Zeilensummen von Null liegt. Ziel ist es daher, diese numerisch so klein wie möglich zu halten.

Numerische Tests

Der gute Punkt ist, dass das Fornberg-Verfahren in dieser Hinsicht recht gut zu sein scheint. Im Bild unten habe ich das Verhalten der exakten, dh analytischen Matrix der ersten Ableitung und der vom Fornberg-Algorithmus abgeleiteten Matrix für eine unterschiedliche Anzahl von Chebyshev-Lobatto-Kollokationspunkten verglichen.

In Anbetracht der Aussage in dem zitierten Artikel impliziert dies wiederum, dass der Fornberg-Algorithmus genauere Ergebnisse für die Ableitung liefert.

Um dies zu beweisen, verwende ich die gleiche Funktion wie in der Arbeit,

\begin{matrix} (2) & f (x) = \frac{1}{1 + x^{2}} . \end{matrix}

$\tag 2 f(x) = \frac{1}{1+x^2}\,.$

\begin{matrix} (3) & E_{n} = max_{i \in {0, \dots, n}} | f^{'} (x_{i}) - \sum_{j = 1}^{n} D_{i j} f (x_{j}) | . \end{matrix}

$\tag 3 E_n = \max_{i\in \{0,\ldots,n\}} \bigg| f^\prime(x_i) - \sum_{j=1}^n D_{ij} f(x_j) \bigg|\,.$

\begin{matrix} (4) & {\tilde{D}}_{j j} = D_{j j} - (\sum_{i = 1}^{n} D_{j i}), for all j . \end{matrix}

$\tag 4 \tilde D_{jj} = D_{jj} - \left(\sum_{i=1}^n D_{ji} \right), \qquad \text{for all}\ j\,.$

Fazit

Zusammenfassend scheint die Fornberg-Methode ziemlich genau zu sein, im Fall sogar um etwa 3 Größenordnungen genauer als die Ergebnisse der analytischen Formeln. Dies sollte für die meisten Anwendungen genau genug sein. Darüber hinaus ist dies bemerkenswert, da Fornberg diese Tatsache nicht explizit in seine Methode einzubeziehen scheint (zumindest wird dies in den beiden Fornberg-Papieren nicht erwähnt). $N=512$

Eine weitere Größenordnung kann für dieses Beispiel durch die einfache Einbeziehung von Gleichung (4) erreicht werden. Da dies ein recht einfacher Ansatz ist und nur einmal für jedes Derivat angewendet wird, sehe ich keinen Grund, ihn nicht zu verwenden.

Die Methode aus dem Baltensperger-Papier, die einen differenzierteren Ansatz zur Bewertung der Summe in Gleichung (1) verwendet, um Rundungsfehler zu reduzieren, ergibt ungefähr die gleiche Größenordnung für den Fehler. Zumindest für dieses Beispiel entspricht es in etwa der oben beschriebenen "Adjusted Fornberg" -Methode.

Davidhigh
quelle

4

Angenommen, Sie versuchen, eine numerische Implementierung einer stetigen Funktion zu unterscheiden, gibt es eine Vielzahl von Methoden:

1) Automatische Differenzierung. Die genaueste und allgemeinste Methode. Das Codieren ist schmerzhaft und erfordert eine Überladung des Bedieners und eine argumentabhängige Suche. Belastet den Benutzer, diese Konzepte zu verstehen. Kämpft auch mit entfernbaren Singularitäten, wie z. B. der Differenzierung von sinc bei . $x= 0$

2) Eine Chebyshev-Transformation. Projizieren Sie Ihre Funktion auf eine Reihe von Chebyshev-Polynomen und unterscheiden Sie die drei Termwiederholungen. Super schnell, sehr genau. Voraussetzung ist jedoch, dass Sie über eine kompakt unterstützte Domäne von Interesse verfügen. Außerhalb der ausgewählten Domäne ist die Wiederholung mit drei Begriffen instabil. $[a,b]$

3) Endliche Differenzierung. In 1D unterschätzt; Siehe Nick Highams Tipps und Tricks zum numerischen Rechnen . Die Idee ist, dass Sie keine Schrittweite auswählen müssen, wenn Sie den Kürzungsfehler und den Rundungsfehler ausgleichen. es kann automatisch ausgewählt werden. In Boost wird diese Idee verwendet, um (standardmäßig) 6/7 der korrekten Ziffern für den Typ wiederherzustellen. (Higham zeigt nur die Idee für den einfacheren Fall von 1/2 der richtigen Ziffern, aber die Idee kann leicht erweitert werden.) Die Koeffizienten stammen aus Fornbergs Tabelle mit gleichem Abstand, aber die Schrittgröße wird unter der Annahme gewählt, dass die Funktion auf 1ULP ausgewertet werden kann Richtigkeit. Der Nachteil ist, dass 2 Funktionsauswertungen erforderlich sind, um die Hälfte der Ziffern des Typs wiederherzustellen, 4, um 3/4 der Ziffern wiederherzustellen, usw. In 1D kein schlechtes Geschäft. In höheren Dimensionen ist es katastrophal.

4) Die komplexe Schrittableitung. Verwenden Sie . Nehmen Sie als Einheitsrundung und dies wird fast jedes Bit korrekt wiederherstellen. Es ist jedoch ein bisschen Betrug, weil es im Allgemeinen schwieriger ist, eine Funktion in der komplexen Ebene zu implementieren, als ihre reale Ableitung von Hand zu codieren. Immer noch eine coole Idee und unter bestimmten Umständen nützlich. $f'(x) \approx \Im(f(x+ih))$ $h$

user14717
quelle

1

Mir ist nicht bekannt, dass jemand Fornbergs Algorithmus verbessert hat (siehe auch sein etwas neueres Papier ). Abgesehen davon scheint es mir nicht richtig zu sein, seinen Algorithmus als Methode zur Berechnung numerischer Ableitungen zu betrachten. Alles, was er getan hat, ist, einen effizienten Algorithmus zur Berechnung von Gewichten für Finite-Differenzen-Methoden abzuleiten . Der Vorteil seiner Methode ist, dass Sie die Gewichte für alle Derivate bis zur gewünschten Ableitung auf einmal erhalten.

Brian Zatapatique
quelle

Ich würde der Aussage nicht zustimmen, "es ist nicht richtig, seinen Algorithmus als einen Weg zur Berechnung numerischer Ableitungen zu betrachten". Sobald Sie die Gewichte für den Punkt ausgewertet haben , können Sie die numerische Ableitung jeder Funktion direkt über .

w_{i}

$w_i$

z

$z$

f (x)

$f(x)$

f^{'} (z) = \sum_{i} w_{i} f (x_{i})

$f^\prime(z) = \sum_i w_i f(x_i)$

Davidhigh

@davidhigh: Wenn Sie Fornbergs Artikel lesen, sprechen sie über die Berechnung von Gewichten für Näherungen mit endlichen Differenzen, nicht über die Berechnung der Näherungen selbst. Natürlich können Sie den Algorithmus auch zum Berechnen der Ableitungen verwenden, aber es ist sinnvoller, die Gewichte zu berechnen, zu speichern und sie dann wiederholt zum Berechnen der ungefähren Ableitungen zu verwenden. Andernfalls berechnen Sie die Gewichte wiederholt, wenn keine Notwendigkeit besteht.

Brian Zatapatique

1

Ein einfacheres Schema

Zusätzlich zu meiner anderen Antwort, die sich mehr mit einer Erweiterung der Fornberg-Methode befasst, werde ich hier die Frage nach einfacheren Alternativen ansprechen.

Dazu skizziere ich ein alternatives Schema, das die Ableitungskoeffizienten der Lagrange-Interpolation direkter erzeugt. Die Implementierung erfordert nur wenige Codezeilen, funktioniert für beliebige Gitter und ist nach meinen ersten Experimenten genauso genau wie die von Fornberg.

Grundlage der Implementierung ist die imaginäre Schrittableitung wobei ist eine Variable in der Reihenfolge der Maschinengenauigkeit. Es ist bekannt, dass die Ableitung im imaginären Schritt Ableitungswerte stabil erzeugt und nicht unter der numerischen Instabilität einer Finite-Differenzen-Implementierung mit leidet .

f^{'} (x) = \frac{1}{h} Im (f (x + i h)),

$f'(x) \ = \ \frac{1}{h}\,\text{Im}\big(f(x+ih)\big)\,,$

h

$h$

h \to 0

$h\rightarrow 0$

Die zweite Zutat ist das Lagrange-Interpolationspolynom auf dem Gitter das mit einer der baryzentrischen Formen ausgewertet wird, z. B. denen Um die Ableitung in komplexen Schritten verwenden zu können, muss sichergestellt werden, dass diese Formeln auch für komplexe Formeln funktionieren Argumente . Darüber hinaus bezeichnen wir für eine gegebene Funktion f (x) und einen Koeffizientenvektor das Interpolationspolynom durch mit $L_i(x)$ $\{x_1,\ldots, x_N\}$

L_{i} (z) = {\begin{cases} 1 & if z = x_{i} \\ \frac{\frac{μ_{i}}{z - x_{k}}}{\sum_{k} \frac{μ_{k}}{z - x_{k}}} & otherwise . \end{cases}

$L_i(z) \ = \ \begin{cases}1 & \text{if }z = x_i\\\ \frac{\frac {\mu_i}{z-x_k}}{ \sum_k\frac {\mu_k}{z-x_k} } & \text{otherwise}\,. \end{cases}$

μ_{i} = \frac{1}{\prod_{k \neq i} (x_{i} - x_{k})}

$\mu_i = \frac{1}{\prod_{k\neq i} (x_i - x_k)}$

z

$z$

f_{i} = f (x_{i})

$f_i = f(x_i)$

(x_{1}, f_{i})

$(x_1,f_i)$

P (x; f) = \sum_{i = 1}^{N} f_{i} L_{i} (x) .

$P(x;f) = \sum_{i=1}^N f_i L_i(x)\,.$

Algorithmus

Der Algorithmus wird im Folgenden skizziert. Es hat die gleichen Eingabe- und Ausgabeparameter wie Fornberg, ist jedoch viel einfacher.

Eingang:

$x$ : ein Gitter mit N verschiedenen Gitterpunkten
$ord$ : abgeleitete Reihenfolge
z: ein Punkt, an dem die Ableitung bewertet werden soll
Möglicherweise: eine Funktion oder ihre Funktionswerte an den Gitterpunkten (nur für Ausgangsvariante 2 erforderlich) $f(x)$ $f_i$

Initialisierung

Initialisieren Sie die Lagrange-Polynome durch baryzentrische Interpolation. $L$
Initialisieren Sie ein Array von Matrizen für $N\times N$ $D^{(o)}$ $o=0, \ldots, ord$
Setze , dh die Einheitsmatrix der Dimension $D^{(0)} := I_N$ $N\times N$
Setze $o := 0$

Algorithmus

Während : $o < ord$

Berechnen Sie mit der komplexen für alles und . Hier bezeichnet die te Zeile von . $D^{(o+1)}_{ik} = CSD \bigg( P\big(x_k; D^{(o)}_{i}\big) \bigg)$ $CSD$ $i$ $k$
$D^{(o)}_{i}$ $i$ $D^{(o)}$
Setze o = o + 1;

Entscheiden Sie, was ausgegeben werden soll :

Der Vektor der endlichen Differenzkoeffizienten am Punkt , wobei . Das macht Fornberg. $d^{(ord)}$ $z$ $d_i = P(z;D^{(ord)}_i)$
Die Interpolationsfunktion zur Ableitung der Bestellung . Dazu müssen Sie eine Funktion bzw. die Funktionswerte bei zum Algorithmus. $p^{(ord)}(x) = \sum_{i=1}^N f(x_i) P(x;D^{(ord)}_i)$ $f^{(ord)}(x)$ $ord$ $f_i$ $x_i$
Eine Metafunktion die die Interpolationsfunktion von Variante 2 zurückgibt, jedoch für eine beliebige Funktion , die an den Gitterpunkten interpoliert werden soll. $\text{diff}(\text{int } ord, \text{function } g)$ $g$

Persönlich mag ich Variante 3. am meisten.

Analyse des Algorithmus

Wie bei Fornberg ist dieser Algorithmus . Ich werde empirischere Ergebnisse in Bezug auf Genauigkeit, Stabilität usw. veröffentlichen, wenn ich die Zeit finde. $\mathcal O(ord \cdot N^2)$

Davidhigh
quelle

0

Gehen Sie wie folgt vor, um die Genauigkeit der numerischen Differenzierung zu erhöhen:

1) Wählen Sie Ihre bevorzugte hochpräzise "Standard" -Methode basierend auf einer Schrittgröße h .

2) Berechnen Sie den Wert der Ableitung mit der in 1) gewählten Methode viele Male mit unterschiedlichen, aber vernünftigen Schrittgrößen h . Jedes Mal können Sie h als Zufallszahl aus dem Intervall (0,5 * H / 10, 1,5 * H / 10) auswählen, wobei H eine geeignete Schrittgröße für die von Ihnen verwendete Methode ist.

3) Durchschnitt der Ergebnisse.

Ihr Ergebnis kann 2-3 Größenordnungen im absoluten Fehler wrt gewinnen. das nicht gemittelte Ergebnis.

https://arxiv.org/abs/1706.10219

F. Jatpil
quelle

3

Willkommen bei SciComp.SE! Diese Antwort wäre noch besser, wenn Sie die Methode kurz zusammenfassen würden.

Christian Clason

2

Auch Ihr Benutzername deutet darauf hin, dass Sie dieses Papier verfasst haben. Bitte lesen Sie unsere Richtlinien zur Eigenwerbung und bearbeiten Sie Ihren Beitrag entsprechend.

Wrzlprmft

Ich denke ehrlich, dass die negative Abstimmung unfair ist, wenn meine Antwort tatsächlich auf eine gültige Antwort hinweist.

F. Jatpil

1

Ich auch ... also nimm meine +1 ... :-)