next up previous contents
Next: 2.3.3 Fl�cheninhalt, Volumen, Schwerpunkt Up: 2.3 Integration Previous: 2.3.1 Oberfl�chenintegrale   Inhalt

2.3.2 Bivariater Romberg-Algorithmus

Nachdem wir die beiden Formen der Oberfl�chenintegrale vorgestellt haben und diese auf eine gemeinsame Form eines bivariaten Integrals gebracht haben, werden wir eine Methode zur numerischen Berechnung solcher bivariaten Integrale vorstellen. Siehe hierzu auch [H�l98] und [KU98]. Zun�chst f�hren wir dazu eine univariate Integrationsformel f�r

If := \int_a^b f(x) dx (2.57)

Abbildung 2.7: Trapezregel
Trapezrege

ein. Zur numerischen Approximation dieses Integrals k�nnen wir f st�ckweise durch Polynome approximieren und die Integrale �ber diese Polynome aufsummieren. Dazu zerlegen wir [a,b] durch die k+1 �quidistanten St�tzstellen x_0, ..., x_k in k Teilintervalle,

a=x_0 < x_1 < ... < x_{k-1} < x_k = b, x_j = a + j (b-a)/k f�r j = 0:k. (2.58)

Approximieren wir, wie in Abbildung 2.7 gezeigt, f auf jedem Teilintervall durch eine lineare Funktion, erhalten wir nach der Integration der linearen Funktionen die Trapezregel T_kf:

If ~T_kf := (b-a)/k (1/2 f(x_0) + \sum_{j=1}^{k-1} f(x_j) + 1/2 f(x_k) ). (2.59)

Der Fehler der Trapezregel T_kf zum exakten Wert l��t sich asymptotisch durch folgendes Theorem beschreiben.

Theorem 2.3 (Euler-Maclaurin-Summationsformel)   Sei \in C^{2m}([a,b]). Dann gilt

If - T_kf = \sum_{j=1}^{m-1} c_{2j} (f^{(2j-1)}(b) - f^{(2j-1)}(a)) h^{2j} + c_{2m} f^{(2m)}(t) (b - a) h^{2m} (2.60)

f�r t\in[a,b] und den von f unabh�ngigen Konstanten

c_{2j} = - p_{2j}(0)/(2j)!, (2.61)

die �ber die Bernoulli Polynome p_j definiert werden:

p_0(x) = 1, p'_j(x) = j p_{j-1}(x) f�r j \in N. (2.62)

Proof. Siehe [H�l98].

Nach (2.60) kann der Fehler der Trapezregel durch

If - T_kf = K_2 h^2+ K_4 h^4 + K_6 h^6 + K_8 h^8 + ... + K_{2m} h^{2m} (2.63)

dargestellt werden. Dabei ist m \in N beliebig, solange f \in C^{2m}([a,b]). Die Konstanten K_{2j} h�ngen dabei nicht von der Schrittweite h, sondern nur von f selbst ab. Damit lassen sich einzelne Fehlerterme �ber die Richardson-Extrapolation eliminieren,

T_i^0 := T_{2^ik} f, T_i^{j+1} := (4^{j+1} T_i^j - T_{i-1}^j)/(4^{j+1} - 1)    mit j = 0:i    f�r i \in N_0. (2.64)

Mit der Trapezregel berechnen wir also zun�chst T_i^0, wobei wir mit k Punkten anfangen und diese pro Schritt in i verdoppeln. Pro Schritt in j verschwindet in T_i^{j+1} jeweils ein Fehlerterm durch Kombination von T_i^j und T_{i-1}^j. Dieser Algorithmus wird auch als Romberg-Algorithmus bezeichnet. Er l��t sich auch durch das folgende Dreieckschema darstellen:

T_1^0; T_2^0, T_1^1; T_3^0, T_2^1 ...    

Aus diesem Algorithmus werden wir ein numerisches Integrationsverfahren f�r bivariate Integrale

If := \iint_A f(x,y) d(x,y) (2.65)

herleiten. F�r unsere Anwendung gen�gt es anzunehmen, da� A ein zwei-dimensionales Intervall [a,b] x [c,d] in R^2 ist. Wir unterteilen beide Intervalle [a,b] und [c,d] in jeweils k_x bzw. k_y Teilintervalle und wenden die Trapezregel auf beide Integrale getrennt an. Mit h_x := (b-a)/(k_x) und h_y := (d-c)/k_y ergibt dies die numerische Integrationsformel
numerische Integrationsformel

Abbildung 2.8: Bivariate Trapezregel
Bivariate Trapezregel

Damit erhalten wir eine bivariate Integrationsformel aus dem Tensorprodukt der univariaten Trapezregel. Wir legen dazu ein zwei-dimensionales Gitter mit den Schrittweiten h_x, h_y �ber das rechteckige Integrationsgebiet und multiplizieren die Funktionswerte an den Knoten mit bestimmten Faktoren. Dabei m�ssen wir zwischen Eckknoten, Randknoten und inneren Knoten unterscheiden. Die Faktoren aus der Formel (2.66) sind in Abbildung 2.8 dargestellt.

Den Fehler f�r (2.66) k�nnen wir durch zweimaliges Anwenden von Theorem 2.3 bestimmen. Danach gilt f�r h = h_x = h_y

If -\int_a^b F(x) dx = \sum_{j=1}^{m-1} ... dx (d - c) h^{2m}, (2.67)
\int_a^b F(x) dx - T_k F = \sum_{j=1}^{m-1} c_{2j} (F^{(2j-1)}(b) - F^{(2j-1)}(a)) h^{2j} + c_{2m} F^{(2m)}(t_x) (b - a) h^{2m} (2.68)

mit t_x \in [a,b] und t_y \in [c,d]. Uns interessiert aber nur die allgemeine Form des Fehlers und dieser enth�lt wieder nur Terme mit h^{2j},

If - T_k F = If - \int_a^b F + \int_a^b F - T_k F = K_2 h^2 + K_4 h^4 + K_6 h^6 + ... + K_{2m} h^{2m} (2.69)

Dabei h�ngen die K_{2j} nur von f ab. Dies entspricht (2.63) und wir k�nnen damit den Romberg-Algorithmus (2.64) auch f�r die bivariate Trapezregel verwenden.

Nun k�nnen wir den allgemeinen Algorithmus zur Berechnung eines Integrals einer G-Spline-Funktion �ber einer G-Spline-Fl�che der Form

\sum_{i\in I} \iint_{[0,1]^2} F_i (2.70)

angeben. Wir gehen davon aus, da� die Fl�che bereits �ber Bézierfl�chenst�cke S_i f�r i \in I gegeben ist. Die Funktion wird auf jedem Bézierfl�chenst�ck S_i durch F_i repr�sentiert. F_i kann dabei wie oben beschrieben �ber ein Bézierkontrollnetz gegeben sein. Es kann aber auch nur von der Oberfl�che selbst abh�ngen. Wir m�ssen nur jeweils die Trapezregel, bzw. die in der Trapezregel verwendete Funktion entsprechend w�hlen. F�r die bisher vorgestellten Integrale geht dies aus den Formeln (2.51) und (2.56) hervor.



Algorithmus 2.1  
integrate
Bivariater Romberg-Algorithmus
1.
Setze das Ergebnis value der Integration auf 0.
2.
Berechne das Integral f�r jedes Bézierkontrollnetz bgs der G-Spline-Fl�che separat �ber den Romberg-Algorithmus:
(a)
Initialisiere die Variablen f�r den Romberg-Algorithmus. Dabei ist t ein Array der L�nge max_iter, welche die maximale Anzahl der Iterationen angibt. n gibt die Anzahl der St�tzstellen in einer Gitterrichtung an und i wird als Iterationsz�hler auf 1 gesetzt.
(b)
Berechne die Trapezregel f�r die Intervalll�nge 1/(n-1) pro Gitterrichtung ( n^2 St�tzstellen) und lege das Ergebnis in t[0] ab.
(c)
Wiederhole Folgendes, solange die gew�nschte Toleranz prec nicht erreicht wurde oder die maximale Anzahl der Iterationen erreicht wurde:
i.
Verdopple n.
ii.
Berechne die Trapezregel f�r die Intervalll�nge 1/(n-1) pro Gitterrichtung ( n^2 St�tzstellen) und lege das Ergebnis in t[i] ab.
iii.
F�r j = i:-1:1 berechne die Romberg-Extrapolation iterativ �ber t[j-1] = (4 ^(j-i+1) * t[j] - t[j-1]) / (4 ^(j-i+1) - 1).
iv.
Zur Bestimmung der Genauigkeit berechne den Betrag der Differenz zwischen dem besten Wert der Extrapolation in t[0] und dem besten Ergebnis der vorherigen Iteration.
v.
Erh�he i um 1.
(d)
Addiere das Ergebnis des Romberg-Algorithmus' zu value.
3.
Das Ergebnis der Integration steht in value.



Es ist zun�chst zu beachten, da� wir das Array t von rechts nach links f�llen. Damit ben�tigen wir nur ein einziges, eindimensionales Array f�r die Extrapolation, in dem aus den alten Werten die neuen berechnet werden. Im ersten Schritt wird das Ergebnis der Trapezregel in ein noch nicht benutztes Element t[i] des Arrays abgelegt und die Extrapolationswerte von diesem Element aus bis t[0] eingetragen.

Die im Algorithmus verwendete Toleranz prec bezieht sich dabei auf das Integral �ber eine einzelne Bézierfl�che. Durch das Aufsummieren der Ergebnisse verringert sich die Toleranz f�r das gesamte Integral. Bei M Bézierfl�chen ist sie nur noch M *   prec.

Die verwendete Trapezregel entspricht nat�rlich der bivariaten Trapezregel (2.66), wobei im Algorithmus die Anzahl der Intervalle k durch die Anzahl der St�tzstellen n ersetzt wurde. Es gilt n = k+1. Je nach Integral kann f durch andere Funktionen ersetzt werden. So verwenden wir z.B. f�r das Oberfl�chenintegral von Vektorfeldern die Funktion <f,n>. Spezielle Integrale stellen wir weiter unten vor.

Abbildung 2.9: Funktionsauswertungen bei der bivariaten Trapezregel
Funktionsauswertungen bei der bivariaten Trapezregel

Weiter ist auch zu beachten, da� wir nach dem ersten Auswerten der Trapezregel f�r die folgenden Auswertungen die Funktion nicht mehr an allen St�tzpunkten berechnen m�ssen. Nachdem wir n jeweils verdoppeln, ben�tigen wir nur noch die St�tzstellen (i 1/(n-1),j 1/(n-1)), bei denen wenigstens eine der beiden Zahlen i und j ungerade ist. F�r die restlichen St�tzstellen k�nnen wir das Ergebnis der vorherigen Trapezregel mit 1/4 multiplizieren. Die Punkte, an denen die Funktion noch ausgewertet werden mu�, sind in Abbildung 2.9 gekennzeichnet. Von den n^2 St�tzstellen wird die Funktion dann nur noch an n^2 - (n/2)^2 = 3/4n^2 St�tzstellen berechnet.

Eine m�gliche Verbesserung des Romberg-Algorithmus' w�re eine adaptive Anpassung an den Approximationsfehler. Wir berechnen dazu zun�chst das Integral f�r ein zwei-dimensionales Integrationsintervall A und sch�tzen zus�tzlich den Fehler des Ergebnisses ab. Ist dieser Fehler gr��er als eine vorgegebene Toleranz s, unterteilen wir das Intervall A in vier gleichgro�e Teilintervalle und berechnen f�r diese Intervalle das Integral getrennt mit der Toleranz s/4. Dies wird rekursiv wiederholt, bis die Fehlerabsch�tzung kleiner als die Toleranz ist. Hierdurch wird die Anzahl der St�tzstellen dem jeweiligen Approximationsfehler auf dem Teilintervall angepa�t. Nachdem wir aber schon das Integrationsgebiet am Anfang in die einzelnen quadratischen Bézierfl�che unterteilen, haben wir den adaptiven Algorithmus nicht implementiert.

next up previous contents
Next: 2.3.3 Fl�cheninhalt, Volumen, Schwerpunkt Up: 2.3 Integration Previous: 2.3.1 Oberfl�chenintegrale   Inhalt
Copyright © 1999-2002 Frank C. Langbein. All rights reserved.
Permission is granted to copy, distribute and/or modify this document under the terms of the GNU Free Documentation License, Version 1.1 or any later version published by the Free Software Foundation.

Contact: webmaster@langbein.org
URI: http://www.langbein.org/fileadmin/research/surfaces/diploma/HTML/node21.html