2009-10-13 33 views
15

Wie wird die Ableitung eines f(x) normalerweise programmgesteuert berechnet, um maximale Genauigkeit zu gewährleisten?Implementieren der Ableitung in C/C++

Ich implementiere die Newton-Raphson Methode, und es erfordert die Ableitung der Ableitung einer Funktion.

+0

Zeigen Sie uns, was Sie bisher getan haben. – Lazarus

+4

willst du, dass ich gefeuert werde :)? – vehomzzz

+5

Die Genauigkeit der Newton-Methode hängt (nur) nicht von der Genauigkeit der Ableitung ab, eine grobe Approximation wird auch funktionieren (im Extremfall erhalten Sie die secant-Methode), aber es wird ein wenig langsamer sein, um die gleiche Genauigkeit zu erhalten (Mehr Iterationen benötigt). – fortran

Antwort

51

Ich stimme mit erikkallen überein, dass (f (x + h) - f (x - h))/2h der übliche Ansatz zur numerischen Approximation von Derivaten ist. Allerdings ist die richtige Schrittgröße h etwas subtil.

Der Näherungsfehler in (f (x + h) - f (x - h))/2h nimmt ab, wenn h kleiner wird, was bedeutet, dass Sie h so klein wie möglich nehmen sollten. Aber wenn h kleiner wird, erhöht sich der Fehler von der Gleitkomma-Subtraktion, da der Zähler fast gleiche Zahlen subtrahieren muss. Wenn h zu klein ist, können Sie bei der Subtraktion viel Genauigkeit verlieren. In der Praxis müssen Sie also einen nicht zu kleinen Wert von h wählen, der die Kombination Approximation Fehler und numerische Fehler minimiert.

Als Faustregel kann man h = SQRT (DBL_EPSILON) verwenden, wobei DBL_EPSILON die kleinste doppelte Genauigkeitszahl e ist, so dass 1 + e! = 1 in Maschinengenauigkeit ist. DBL_EPSILON ist etwa 10^-15, also könnten Sie h = 10^-7 oder 10^-8 verwenden.

Weitere Einzelheiten finden Sie unter notes zum Auswählen der Schrittgröße für Differentialgleichungen.

+0

+1 Sehr interessant, in es. – vehomzzz

+4

+1 für Erwähnung und Erklärung des Stepsize-Dilemmas – sellibitze

+3

Ich denke, Ihre Faustregel geht davon aus, dass Sie eine Regel erster Ordnung verwenden, um die Ableitung zu approximieren. Die zentrale Differenzregel, die Sie erwähnen, ist jedoch die zweite Ordnung, und die entsprechende Faustregel lautet h = EPSILON^(1/3), was ungefähr 10^(- 5) ist, wenn doppelte Genauigkeit verwendet wird. –

5
fprime(x) = (f(x+dx) - f(x-dx))/(2*dx) 

für einige kleine dx.

+2

Wie klein? danke – vehomzzz

+1

Numerische Rezepte hat einige Kommentare dazu http://books.google.co.uk/books?id=1aAOdzK3FegC&printsec=frontcover&dq=related:ISBN0521437202#v=onepage&q=newton-Raphson&f=false – Mark

10

Newton_Raphson nimmt an, dass Sie zwei Funktionen f (x) und seine Ableitung f '(x) haben können. Wenn Sie die Ableitung nicht als Funktion verfügbar haben und die Ableitung von der ursprünglichen Funktion schätzen müssen, sollten Sie einen anderen Wurzelfindungsalgorithmus verwenden.

Wikipedia root finding gibt mehrere Vorschläge wie jede numerische Analyse Text.

5

Was wissen Sie über f (x)? Wenn Sie nur f als Blackbox haben, können Sie nur die Ableitung numerisch approximieren. Aber die Genauigkeit ist normalerweise nicht so gut.

Sie können viel besser tun, wenn Sie den Code berühren können, der f berechnet. Versuchen Sie "automatic differentiation". Da sind einige nette Bibliotheken dafür verfügbar. Mit ein wenig Bibliothek Magie können Sie Ihre Funktion leicht in etwas konvertieren, das die Ableitung automatisch berechnet. Ein einfaches C++ - Beispiel finden Sie in der deutschen Diskussion source code in this.

2

Zusätzlich zu John D. Cooks Antwort oben ist es wichtig, nicht nur die Fließkomma-Genauigkeit, sondern auch die Robustheit der Funktion f (x) zu berücksichtigen. Zum Beispiel ist es im Finanzbereich ein häufiger Fall, dass f (x) tatsächlich eine Monte-Carlo-Simulation ist und der Wert von f (x) ein gewisses Rauschen aufweist. Die Verwendung einer sehr kleinen Schrittgröße kann in diesen Fällen die Genauigkeit der Ableitung stark beeinträchtigen.

8

alt text

alt text

1) Erster Fall:

alt text

alt text - relative Rundungsfehler, etwa 2^{- 16} für Doppel- und 2^{- 7} für den Schwimmer.

Wir können Gesamtfehler berechnen:

alt text

Angenommen, Sie doppelte floating Betrieb verwenden. Somit ist der optimale Wert von h 2sqrt (DBL_EPSILON/f '' (x)). Sie wissen nicht f '' (x). Aber Sie müssen diesen Wert schätzen. Wenn zum Beispiel f '' (x) ungefähr 1 ist, dann ist der optimale Wert h 2^{- 7}, aber wenn f '' (x) ist ungefähr 10^6 dann der optimale Wert von h ist 2^{- 10}!

2) Zweiter Fall:

alt text

Hinweis, dass die zweite Näherungsfehler zu 0 tendiert schneller als erstes. Aber wenn f ‚‘ '(x) ist sehr lagre dann die erste Option ist mehr bevorzugt:

alt text

Beachten Sie, dass im ersten Fall h e proportional ist, aber im zweiten Fall h ist proportional zu e^{1/3}. Für Double Floating-Operationen ist e^{1/3} 2^{- 5} oder 2^{- 6}. (Ich nehme an, dass f '' '(x) etwa 1 ist).


Welcher Weg ist besser? Es ist unbekannt, wenn Sie nicht wissen, f '' (x) und f '' '(x) oder Sie können diese Werte nicht schätzen. Es wird angenommen, dass die zweite Option vorzuziehen ist. Aber wenn Sie wissen, dass f '' '(x) sehr groß ist, verwenden Sie zuerst.

Was ist der optimale Wert von h? Angenommen, dass f '' (x) und f '' '(x) ungefähr 1 sind. Nehmen wir auch an, dass wir Double-Floating-Operationen verwenden. Dann ist h im ersten Fall ungefähr 2^{- 8}, im ersten Fall ist h ungefähr 2^{- 5}. Korrigiere diese Werte, wenn du f '' (x) oder f '' '(x) kennst.

+0

Epsilon sollte 2^-53 für das Doppelte und 2^-24 für das Float sein (das ungefähr 10^-16 und 10 ist^-7, jeweils). –

+0

Epsilon ist ** relative ** Rundungsfehler (nicht absolut). Es ist immer etwa 10^{- 16} für Double und 10^-7 für Float –

+0

Ja, ich weiß. In Ihrer Antwort sagen Sie "Epsilon - relativer Rundungsfehler, etwa 2^{- 16} für Double und 2^{- 7} für Float", was eindeutig falsch ist. Der relative (Vorwärts-) Rundungsfehler ist auch ** nicht ** immer auf dieser Skala, sondern eher der Rückwärtsfehler. Der Vorwärtsfehler kann sehr viel größer sein, wenn eine Löschung auftritt, wie es hier wahrscheinlich passieren wird. –

4

Sie wollen auf jeden Fall John Cooks Vorschlag für die Auswahl von h berücksichtigen, aber Sie möchten in der Regel keine zentrierte Differenz verwenden, um die Ableitung zu approximieren. Der Hauptgrund dafür ist, dass es eine zusätzliche Funktion Auswertung kostet, wenn Sie einen Vorwärts-Differenz zu verwenden, dh

f'(x) = (f(x+h) - f(x))/h 

Dann werden Sie Wert von f erhalten (x) kostenlos, weil Sie es bereits für Newton berechnen müssen Methode. Dies ist nicht so eine große Sache, wenn Sie eine skalare Gleichung haben, aber wenn x ein Vektor ist, dann ist f '(x) eine Matrix (die Jacobi), und Sie müssen n zusätzliche Funktionsauswertungen durchführen, um es anzunähern Verwenden des zentrierten Differenzansatzes.

1

In der Regel beeinflusst das Signalrauschen die Ableitungsqualität mehr als alles andere. Wenn Sie Rauschen in Ihrem f (x) haben, ist Savtizky-Golay ein exzellenter Glättungsalgorithmus, der oft verwendet wird, um gute Ableitungen zu berechnen. Kurz gesagt, SG passt ein Polynom lokal an Ihre Daten an, das Polynom kann dann zur Berechnung der Ableitung verwendet werden.

Paul

Verwandte Themen