Numerische Integration

Bei vielen Problemen des naturwissenschaftlichen Rechnens treten Integrale auf, die nicht in expliziter Form dargestellt werden können, sei es, daß kein geschlossener Ausdruck für eine Stammfunktion existiert (z.B. òe^-x2 dx), oder daß der Integrand selbst nicht in geschlossener Form bekannt ist (z.B. weil er nur an diskreten Stellen aus Messungen oder Simulationsrechnungen bestimmt wurde). In all diesen Fällen ist man auf Näherungsverfahren angewiesen. Unter numerischer Integration versteht man die näherungsweise Berechnung von bestimmten Integralen

I(f) : = ó õ b a  f(x)  dx

(1)

mit einem (endlichen) Intervall [a, b] als Integrationsbereich. Die Formeln zur näherungsweisen Berechnung von I(f) heißen Integrations- oder Quadraturformeln . Die einfachste Approximation von I(f) mit Hilfe von Riemann-Summen

ó õ b a  f(x) dx   »    (b - a) n n-1 å i=0  f(xi)

(2)

und einer äquidistanten Zerlegung x_i = a + i (b - a) / n von [a, b] konvergiert jedoch für die meisten praktischen Anwendungen zu langsam. Die Situation läßt sich z.B. durch folgende Überlegung verbessern: Der Integrand f wird an vorgegebenen Stützpunkten durch Polynome interpoliert und diese werden dann analytisch integriert. Diese Vorgangsweise führt zu sogenannten interpolatorischen Quadraturformeln:

Einfache Newton-Cotes-Formeln

Das Integrationsintervall [a, b] sei durch (n + 1) äquidistante Stützstellen

xi : = a + i h        (i = 0, 1, ¼, n)

(3)

mit x₀ = a und x_n = b in n Teilintervalle der Breite (Schrittweite) h : = (b - a) / n unterteilt. Sei P_01¼n(x) das zugehörige Interpolationspolynom, d.h.

P01¼n(xi) = f(xi) = fi        (i = 0, 1, ¼, n)   .

(4)

n = 1 : Lineares Interpolationspolynom P₀₁ durch 2 Stützpunkte (Trapezregel).

ó õ b a  f(x) dx » ó õ x1 x0  P01(x) dx =  h 2 (f0 + f1)

(5)

Diese Formel heißt Trapezregel , weil der Näherungswert von I(f) die Fläche des Trapezes mit den Ecken (x₀,0), (x₁,0), (x₁,f₁), (x₀,f₀) ist.

Abbildung 1: Trapezregel.

Abbildung 2: Simpsonregel.

n = 2 : Quadratisches Interpolationspolynom P₀₁₂ durch 3 Stützpunkte (Simpsonregel).

ó õ b a  f(x) dx » ó õ x2 x0  P012(x) dx =  h 3 (f0 + 4 f1 + f2)

(6)

Zur Herleitung der Simpsonregel wählen wir folgende Darstellung für das quadratische Interpolationspolynom P₀₁₂(x):

P012(x) = a (x - x0)2 + b (x - x0) + c   .

(7)

Mit diesem Ansatz und den Bedingungen P₀₁₂(x_i) = f_i für i = 0, 1, 2 ergibt sich ein lineares Gleichungssystem für die unbekannten Koeffizienten a, b und c:

c = f0 (x1 - x0)2 a + (x1 - x0) b + c = f1 (x2 - x0)2 a + (x2 - x0) b + c = f2   .

Wegen (x₁ - x₀) = h bzw. (x₂ - x₀) = 2 h ist das gleichbedeutend mit

h  a + b =  f1 - f0 h 2 h  a + b =  f2 - f0 2 h   ,

woraus man sofort die Lösungen für a und b abliest:

a =  f0 - 2 f1 + f2 2 h2 b =  - 3 f0 + 4 f1 - f2 2 h   .

Mit (7) erhält man damit den gewünschten Näherungswert für I(f):

ó õ b a  f(x) dx » ó õ x2 x0  P012(x) dx = ó õ x2 - x0 0  (a u2 + b u + c)  du = é ë  a 3 u3 +  b 2 u2 + c u ù û 2 h 0  =  h 3 { ( 4 f0 - 8 f1 + 4 f2) + ( -9 f0 + 12 f1 - 3 f2) + 6 f0 } =  h 3 ( f0 + 4 f1 + f2 )   .

Zusammengesetzte Newton-Cotes-Formeln

Im Prinzip sind der Erhöhung des Grades n der Interpolationspolynome keine Grenzen gesetzt. Für n = 8 und n ³ 10 werden die Quadraturformeln aber numerisch instabil. Außerdem können Interpolationspolynome höheren Grades an den Intervallenden stark oszillieren. Es ist daher besser, das Grundintervall [a, b] in N gleichgroße Teilintervalle zu unterteilen, die Newton-Cotes-Formeln auf jedes Teilintervall anzuwenden und die einzelnen Beiträge dann zu addieren:

Zusammengesetzte Trapezregel

Abbildung 3: Zusammengesetzte Trapezregel.

Für die Unterteilung x_i = x₀ + ih, (i=0,1,¼,N), von [a,b] mit x₀=a, x_N=b und h=(b-a)/N erhält man aus der Trapezregel (5) im Teilintervall [x_i,x_i+1] den Näherungswert

Ii  :=   h 2 ( fi + fi+1 )           i = 0, 1, ¼ , (N - 1)   .

Für das gesamte Intervall [a,b] ergibt sich damit der Näherungswert

ó õ b a  f(x) dx » N-1 å i=0   h 2  (fi + fi+1) =  h 2  [ (f0 + f1) + (f1 + f2) + ¼ + (fN-2 + fN-1) + (fN-1 + fN) ] =  h 2  [ f0 + 2 f1 + ¼ + 2 fN-1 + fN ]   ,

ó õ b a  f(x) dx    »    h é ë  1 2 f0 + f1 + ¼ + fN-1 +  1 2 fN ù û   =:   T(h)

(8)

Die hier auftretende Summe T(h) heißt Trapezsumme zur Schrittweite h.

Zusammengesetzte Simpsonregel

Für eine gerade Anzahl N von Teilintervallen von [a, b] ergibt sich unter Benutzung der Simpsonregel (6) auf jedem der Teilintervalle [x_2i, x_2i+1, x_2i+2] der Näherungswert

Ii  :=   h 3 ( f2i + 4 f2i+1 + f2i+2 )           i = 0, 1, ¼ , (N/2 - 1)   .

Durch Summation über alle Teilintervalle erhält man für das gesamte Intervall [a, b]

ó õ b a  f(x) dx    »     h 3 [ f0 + 4 f1 + 2 f2 + 4 f3 + ¼ + 2 fN-2 + 4 fN-1 + fN ]   =:   S(h)

(9)

Fehlerabschätzungen

Wir betrachten für die Trapezregel zunächst ein einzelnes Intervall [-[ h/2], [ h/2]]; dort sei der Integrand f(x) um x = 0 in eine Potenzreihe entwickelbar:

f(x) = a0 + a1 x + a2 x2 + a3 x3 + a4 x4 + ¼

Damit ist der Quadraturfehler (Verfahrensfehler, truncation error) der Trapezregel

R(h)  := T(h) - I(f)  =   h 2 é ë f ( -  h 2 ) + f (  h 2 ) ù û - ó õ h/2 -h/2  f(x) dx   .

Einsetzen der Potenzreihe in die rechte Seite der letzten Gleichung ergibt

R(h)  =   h 2 é ë 2a0 + 2a2  h2 4 + 2a4  h4 16 + ¼ ù û - é ë a0x +  1 2 a1x2 +  1 3 a2x3 +  1 4 a3x4 +  1 5 a4x5 + ¼ ù û +h/2 -h/2  = æ è a0h +  1 4 a2h3 +  1 16 a4h5 + ¼ ö ø - æ è a0h +  2 3 a2  h3 8 +  2 5 a4  h5 32 + ¼ ö ø =  1 4 a2h3 -  1 12 a2h3 + O(h5)   =    h3 12 (2a2) + O(h5)   =    h3 12 f¢¢(0) + O(h5)   .

In dieser Entwicklung des Verfahrensfehlers R(h) der Trapezregel treten nur ungerade Potenzen von h auf. Eine genauere Rechnung liefert

R(h) =  h3 12  f¢¢(x)

(10)

für eine Zwischenstelle x Î (-[ h/2], [ h/2]).

Für N Teilintervalle bei der zusammengesetzten Trapezregel ergibt sich daraus der Verfahrensfehler zu N ·R(h) , also wegen N = (b - a) / h

T(h) - I(f)  =   (b - a) h2 12  f¢¢(x)

(11)

für ein x Î (a, b). Der Fehler bei der zusammengesetzten Trapezregel geht also bei Verkleinerung der Teilintervalle von [a, b] mit h² gegen Null (Verfahren 2. Ordnung ).

Analog erhält man für die zusammengesetzte Simpsonregel den Quadraturfehler zu

S(h) - I(f)  =   (b - a) h4 180  f(4)(x)

(12)

für eine Zwischenstelle x Î (a, b). Die zusammengesetzte Simpsonregel ist also ein Verfahren 4. Ordnung und damit wesentlich genauer als die zusammengesetzte Trapezregel.

Anmerkung:
In der Praxis benützt man zur Kontrolle des Quadraturfehlers oft folgendes Kriterium für die Anzahl N der Teilintervalle von [a, b]: Man gibt eine Schranke e für die gewünschte relative Genauigkeit vor und berechnet durch fortgesetzte Halbierung der Schrittweite h bzw. Verdopplung der Anzahl der Teilintervalle N, d.h. h_n = (b - a) / 2ⁿ bzw. N = 2ⁿ für n = 0, 1, ¼, n_max , eine Folge I_n von Näherungswerten für I(f), bis entweder

| In - In-1 | < e |In-1|

oder n > n_max . Die zweite Bedingung erweist sich als nützlich, um bei einer zu restriktiven Vorgabe von e das Auftreten von Endlosschleifen zu vermeiden.

Programmierung der Simpsonregel mit Hilfe der Trapezregel

Wir betrachten für eine gerade Anzahl N von Teilintervallen von [a, b] folgende Linearkombination der Trapezsumme T(h) = h ( [ 1/2]f₀ + f₁ + ¼+ f_N-1 + [ 1/2]f_N ) zur Schrittweite h = (b - a) / N:

4 3  T(h) =  h 3 ( 2 f0 + 4 f1 + 4 f2 + 4 f3 + 4 f4 + ¼ + 2 fN )  1 3  T(2h) =  2h 3 æ è  1 2 f0 + f2 + f4 + ¼ +  1 2 fN ö ø =  h 3 ( f0 + 2 f2 + 2 f4 + ¼ + fN )  4 3  T(h) -  1 3  T(2h) =  h 3 ( f0 + 4 f1 + 2 f2 + 4 f3 + 2 f4 + ¼ + fN )   =  S(h)

Man erhält also durch eine Halbierung der Schrittweite (von 2h auf h) und Bildung einer Linearkombination aus den Trapezsummen T(2h) und T(h) die Simpsonregel in der Form

S(h)  =   4 3  T(h)  -   1 3 T(2h)

(13)

und damit eine Erhöhung der Fehlerordnung von h² auf h⁴.

Gauß-Quadratur

In diesem Abschnitt sollen Integrale des Typs

I(f)  :=  ó õ b a  w(x) f(x)  dx

(14)

mit einer für alle x Î [a, b] positiven Gewichtsfunktion w(x) > 0 durch Summen der Form

I(f)   »   n å i=1  wi f(xi)

(15)

approximiert werden. Die Stützstellen (Knoten) x_i und Gewichte w_i werden hier mit dem Index 1 beginnend durchnummeriert.

Bei den Newton-Cotes-Formeln sind n äquidistante Stützstellen x₁, ¼, x_n vorgegeben. Dazu werden n Integrationsgewichte so bestimmt, daß die resultierenden Quadraturformeln exakt sind für alle Polynome bis zum Grad £ (n - 1). Die Idee bei der Gauß-Quadratur besteht darin, die Einschränkung von äquidistant vorgegebenen Stützstellen aufzugeben und durch freie Wahl von 2 n Zahlen x_i und w_i als Stützstellen und Integrationsgewichte ( i = 1 ,¼, n ) Polynome möglichst hohen Grades exakt zu integrieren. Betrachtet man die Koeffizienten eines Polynoms als Parameter, so haben Polynome vom Grad (2 n - 1) auch 2 n Parameter. Man kann also erwarten, daß bei geeigneter Wahl der x_i und w_i Polynome vom Grad £ (2 n - 1) exakt integriert werden.

Als Beispiel betrachten wir den Fall n = 2 und das Integrationsintervall [-1, 1]. Es sollen Gewichte w₁, w₂ und Stützstellen x₁, x₂ so bestimmt werden, daß die Integrationsformel

ó õ +1 -1  f(x)  dx   »  w1 f(x1) + w2 f(x2)

(16)

das exakte Ergebnis liefert, wenn f(x) ein Polynom vom Grad 2 ·2 - 1 = 3 (oder weniger) darstellt, also insbesondere für f(x) = 1, x, x² und x³. Daraus gewinnt man die 4 benötigten Gleichungen für die 4 Unbekannten w₁, w₂, x₁, x₂:

w1 · 1  +  w2 · 1 = ó õ 1 -1  1  dx    =   2 w1 · x1  +  w2 · x2 = ó õ 1 -1  x  dx    =   0 w1 · x12  +  w2 · x22 = ó õ 1 -1  x2  dx =   2 / 3 w1 · x13  +  w2 · x23 = ó õ 1 -1  x3  dx =   0

Das ist ein System von 4 nichtlinearen Gleichungen. Aus (18) und (20) erhält man x₁² = x₂² und wegen x₁ ¹ x₂ also x₁ = - x₂. Damit liefert (18) w₁ = w₂, und wegen (17) ist dann w₁ = 1. Aus Gleichung (19) ergibt sich damit 2 x₁² = 2/3 oder x₁ = 1 / Ö3, also insgesamt

ó õ +1 -1  f(x)  dx   »   1 · f æ è -  Ö3 3 ö ø  +   1 · f æ è +  Ö3 3 ö ø

(17)

Im Prinzip könnte man mit diesem Verfahren auch die Stützstellen und Gewichte für n > 2 bestimmen, allerdings erhält man für die Stützstellen x_i wieder ein nichtlineares Gleichungssystem, dessen Lösung schwierig ist.

Man betrachtet daher im allgemeinen für auf [a, b] stetige Funktionen g und h das Skalarprodukt (g, h) bezüglich der Gewichtsfunktion w,

(g, h)  :=  ó õ b a  w(x) g(x) h(x)  dx   ,

(18)

und konstruiert einen Satz von orthogonalen Polynomen p₀(x), ¼, p_n(x), d.h. 

(pi, pk) = 0                für       i ¹ k   , (pi, pk) = ck > 0        für       i = k   ,

deren höchste Potenz von x jeweils den Koeffizienten 1 hat, d.h. die Polynome sollen normiert sein:

p0(x) = 1   , : pn(x) = 1 · xn + an-1 xn-1 + ¼ + a0   .

Diese Polynome heißen die zur Gewichtsfunktion w(x) und zum Intervall [a, b] gehörigen Orthogonalpolynome .

Man kann zeigen: Die Nullstellen { x₁,¼,x_n } von p_n(x) sind reell und einfach und liegen alle im offenen Intervall (a, b). Sie sind damit als Stützpunkte für eine Quadraturformel geeignet. Wählt man als Stützpunkte die Nullstellen { x₁,¼,x_n } von p_n(x) und als Integrationsgewichte { w₁,¼,w_n } die Lösung des linearen Gleichungssystems

æ ç ç ç ç ç è p0(x1) ¼ p0(xn) p1(x1) ¼ p1(xn) : : pn-1(x1) ¼ pn-1(xn) ö ÷ ÷ ÷ ÷ ÷ ø æ ç ç ç ç ç è w1 w2 : wn ö ÷ ÷ ÷ ÷ ÷ ø   =   æ ç ç ç ç ç è (p0,p0) 0 : 0 ö ÷ ÷ ÷ ÷ ÷ ø   ,

(19)

dann sind alle Gewichte w_i > 0 und die Quadraturformel (15) ist exakt für alle Polynome p(x) bis zum Grad £ (2 n - 1) , d.h. für diese Polynome gilt

ó õ b a  w(x) p(x) dx  =  n å i=1  wi  p(xi)   .

(20)

Die Stützstellen x_i und Integrationsgewichte w_i hängen vom Intervall [a, b] und von der Gewichtsfunktion w(x) ab. Für die praktisch wichtigste Gewichtsfunktion w(x) º 1 und das Intervall [-1, 1] stammen die Resultate von Gauß. Die zugehörigen Orthogonalpolynome p_k(x), k=0,1,¼,n sind bis auf einen Normierungsfaktor gerade die Legendre-Polynome P_k(x) und die so erhaltenen Quadraturformeln heißen Gauß-Legendre-Formeln :

ó õ 1 -1  f(x) dx   »   n å i=1  wi  f(xi)

(21)

Abbildung 4: Legendre-Polynome Pk(x).

Die P_k(x) sind so gewählt, daß P_k(1)=1. Die Nullstellen der Legendre-Polynome sind alle verschieden, liegen im offenen Intervall (-1, 1) und sind bezüglich des Ursprungs symmetrisch. Sie dienen als Stützstellen x_i in den Gauß-Legendre-Formeln. Man findet sie, ebenso wie die zugehörigen Integrationsgewichte w_i, in Tabellenwerken (bzw. mit Mathematica).

Die Festlegung des Integrationsintervalls auf [-1, 1] in (24) stellt keine Einschränkung dar, denn das Integral

I = ó õ b a  f(t) dt

läßt sich mit der Variablensubstitution

t =  b - a 2  x  +    a + b 2 dt =  b - a 2  dx

in ein Integral über das Intervall [-1, 1] transformieren:

I =  (b - a) 2 ó õ 1 -1  f æ è  b - a 2  x +  a + b 2 ö ø dx   .

Insgesamt erhalten damit die Gauß-Legendre-Formeln die Gestalt

ó õ b a  f(t) dt    »     (b - a) 2   n å i=1  wi  f æ è  b - a 2  xi  +    a + b 2 ö ø   .

(22)

Der Vorteil der Gaußschen Quadratur liegt darin, daß sie bei gleichem Rechenaufwand (gemessen an der Zahl der Funktionsauswertungen) die genauesten Resultate liefert, verglichen z.B. mit den Newton-Cotes-Formeln. Sie wird daher häufig bei Problemen mit aufwendigen Funktionsauswertungen sowie bei der Approximation von Mehrfachintegralen verwendet. (Bei letzteren wächst die Anzahl der Funktionsauswertungen mit einer Potenz der Anzahl der auszuwertenden Integrale.) Als Nachteil der Gaußschen Quadraturformeln erweist sich, daß die Stützstellen x_i und Gewichte w_i im allgemeinen nur numerisch bekannt sind; außerdem ändert sich bei einer Änderung von n auch die Lage der Stützstellen, sodaß die einmal für ein n berechneten Funktionswerte nicht für andere Werte von n wiederverwendet werden können.

Das folgende Beispiel soll den Genauigkeitsgewinn bei Verwendung der Gauß-Quadratur demonstrieren:

ó õ 1 -1  ex  dx  = e1 - e-1  = 2.3504   .

Mit der einfachen Trapezregel erhält man den Näherungswert

ó õ 1 -1  ex  dx   »    2 2 · ( e1 + e-1 )  = 3.0862   ,

wogegen eine Gauß-Legendre-Integration mit 2 Stützstellen x₁ = -0.57755, x₂ = 0.57755 und den Gewichten w₁ = 1, w₂ = 1 (vgl. (21)) ein wesentlich genaueres Resultat liefert:

ó õ 1 -1  ex  dx   »   1 · ( e-0.57755 + e0.57755 )  = 2.3429   .