Ing Mathematik: Numerische Berechnung bestimmter Integrale

Zurück zu Approximation und Interpolation | Hoch zu Gesamtinhaltsverzeichnis | Vor zu Näherungsweises Lösen gewöhnlicher Dgl.

Gesucht ist $I=\int _{a}^{b}f(x)dx$. Mittels Quadraturverfahren kann man dieses Integral numerisch lösen.

Rechteckformel

Die Fläche unter der Kurve $f$ zwischen $a,b$ ist näherungsweise das rot eingefärbte Rechteck. Als Formel geschrieben:

$I\approx f\left({\frac {a+b}{2}}\right)\cdot (b-a)$.

Das ist nicht sehr genau. Besser ist es, das Rechteck in mehrere kleinere Rechtecke aufzusplitten (= summierte Rechteckregel).

$h={\frac {b-a}{n}},\quad x_{i}=a+i\cdot h$

$I\approx h\sum _{i=0}^{n-1}f\left(x_{i}+{\frac {h}{2}}\right)$.

Python-Programm für eine Parabel:

import numpy as np

def funktion(x):
    return x**2 

def rechteckregel(a, b, n): 
    I = 0.
    h = (b-a) / n

    for i in np.arange(0, n):
        xi = a + i*h
        I += h * funktion(xi + h/2)

    return I

A1 = rechteckregel(-1, 1, 10)
A2 = rechteckregel(-1, 1, 100)
A3 = rechteckregel(-1, 1, 1000)

print("A1 = ", A1)
print("A2 = ", A2)
print("A3 = ", A3)
   

Ausgabe:

A1 =  0.6600000000000001
A2 =  0.6665999999999999
A3 =  0.6666659999999996

Man sieht, dass die Folge nur sehr langsam konvergiert. Es sind für die praktsche Anwendung bessere Verfahren notwendig.

Python-Programm für die eingefärbte Fläche unter der Parabel:

import numpy as np
import matplotlib.pyplot as plt

x = np.arange(-1., 1.01, 0.01)

y = x**2

plt.plot(x, y, color='black')
plt.fill_between(x, y, where=y>=0)

plt.show()

Ausgabe:

Siehe auch Knorrenschild: Seite 128f

Trapezformel

Wiederum wird die rot eingefärbte Fläche berechnet.

$I\approx \left({\frac {f(a)+f(b)}{2}}\right)\cdot (b-a)$

Es sei wiederum $h={\frac {b-a}{n}},\quad x_{i}=a+i\cdot h$ mit $n\in \mathbb {N}$

$I\approx h\cdot \left({\frac {f(a)+f(b)}{2}}+\sum _{i=1}^{n-1}f(x_{i})\right)$.

Python-Programm für eine Parabel:

import numpy as np

def funktion(x):
    return x**2

def trapezregel(a, b, n): 
    h = (b-a) / n
    I = 0.

    for i in np.arange(1, n):
        xi = a + i*h
        I +=  h * funktion(xi)
        
    I += h * (funktion(a)+ funktion(b))/2.
    return I

A1 = trapezregel(-1, 1, 1)
A2 = trapezregel(-1, 1, 100)
A3 = trapezregel(-1, 1, 1000)
print("A1 = ", A1)
print("A2 = ", A2)
print("A3 = ", A3)

Ausgabe:

A1 =  2.0
A2 =  0.6668000000000002
A3 =  0.6666679999999993

Siehe auch Knorrenschild: Seite 127ff, Hanke-Bourgeois: Seite 317ff und Burg, Haf, Wille, Meister: Seite 331f

Simpsonformel

Bei der Verwendung der Simpson-Regel muss die Funktion im gegebenen Intervall vierfach stetig differenzierbar sein. Die Funktion wird durch ein Polynom zweiter Ordnung angenähert.

Es ergibt sich folgende Formel:

$I\approx {\frac {h}{3}}\cdot \left(f(x_{0})+4f(x_{1})+2f(x_{2})+4f(x_{3})+2f(x_{4})+\dots +4f(x_{n-1})+f(y_{n}))\right)$.

mit $h,x_{i}$ wie oben (siehe z.B. Rechteck- oder Trapez-Verfahren). $n$ sei gerade.

Python-Programm:

import numpy as np

def funktion(x):
    return x**2

def simpson(a, b, n): 
    if n%2 == 0:
        h = (b-a) / n
        I = 0.

        for i in np.arange(0, n+1):
            xi = a + i*h
            if i % 2 == 0:
                if i == 0 or i == n:
                    I += h / 3. * funktion(xi)
                else:
                    I +=  h / 3. * 2 * funktion(xi)
            else:
               I +=  h / 3. * 4 * funktion(xi)
        return I
    else:
        print ("Abbruch: n nicht gerade!")
        return


A1 = simpson(-1, 1, 2)
A2 = simpson(-1, 1, 10)
A3 = simpson(-1, 1, 100)
print("A1 = ", A1)
print("A2 = ", A2)
print("A3 = ", A3)

Ausgabe:

A1 =  0.6666666666666666
A2 =  0.6666666666666667
A3 =  0.6666666666666667

Siehe auch  Simpsonregel, Burg, Haf, Wille, Meister: Seite 332f, Knorrenschild: Seite 130 und Hanke-Bourgeois: Seite 326ff.

Romberg-Verfahren

Werner Romberg (deutscher Mathematiker, 1909-2003)

Das Romberg-Verfahren wurde von W. Romberg um 1955 entwickelt. Es ist eine Verbesserung der Trapezregel durch Extrapolation. Das SciPy-Paket beinhaltet eine romberg-Funktion. Diese ist allerdings als deprecated markiert und soll nicht mehr verwendet werden. Stattdessen gibt es die Methode quad, die nachfolgend auch kurz beschrieben wird.

Siehe  Romberg-Integration, Knorrenschild: Seite 140ff, Hanke-Bourgeois: Seite 332ff.

Allgemeine Quadratur mit Python und SciPy

Python und insbesondere das Python-Paket SciPy bieten natürlich wieder eine vorgefertigte Möglichkeit der numerischen Integration. Und zwar die Funktion quad. Auch Mehrfachintegrale lassen sich lösen. Dafür gibt es die z.B. die Funktion dblquad. Dies sei hier aber nicht vertieft. Siehe ggf. die SciPy-Dokumentation [1].

Python:

from scipy import integrate

y = lambda x: x**2

I, err = integrate.quad(y, -1, 1)
print(I)

Ausgabe:

0.6666666666666666

Gedruckte Werke (auszugsweise)

• Burg, Haf, Wille, Meister: Höhere Mathematik für Ingenieure, Band I: Analysis. 9. Auflage, Vieweg+Teubner, 2011, ISBN 978-3-8348-1218-6
• Hanke-Bourgeois: Grundlagen der Numerischen Mathematik und des Wissenschaftlichen Rechnens. 3. Aufl., Vieweg+Teubner, 2009, ISBN 978-3-8348-0708-3
• Knorrenschild: Numerische Mathematik, Eine beispielorientierte Einführung. 6. Aufl., Hanser, 2017, ISBN 978-3-446-45161-2

Zurück zu Approximation und Interpolation | Hoch zu Gesamtinhaltsverzeichnis | Vor zu Näherungsweises Lösen gewöhnlicher Dgl.