Romberg-Integration

03.09.2008, 21:47

tigerbine

Romberg-Integration

So, nun ein neues Kapitel, aber immer noch auf der Suche nach dem heiligen Integral. Die Idee ist in den meisten Büchern auch schnell beschrieben. Man nehme eine summierte Newton-Cotes-Formel, gerne die Trapezregel, und bestimme zu verschiedenen Maschenweiten h (oder eben immer mehr Teilintervallen), den Näherungswert. Sollte die zu integrierende Funktion dann nur genügend oft stetig differenzierbar sein, konvergiert die Integrationsformel für h->0 gegen das Integral.

Fragen zur Euler-Maclaurenschen Summenformel werde ich später stellen.

Es scheint auch einsichtig die Maschenweiten durch Halbierung zu verfeinern. Diese wird dann als Romberg-Folge bezeichnet:

$\begin{eqnarray*} h_i:=\frac{1}{2^i} \end{eqnarray*}$

Bezüglich dieser bestimmt man dann das Ergebnis der Summierten Trapezregel. Als Beispielintegral sei gewählt:

$\begin{eqnarray*} \int_1^2~\frac{1}{x}~dx = \ln(2) \approx 0.693147 \end{eqnarray*}$

Bestimmt man die ersten Näherungen, so erhält man:

$\begin{eqnarray*} T(1) ~= 0.750000 \end{eqnarray*}$

$\begin{eqnarray*} T(\frac{1}{2})~=0.708333 \end{eqnarray*}$

$\begin{eqnarray*} T(\frac{1}{4})~= 0.697024 \end{eqnarray*}$

$\begin{eqnarray*} T(\frac{1}{8})~= 0.694122 \end{eqnarray*}$

$\begin{eqnarray*} T(\frac{1}{16})= 0.693391 \end{eqnarray*}$

Nun könnte man h sicher immer weiter verkleinern. Die Extrapolationsverfahren verfolgen allerdings einen anderen Ansatz. Man legt durch die Punkte (Maschenweite/Näherungswert) ein Interpolationspolynom und wertet dieses an der Stelle x=0 aus. Bei der Umsatzung müssen dann zwei Fragen geklärt werden:

Wie bestimmt man das IP, die Knoten liegen ja in der falschen Reihenfolge vor und können auch nicht umsortiert werden. Eigentlich sollte es egal sein, aber siehe unten, ich erhalte immer etwas anderes als in einem vorgerechneten Beispiel.
Wieso bringt dieser Ansatz überhaupt erfolgt, bei der Polynominterpolation wurde ja gesagt, dass sie eigentlich nicht zu Extrapolation geeignet ist. ([WS] Polynominterpolation-Theorie)

Zur Frage 1. Aus dem Workshop Polynominterpolation ist das Neville-Schema zur Bestimmung des IPs bekannt. Da hier nur der Funktionswert an einer bestimmten Stelle, nämlich x=0, interessiert, können auch die konkreten Funktionswerte gleich eingetragen werden.

$\begin{eqnarray*} p_{i,k}(x) = \frac{(x-x_i)\cdot p_{i+1,k-1}(x) - (x-x_{i+k})\cdot p_{i,k-1}(x)}{x_{i+k}-x_i} \end{eqnarray*}$

$\begin{eqnarray*} x_0 \quad ||~ y_0 = p_{0,0} \quad | ~p_{0,1} \quad |~ p_{0,2} \quad |~p_{0,3} \quad | \end{eqnarray*}$

$\begin{eqnarray*} x_1 \quad ||~ y_1 = p_{1,0} \quad | ~p_{1,1} \quad |~ p_{1,2} \quad | \end{eqnarray*}$

$\begin{eqnarray*} x_2 \quad ||~ y_2 = p_{2,0} \quad | ~p_{2,1} \quad | \end{eqnarray*}$

$\begin{eqnarray*} x_3 \quad ||~ y_3 = p_{3,0} \quad | \end{eqnarray*}$

Im ersten Schritt bestimmt man also die IPs durch 2 benachbarte Knoten und wertet diese in x=0 aus. Man erhält:

Nun komme ich hier schon nicht auf die Werte wie in diesem PDF - Seite 3. Im Beispiel weicht der Wert für h=0.5 schon ab. Es könnte auch einfach eine 0 in der Datei fehlen oder ich rechne hier mal wieder Mist zusammen. Wäre super, wenn das mal jemand mal nachrechnen könnte.

Danke Wink

04.09.2008, 00:02

tigerbine

Auf diesen Beitrag antworten »

Falls der viele Text abschreckt, wäre jemand so nett mir mit der Summierten Trapezregel (2 Teilintervalle) die Näherung für

$\begin{eqnarray*} \int_1^2~\frac{1}{x}~dx \end{eqnarray*}$

zu bestimmen. Augenzwinkern

04.09.2008, 14:28

Huggy

Auf diesen Beitrag antworten »

Seltsame Bitte! Die hast du doch oben schon hingeschrieben.

04.09.2008, 18:35

tigerbine

Auf diesen Beitrag antworten »

Was ist seltsam? Ja ich habe das schon berechnet, aber ich komme doch nicht auf den gleichen Wert wie in der PDF Datei.

ICH erhalte:
$\begin{eqnarray*} T(\frac{1}{2})~=0.708333 \end{eqnarray*}$

PDF enthält:
$\begin{eqnarray*} T(\frac{1}{2})~=0.78333 \end{eqnarray*}$

Wer hat nun Recht?

04.09.2008, 19:52

Huggy

Auf diesen Beitrag antworten »

Du hast recht. Das kann in der PDF-Datei nur ein Schreibfehler sein.

04.09.2008, 20:30

tigerbine

Auf diesen Beitrag antworten »

Ok, dann kann ich ja weiter rechnen.

15.09.2008, 02:43

tigerbine

Auf diesen Beitrag antworten »

so, ich habe weitergerechnet. Augenzwinkern

Verschiedene Bücher nähern sich dem Schema auf verschiedene Arten. Sie enden aber alle beim gleichen Schaubild. In Anlehungen an das Neville Schema, habe ich mich für eine obere Dreiecksmatrix entschieden, aber es ist imho Sache des Betrachters, wie er die weiteren Indizes vergibt. Haupstsache die Berechnungsformel wird entsprechend angepasst. Für obiges Beispiel ergibt sich:

Zitat:

Es wird ein Integral durch Romberg-Extrapolation approximiert.

Beachte: Der Datensatz hat die Form
Knoten: h_0,...,h_n
Integralwerte: i_0,...,i_n
Iterationen: n
Funktion in rombergf.m anlegen

--------------------------------------------------------------------------------------------

I_Matlab=3.059117e+000

Neville-Romberg Schema - Funktionwerte bei h=0
=====================================

IR =
1.00000000000000 3.20640493896219 3.06066345535987 3.05914424200495 3.05911683681869
0.50000000000000 3.09709882626045 3.05923919283964 3.05911726502473 0
0.25000000000000 3.06870410119484 3.05912488551316 0 0
0.12500000000000 3.06151968943358 0 0 0

Fehler = 0.000000

In der ersten Spalte steht die Maschenweite, in der zweiten der Wert der zugehörigen summierten Trapezregel. Diese kann man nun in ein h-I - Koordinatensystem eintragen. Man könnte nun nach dem Schema von Neville die IPs durch diese Punkte bestimmen und bei h=0 auswerten und in das Schema eintragen. Das wird jedoch bei romberg nicht gemacht.

Warum? Die Idee wird ja beibehalten, nur werden wohl nur Polynome mit geraden Exponenten verwendet. Die Daten sollten dennoch reichen, da man ja die Koeffizienten der ungeraden Potenzen einfach 0 setzten kann. Nur wiso wählt man gerade Polynome? Es wird immer wieder auf die Euler Meclaurin-Formel angespielt, aber irgendwie check ich es nicht, was die hier mit dem IP zu tun hat.

Wer kann helfen?

Danke,
tigerbine Wink

15.09.2008, 17:51

tigerbine

Auf diesen Beitrag antworten »

Mehr information
Mit der Euler-Maclaurinschen Summenformel zeigt man, dass die Trapezsumme T zur Maschenweite h folgende asymptotisch Entwicklung besitzt:

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

Dabei ist f eine 2m+2-mal stetig differenzierbare Funktion.

Da I(f) - Das gesuchte Integral- ja unabhängig von h ist, also eine Konstante , ist der der Wert T(h) den man nun gegenüber h in das Koordinatensystem einträgt, der Funktionswert eines Polynoms in h, mit lauter geraden Koeffizienten.

Legt man deswegen durch die Punkte ein IP mit nur geraden Exponenten?

17.09.2008, 20:09

tigerbine

Auf diesen Beitrag antworten »

RE: Mehr information
push... Ups