[WS] Extrapolation

25.09.2008, 22:01

tigerbine

[WS] Extrapolation

Gliederung

Allgemeine Idee (Richardson Extrapolation)

Neville-Schema

Romberg-Integration mit Trapezregel

Differenzieren mit h-Methode - Teil 1

Differenzieren mit der h-Methode - Teil 2

L'Hospital

25.09.2008, 22:10

tigerbine

Auf diesen Beitrag antworten »

Allgemeine Idee (Richardson Extrapolation)
Gegeben ist ein zu berechnender Wert einer Funktion. Nennen wir in Allgemein F. Das kann zum Beispiel ein bestimmes Integral sein. In Abhängigkeit von der Schrittweite h haben wir eine Näherungsformel $\begin{eqnarray*} F_h \end{eqnarray*}$ .

Nun treffen wir die Annahme, dass der Fehler eine bestimmte Gestalt hat.

$\begin{eqnarray*} F_h(f)-F(f) = c_1h^{k_1} + c_2h^{k_2} + c_3h^{k_3} + c_4h^{k_4} + ...~,\quad k_1 < k_2 \end{eqnarray*}$

Dabei habe diese Formel die Ordnung k1, d.h. wir können den Approximationsfehler mit $\begin{eqnarray*} O\left(h^{k_1}\right) \end{eqnarray*}$ beziffern (Erläuterung,). Dabei steht dies verkürzend für eine Summe vonTermen der Ordnung $\begin{eqnarray*} h^{k_1} \end{eqnarray*}$ und höher. (Siehe PDF im Anhang.)

Nun nehmen wir eine zweite Formel, mit der halben Schrittweite

$\begin{eqnarray*} F_{h/2}(f)-F(f) = c_1\left(\frac{1}{2}h\right)^{k_1} + c_2\left(\frac{1}{2}h\right)^{k_2} + c_3\left(\frac{1}{2}h\right)^{k_3} + c_4\left(\frac{1}{2}h\right)^{k_4} + ...~,\quad k_1 < k_2 \end{eqnarray*}$ (*)

Aus diesen beiden konstruiert man nun eine neue Formel.

$\begin{eqnarray*} \hat{F}_{h/2} &&=\frac{2^{k_1}F_{h/2}-F_h}{2^{k_1}-1} \quad(**) \end{eqnarray*}$

Welche Ordnung die neue Formel?

$\begin{eqnarray*} \frac{2^{k_1}F_{h/2}-F_h}{2^{k_1}-1} && = \frac{2^{k_1}\left[c_1\left(\frac{1}{2}h\right)^{k_1} + c_2\left(\frac{1}{2}h\right)^{k_2} + c_3\left(\frac{1}{2}h\right)^{k_3} + c_4\left(\frac{1}{2}h\right)^{k_4}+ \ldots\right]-\left[c_1h^{k_1} + c_2h^{k_2} + c_3h^{k_3} + c_4h^{k_4} + \ldots\right]} {2^{k_1}-1} \\&& =\frac{2^{k_1}\left[c_2(\frac{1}{2}h)^{k_2} + c_3\left(\frac{1}{2}h\right)^{k_3} + c_4\left(\frac{1}{2}h\right)^{k_4}+ \ldots\right]-\left[ c_2h^{k_2} + c_3h^{k_3} + c_4h^{k_4} + \ldots\right]} {2^{k_1}-1} \\&&=d_2h^{k_2} + d_3h^{k_3} + d_4h^{k_4} + ... \end{eqnarray*}$

Hier können wir nun den Approximationsfehler mit $\begin{eqnarray*} O\left(h^{k_2}\right) \end{eqnarray*}$ beziffern, die Formel hat also in obigem Sinne eine höherer Ordnung.

Was bedeutet dies?

Wählen wir als konkretes Beispiel $\begin{eqnarray*} k_1=2 \end{eqnarray*}$ und $\begin{eqnarray*} k_2=4 \end{eqnarray*}$ . Der Approximationsfehler der ersten Formel hat dann die Größenordnung $\begin{eqnarray*} O(h^2) \end{eqnarray*}$ , das bedeutet, halbiert man h, so viertelt sich der maximale Fehler dieses Verfahrens. Die durch den Quotienten gewonne Formel leistet, da ihr Approximationsfehler von der Größe $\begin{eqnarray*} O\left(h^4\right) \end{eqnarray*}$ bei der Halbierung von h schon eine Sechzehntelung.

Wie kann man es nutzen?

Da es nicht immer möglich ist, die benötigte Schrittweite h a priori zu bestimmen, wird man sich bei der Umsetzung von Verfahren für ein Abbruch Kriterium entscheiden müssen. Welche Möglichkeiten bieten sich einem? Zunächst einmal könnte man mit einer Schrittweite h begonnen, den Wert der Näherung einer Prüfung unterziehen und falls sie nicht bestanden wird, die Berechnung mit der halbierten Schrittweite wiederholen, so lange bis die Abbruchbedingung erfüllt ist.

Sollte der Algorithmus nicht schon direkt für erst gewähltes h abbrechen, so wird man sicher zunächst einen zweiten Durchlauf starten. Wenn immer noch nicht abgebrochen wird, so bietet uns obige Formel (*) aus den bereits vorhandenen Daten (sofern man sie gespeichert hat) mit geringem Aufwand eine bessere Approximation zu bestimmen.

Diese Idee kann man in einem Dreiecksschema umsetzten und wir werden im folgenden den Bezug zum Neville-Schema, welches im [WS] Polynominterpolation-Theorie vorgestellt wird aufzeigen.

25.09.2008, 22:11

tigerbine

Auf diesen Beitrag antworten »

Neville Schema
Wie hängen nun die obige Formel und die im Nevilleschema verwendete Formel zusammen? Zunächst bestimmt man sich die Näherungen der gesuchten Größe bezüglich h und h/2. Da man sich am Ende für den Wert für h=0 interessiert, könnte man durch diese beiden Punkte ein Interpolationspolynom legen und bei 0 auswerten.

Dies soll mit dem "klassischen" Neville-Schema geschehen. ([WS] Polynominterpolation-Beispiele) Wir haben die Daten:

$\begin{eqnarray*} (h /F(h)=p_{0,0}(x)=p_{0,0}(0)),~(0.5h / F(0.5h)=p_{1,0}(x)=p_{1,0}(0)) \end{eqnarray*}$

Dadurch legen wir ein Interpolationspolynom wie folgt

$\begin{eqnarray*} p_{0,1}(x) &&= \frac{(x-h)\cdot p_{1,0}(x) - (x-0.5h)\cdot p_{0,0}(x)}{0.5h-h} \end{eqnarray*}$

Ausgewertet werden soll bei x=0. Also:

$\begin{eqnarray*} p_{0,1}(0) &&=\frac{(-h)\cdot F(0.5h)- (-0.5h)\cdot F(h)}{h(0.5-1)}\\ &&=\frac{(-h)\cdot F(0.5h)- 0.5(-h)\cdot F(h)}{(-h)(1-0.5)}\\ &&=\frac{0.5(-h)\cdot( 2^1F(0.5h)- F(h))}{0.5(-h)(2^1-1)}\\&&=\frac{2^1F(0.5h)- F(h)}{2^1-1} \end{eqnarray*}$

Das entspricht also der Formel (**) für $\begin{eqnarray*} k_1=1 \end{eqnarray*}$

26.09.2008, 17:50

tigerbine

Auf diesen Beitrag antworten »

Romberg-Verfahren
Mit der Euler-Maclaurinschen Summenformel kann man zeigen, dass die Trapezsumme T zur Maschenweite h folgende asymptotisch Entwicklung besitzt:

$\begin{eqnarray*} T(h) = I(f) + r_2h^2+r_4h^4+...+r_{2m}h^{2m} + O(h^{2m+2}) \end{eqnarray*}$

Dabei ist f eine 2m+2-mal stetig differenzierbare Funktion. Es treten hier also nur gerade Potenzen auf. Es ergeben sich also in Anlehnung an obige Notation die Potenzen:

$\begin{eqnarray*} k_1=2 \cdot 1,~k_2=2\cdot 2,~k_3=2\cdot 3,~...,k_n=2\cdot n \end{eqnarray*}$

Somit können wir obige Rekursionsformel (**) wiederholen und in ein Schema analog zur Nevilledarstellung eintragen.

$\begin{eqnarray*} h/1 \quad ||~ y_0 = T_{0,0} \quad | ~T_{0,1} \quad |~ T_{0,2} \quad |~T_{0,3} \quad | \end{eqnarray*}$

$\begin{eqnarray*} h/2 \quad ||~ y_1 = T_{1,0} \quad | ~T_{1,1} \quad |~ T_{1,2} \quad | \end{eqnarray*}$

$\begin{eqnarray*} h/4 \quad ||~ y_2 = T_{2,0} \quad | ~T_{2,1} \quad | \end{eqnarray*}$

$\begin{eqnarray*} h/8 \quad ||~ y_3 = T_{3,0} \quad | \end{eqnarray*}$

$\begin{eqnarray*} T_{i,k}=\frac{2^{2k}T_{i+1,k-1}-T_{i,k-1}}{2^{2k}-1} \end{eqnarray*}$

Wie kann man das wieder in Bezug zu den Interpolationspolynomen stellen?

Wir ersetzen x:=h² und erhalten somit

$\begin{eqnarray*} x_0=1,~x_1=\frac{1}{4},~x_2=\frac{1}{16},~x_3=\frac{1}{64} \end{eqnarray*}$

Dann kann man auch wieder das klassische Neville-Schema anwenden, die interpolierenden Polynome sind gerade in h. Zu beachten ist, dass man die summierte Trapezregel jedoch für die Schrittweiten 1, 1/2, 1/4, 1/8 berechnet!

Beispielrechnung für $\begin{eqnarray*} \int_1^2~\frac{e^x}{x}~dx \end{eqnarray*}$

code:

1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:

Es wird ein Integral durch Romberg-Extrapolation approximiert.
 
Funktion in rombergf.m anlegen
Intergralanfang eingeben. a = 1
Intergralende eingeben.   b = 2
Genauigkeit           eps = 0.01
 
--------------------------------------------------------------------------------------------
 
I_Matlab=3.059117e+000 
 
Neville-Romberg Schema  - Funktionswerte bei h=0
===============================================
 
IR =
    1.0000    3.2064    3.0607    3.0591    3.0591
    0.5000    3.0971    3.0592    3.0591         0
    0.2500    3.0687    3.0591         0         0
    0.1250    3.0615         0         0         0
 
Fehler = 0.000000

29.09.2008, 16:46

tigerbine

Auf diesen Beitrag antworten »

Differenzieren mit der h-Methode - Teil 1
Sei die Funktion im Punkt $\begin{eqnarray*} x_0 \end{eqnarray*}$ differenzierbar, so existiert der Grenzwert:

$\begin{eqnarray*} \lim_{h \to 0}~ \frac{f(x_0+h)-f(x_0)}{h} \end{eqnarray*}$

Dabei ist der Bruch der Differenzenquotient und man kann ihn als Steigung der Sekanten durch die Punkte $\begin{eqnarray*} (x_0,f(x_0)) (x_0+h,f(x_0+h)) \end{eqnarray*}$ auffassen. Diese können wir also als Funktion in h auffassen und interessieren uns für den Wert an der Stelle h=0.

$\begin{eqnarray*} S(h) = \frac{f(x_0+h)-f(x_0)}{(x_0+h)-x_0} =\frac{f(x_0+h)-f(x_0)}{h} \end{eqnarray*}$

Schnell ist also die Idee formuliert. Man bestimmt für verschiedene Werte von h den Funktionswert von S und legt durch die Daten das Interpolationspolynom und wertet es an der Stelle 0 aus. Dabei ist es sinnvoll die Werte durch eine Nullfolge zu bestimmen, zum Beispiel wieder durch die Rombergfolge:

$\begin{eqnarray*} h_0=1, ~h_n=\frac{1}{2^n} \end{eqnarray*}$

Als Beispiel wollen wir uns die Funktion $\begin{eqnarray*} f(x)=\sin(x) \end{eqnarray*}$ anschauen und die Ableitung im Punkt $\begin{eqnarray*} x_0=0 \end{eqnarray*}$ approximieren. Theoretisch würden wir erhalten $\begin{eqnarray*} f'(0)=\cos(0)=1 \end{eqnarray*}$ . Dann lautet die Funktion S:

$\begin{eqnarray*} S(h) = \frac{\sin(0+h)-\sin(0)}{(0+h)-0} = \frac{\sin(h)-0}{h} = \frac{\sin(h)}{h} \end{eqnarray*}$

Wählen wir die Knoten (bzgl. h), so erhalten wir den Startdatensatz:

$\begin{eqnarray*} h_0 =\frac{1}{8},~h_1=\frac{1}{16},~h_2=\frac{1}{32} \end{eqnarray*}$

$\begin{eqnarray*} S(h_0) =8\cdot \sin(\frac{1}{8}),~S(h_1)=16\cdot \sin(\frac{1}{16}),~S(h_2)=32 \cdot \sin(\frac{1}{32}) \end{eqnarray*}$

code:
1: 2: 3: 4: 5: 6: 7: 8: 9: 10:	Neville Schema - Funktionswerte bei 0 ===================================== NW = 0.1250 0.9974 1.0013 1.0004 0.0625 0.9993 1.0007 0 0.0313 1.0000 0 0

[attach]8737[/attach]

29.09.2008, 18:15

tigerbine

Auf diesen Beitrag antworten »

Differenzieren mit der h-Methode - Teil 2
Ein alternativer Weg wäre die Verwendung des zentralen Differenzenquotienten. Wieder erkennen wir die zugehörige Sekantensteigung und im Grenzwertfall die Ableitung.

$\begin{eqnarray*} Z(h)=\frac{f(x_0+h) - f(x_0-h)}{(x_0+h)-(x_0-h)}= \frac{f(x_0+h) - f(x_0-h)}{2h} \end{eqnarray*}$

Was sind die Vorteile dieser Darstellung?

Untersuchen wir die Funktion Z ein wenig. Es ist:

$\begin{eqnarray*} Z(-h)=\frac{f(x_0+(-h)) - f(x_0-(-h))}{2(-h)} =\frac{f(x_0-h) - f(x_0+h)}{-2h} =\frac{f(x_0+h) - f(x_0-h)}{2h} = Z(h) \end{eqnarray*}$

D.h. Z ist eine Gerade Funktion in h. Daher macht es Sinn, auch nur gerade Polynome (nur gerade Potenzen) zur Interpolation zu verwenden.
(Die entspricht dem Vorgehen, welches wir schon bei der Romberg-Integration gesehen haben)

Wenden wir uns wieder dem Beispiel zu. Der Sinus ist eine ungerade Funktion, daher lautet die Funktion:

$\begin{eqnarray*} Z(h) = \frac{\sin(0+h)-\sin(0-h)}{2h} = \frac{\sin(h)-\sin(-h)}{2h} =\frac{\sin(h)+\sin(h)}{2h} = \frac{\sin(h)}{h} \end{eqnarray*}$

(Dass sich S und Z hier gleichen liegt daran, dass die Funktion f hier ungerade ist)

Die Interpolationspolynome haben dann bezüglich der Variablen h die Gestalt:

$\begin{eqnarray*} p_{i,k}(h) = \sum_{j=0}^{k}~a^{(i)}_kh^{2k} \end{eqnarray*}$

Um das Neville-Schema übernehmen zu können, schreiben wir mit der Substitution $\begin{eqnarray*} x:=h^2 \end{eqnarray*}$

$\begin{eqnarray*} p_{i,k}(x) = \sum_{j=0}^{k}~a^{(i)}_kx^{k} \end{eqnarray*}$

Wieder wählen wir die Knoten (bzgl. h), so erhalten wir den Startdatensatz:

$\begin{eqnarray*} h_0 =\frac{1}{8},~h_1=\frac{1}{16},~h_2=\frac{1}{32} \end{eqnarray*}$

$\begin{eqnarray*} Z(h_0) =8\cdot \sin(\frac{1}{8}),~Z(h_1)=16\cdot \sin(\frac{1}{16}),~Z(h_1)=32 \cdot \sin(\frac{1}{32}) \end{eqnarray*}$

Das ergibt bzgl. x (beachte, die Funktionswerte werden beibehalten!)

$\begin{eqnarray*} x_0 =\frac{1}{64},~x_1=\frac{1}{256},~x_2=\frac{1}{1024} \end{eqnarray*}$

Das Neville Schema liefert bzgl. x:

code:
1: 2: 3: 4: 5: 6: 7: 8: 9: 10:	Neville Schema - Funktionwerte bei 0 ===================================== NW = 0.0156 0.9974 1.0000 1.0000 0.0039 0.9993 1.0000 0 0.0010 0.9998 0 0

Über die Wahl einer Abbruchbedingung kann man sich im anhängenden PDF informieren. Sprich wann liefern kleine h Werte praktisch aufgrund von Rechenungenauigkeiten (Maschinenzahlen) keine besseren Ergebnisse mehr.

02.10.2008, 18:16

tigerbine

Auf diesen Beitrag antworten »

L'Hospital
Gesucht ist:

$\begin{eqnarray*} \lim_{x \to 0+} \frac{\cos(x)-1}{\sin(x)} ~(=0) \end{eqnarray*}$

Man setzt:

$\begin{eqnarray*} L(x):= \frac{\cos(x)-1}{\sin(x)} \end{eqnarray*}$

Und bestimmt einen Startdatensatz:

$\begin{eqnarray*} x_0=\frac{1}{8},~x_0=\frac{1}{16},x_0=\frac{1}{32} \end{eqnarray*}$

$\begin{eqnarray*} L(\frac{1}{8}) = \frac{\cos(\frac{1}{8})-1}{\sin(\frac{1}{8})},~L(\frac{1}{16}) = \frac{\cos(\frac{1}{16})-1}{\sin(\frac{1}{16})},~L(\frac{1}{32}) = \frac{\cos(\frac{1}{32})-1}{\sin(\frac{1}{32})} \end{eqnarray*}$

Das Neville-Schema liefert:

code:
1: 2: 3: 4: 5: 6: 7: 8: 9:	Neville Schema - Funktionswerte bei 0 ===================================== NW = 0.1250 -0.0626 0.0001 -0.0000 0.0625 -0.0313 0.0000 0 0.0313 -0.0156 0 0

Neue Frage »

Antworten »

[WS] Extrapolation

Verwandte Themen