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 xi = 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 x0 = a und xn = b in n Teilintervalle der Breite (Schrittweite) h : = (b - a) / n unterteilt. Sei P01n(x) das zugehörige Interpolationspolynom, d.h.
P01n(xi) = f(xi) = fi        (i = 0, 1, , n)   .
(4)


n = 1 : Lineares Interpolationspolynom P01 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 (x0,0), (x1,0), (x1,f1), (x0,f0) ist.

trapez.png

Abbildung 1: Trapezregel.

simpson.png

Abbildung 2: Simpsonregel.


n = 2 : Quadratisches Interpolationspolynom P012 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 P012(x):
P012(x) = a (x - x0)2 + b (x - x0) + c   .
(7)

Mit diesem Ansatz und den Bedingungen P012(xi) = fi 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 (x1 - x0) = h bzw. (x2 - x0) = 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

ztrapez.png

Abbildung 3: Zusammengesetzte Trapezregel.

Für die Unterteilung xi = x0 + ih, (i=0,1,,N), von [a,b] mit x0=a, xN=b und h=(b-a)/N erhält man aus der Trapezregel (5) im Teilintervall [xi,xi+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 [x2i, x2i+1, x2i+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 h2 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. hn = (b - a) / 2n bzw. N = 2n für n = 0, 1, , nmax , eine Folge In von Näherungswerten für I(f), bis entweder
| In - In-1 | < e |In-1|

oder n > nmax . 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]f0 + f1 + + fN-1 + [ 1/2]fN ) 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 h2 auf h4.




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) xi und Gewichte wi werden hier mit dem Index 1 beginnend durchnummeriert.

Bei den Newton-Cotes-Formeln sind n äquidistante Stützstellen x1, , xn 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 xi und wi 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 xi und wi 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 w1, w2 und Stützstellen x1, x2 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, x2 und x3. Daraus gewinnt man die 4 benötigten Gleichungen für die 4 Unbekannten w1, w2, x1, x2:
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 x12 = x22 und wegen x1 x2 also x1 = - x2. Damit liefert (18) w1 = w2, und wegen (17) ist dann w1 = 1. Aus Gleichung (19) ergibt sich damit 2 x12 = 2/3 oder x1 = 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 xi 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 p0(x), , pn(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 { x1,,xn } von pn(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 { x1,,xn } von pn(x) und als Integrationsgewichte { w1,,wn } 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 wi > 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 xi und Integrationsgewichte wi 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 pk(x), k=0,1,,n sind bis auf einen Normierungsfaktor gerade die Legendre-Polynome Pk(x) und die so erhaltenen Quadraturformeln heißen Gauß-Legendre-Formeln :

1

-1 
f(x) dx     n

i=1 
wi  f(xi)
(21)



Tabelle 1: Orthogonalpolynome pk(x) und Legendre-Polynome Pk(x).

p0(x) = 1
      P0(x) = 1
p1(x) = x
      P1(x) = x
p2(x) = x2 -  1

3
      P2(x) =  3

2
·p2(x) =  1

2
(3 x2 - 1)
p3(x) = x3 -  3

5
x
      P3(x) =  5

2
·p3(x) =  1

2
(5 x3 - 3 x)
      Pn+1(x) =  2 n + 1

n + 1
x Pn(x) -  n

n + 1
Pn-1(x)



legendre.png
Abbildung 4: Legendre-Polynome Pk(x).

Die Pk(x) sind so gewählt, daß Pk(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 xi in den Gauß-Legendre-Formeln. Man findet sie, ebenso wie die zugehörigen Integrationsgewichte wi, 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 xi und Gewichte wi 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 x1 = -0.57755, x2 = 0.57755 und den Gewichten w1 = 1, w2 = 1 (vgl. (21)) ein wesentlich genaueres Resultat liefert:

1

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



File translated from TEX by TTH, version 3.06.
On 1 Apr 2004, 18:24.