Numerische Mathematik, Teil 3: Integration und Differentiation

In dieser Folge befassen wir uns mit der Differential- und Integralrechnung. Dazu werden einige Begriffe erklärt und Beispielprogramme zu grundlegenden Algorithmen aufgeführt.

Was eine Funktion ist, wissen wir bereits aus der ersten Folge. Funktionen haben spezielle Eigenschaften. Stetige Funktionen haben einen Graphen, den man quasi mit einem Linienzug zeichnen kann, d.h. wenn man ihn zeichnen wollte, bräuchte man den Stift nicht abzusetzen. Ist eine Funktion stetig differenzierbar, verläuft der Graph so schön glatt, daß keine “Knicke” auftauchen. Wird eine Funktion f differenziert, erhält man die erste Ableitung dieser Funktion, genannt f. Auch die Ableitung f kann über die Eigenschaften Stetigkeit und Differenzierbarkeit verfügen. Ist f differenzierbar, kommt man zur zweiten Ableitung f'. Sofern also Ableitungen nullten, ersten, zweiten,..., n-ten Grads differenzierbar sind, können wieder Ableitungen einer höheren Stufe gebildet werden.

Es gibt Funktionen, die einmal, aber nicht zweimal stetig differenzierbar sind, auch welche, die nicht einmal stetig sind. Nehmen wir als Beispiel die Betragsfunktion f(x)=/x/. Die Funktion ist bei x=0 nicht differenzierbar, da der zugehörige Graph dort einen “Knick” hat.

Ist eine Funktion f stetig, hat sie eine Stammfunktion F, d. h. die Ableitung von F ist wieder f Auch wenn die Existenz einer Stammfunktion bekannt ist, kann es unter Umständen sehr schwierig sein, diese zu bestimmen. Aber mit einer Stammfunktion F von f kann man eine Funktion f integrieren, d. h. den Wert des Integrals I von f über einem Intervall [a;b] angeben. Zudem gibt es auch noch nichtstetige Funktionen, die trotzdem integrierbar sind. Da für viele Probleme nur ein Näherungswert eines Integrals benötigt wird, greifen in diesen Fällen unsere Methoden der Numerischen Mathematik!

Ableitungen und Integrale haben neben theoretischen auch anschauliche Bedeutungen. Ein Integral von /über einem Intervall gibt die Fläche an, die der Graph von/in diesem Intervall und die (x-)Koordinatenachse umspannen. Dabei werden Flächen unterhalb der Koordinatenachse negativ, die Flächen oberhalb positiv gezählt. Die erste Ableitung einer Funktion in einem Punkt x gibt die Steigung des Graphen der Funktion in diesem Punkt an, die zweite sagt etwas über das Krümmungsverhalten aus.

Der Integralbegriff, der hier zugrunde liegt, stammt von Bernhard Riemann (geb. 17.9.1826 in Breselenz, gest. 20.7. 1866 in Salesca). Er wird in der Regel über Treppenfunktionen eingeführt. Hierzu werden die Koordinatenachse in kleine Teilintervalle unterteilt und dann die Flächen von passenden Rechtecken zwischen Graph und x-Achse berechnet (siehe Bild 2). Näherungswerte für den korrekten Flächeninhalt erhält man damit durch geeignete Verkleinerung der Intervallängen zwischen den Intervallgrenzen. Aber wir wollen die Theorie ja nicht zu weit treiben, kommen wir zur Numerik.

In der letzten Folge haben wir Interpolationsformeln kennengelernt. Darunter war auch die Lagrangesche. Integriert man diese explizit, so erhält man einen Satz von sehr bekannten Näherungsformeln für den Integralwert. Die explizite Integration bereitet keine größeren Schwierigkeiten, da die Lagrangesche Interpolationsformel ein Polynom ist. Die entstehenden Formeln heißen Newton-Cotes-Formeln.

Die einzelnen Formeln, die sich nach der Anzahl der Stützstellen des Interpolationsproblems richten, haben wiederum Namen. So entsteht aus nur einer Stützstelle die bekannte Trapezregel, aus zwei Stützstellen die Simpson-Regel (nach Thomas Simpson, geb. 20. 8. 1710 in Market-Bosworth und gest. 14. 5. 1761 ebenda), aus drei Stützstellen die 3/8-Regel, aus vier Stützstellen die Milne-Regel und aus sechs Stützstellen die Weddle-Regel.

Bei den Formeln von Newton-Cotes gelten die Intervallgrenzen immer als Stützstellen. Bei einer anderen Gruppe von Integrationsformeln wie z. B. den Formeln von Euler und Maclaurin ist dies nicht der Fall.

Als Beispiel zu den Newton-Cotes-Formeln wurde für das Listing 1 die Simpson-Regel benutzt. Allerdings wurde das Intervall weiter unterteilt, um eine höhere Genauigkeit zu erreichen. Damit haben wir die sogenannte zusammengesetzte Simpson-Regel.

Berechnet man nicht wie bei den Newton-Cotes-Formeln die Stützstellen und damit die Gewichte, sondern wählt Stützstellen und Gewichte so, daß eine maximale Genauigkeit erreicht wird, kommt man zu den Integrationsformeln von Gauß.

Als Stützstellen ergeben sich dann gerade die Nullstellen eines passenden Legendre-Polynoms. Auf diese Berechnung wollen wir hier jedoch nicht näher eingehen.

Die Integrationsformeln von Gauß sind also genauer als die Formeln von Newton-Cotes. Sie haben aber einen Nachteil, der gerade noch nicht klar wurde: Die Berechnung der Stützstellen ist aufwendiger. Außerdem erhält man nicht ganz so “schöne” Werte für die Koeffizienten, wie es z.B. bei den Newton-Cotes-Formeln der Fall ist.

Das zweite Beispielprogramm zeigt die Implementation dreier Gaußscher Integrationsformeln. An den Ergebnissen erkennt man jedoch, daß auch in unserem Beispiel die Formeln von Gauß genauer sind als die von Newton-Cotes.

Geht man den anderen Weg und differenziert das Interpolationspolynom von Lagrange, erhält man Formeln für die numerische Differentiation. Das Interpolationspolynom kann man theoretisch beliebig oft ableiten und erhält dann auch Formeln für höhere Ableitungen, die jedoch mit äußerster Vorsicht zu genießen sind, da sie meistens nicht genauer werden.

Man überlege sich dazu nur, daß man die Funktion sin(x) mit einem Interpolationspolynom lokal nähert. Beliebig hohe Ableitungen werden nie konstant gleich 0, aber wenn man das Interpolationspolynom ausreichend oft differenziert, wird es irgendwann konstant gleich 0.

Trotzdem erhält man zumindest für die erste und zweite Ableitung interessante Formeln. Sie heißen je nach Anzahl ihrer Stützstellen 3-, 4-, 5-Punkte-Formel,... . (Die Formeln finden sich in Listing 3.)

Das dritte Beispielprogramm wertet drei Formeln dieser Art aus. Im Vergleich mit dem im folgenden beschriebenen Romberg-Verfahren erhält man jedoch meistens schlechtere Werte. Das Romberg-Verfahren hat einen Vorteil, der hier vorweggenommen werden soll: Die Genauigkeit ist nicht nur höher, sie kann vorher bestimmt werden.

Ähnlich den Interpolationsverfahren von Newton und Neville aus der letzten Folge unserer Serie kann man auch beim Romberg-Verfahren ein Rechenschema angeben. Es ist in Bild 7 zu sehen. Der Index, der oben rechts am D steht, bezeichnet den Schritt, der Index unten rechts, der wievielte Wert des passenden Schritts vorliegt. Ausgehend von der Anzahl der Werte pro Schritt wird die Schrittweite hj gewählt. Sie wird für jeden folgenden Wert eines Schritts halbiert.

Das Rechenverfahren wird solange fortgesetzt, bis die gewünschte Genauigkeit erreicht ist. Das ist genau dann der Fall, wenn - wie in Listing 4 zu sehen ist -abs(d[1]-d[0]) unterhalb einer vorgegebenen Obergrenze liegt. Allerdings sollte man eine Maximalzahl an Schritten vorgeben, denn auch beim Romberg-Verfahren kann es passieren, daß Konvergenz nicht erreichbar ist.

Der einzige Nachteil des Verfahrens: Die Programmierung, der Speicherbedarf und die Rechenzeiten sind größer als bei den einfachen n-Punkte-Formeln.

Nun sind wir schon am Ende der dritten Folge angelangt. In der nächsten wenden wir uns einem speziellen Integrationsproblem zu: der Berechnung einer Lösung einer gewöhnlichen Differentialgleichung erster Ordnung.

Auf bald...

Dipl.-Math. Dietmar Rabich

Literatur:

[1] Einführung in die Numerische Mathematik I, J. Stoer, Springer Berlin/Heidelberg/New York/Tokyo, 4. Aufl. 1983, S. 106ff

[2] Formelsammlung zur Numerischen Mathematik mit BASIC-Programmen, G. Engeln-Müllges/F. Reutter, Bibliographisches Institut Mannheim/Wien/Zürich, 1. Aufl. 1983, S. 186ff

[3] Numerische Mathematik, H. R. Schwarz, Teubner Stuttgart, 1. Aufl. 1986, S. 319ff

[4] Numerische Methoden, A. Björck/G. Dahlquist, Oldenbourg München! Wien, 2. Aufl. 1979, S. 213ff

[5] Methode der Numerischen Mathematik, W. Böhm/G. Gose/ J. Kahmann, Vieweg Braunschweig/Wiesbaden, 1. Aufl. 1985, S. 122ff

[6] Erfolgreich programmieren mit C, J. A. Illik, Sybex Düsseldorf/San Francisco/Paris/London, 4. Aufl. 1987

[7] Programmieren in C, B. W. Kernighan/D. M. Ritchie, Hanser München/Wien, 1. Aufl. 1983

[8] PASCAL für Anfänger. H. Schauer, Oldenbourg Wien/München, 4. Aufl. 1982

[10] PASCAL für Fortgeschrittene. H. Schauer, Oldenbourg Wien/München, 2. Aufl. 1983

/************************************************/ 
/* Beispielprogramm zur Integration mit der     */
/* summierten Simpsonregel                      */
/*----------------------------------------------*/
/* Entwickelt mit Turbo C.  15.02.1989          */
/************************************************/

/*---------- + ------------------*/
/* Listing 1 / by D. Rabich      */
/*   (c) MAXON Computer GmbH     */

/* Ein-/Ausgabe-Routinen importieren */
# include <stdio.h>

/* Grenze zur Berücksichtigung eventueller Ungenauigkeiten */
# define FLOATFEHLER 1.0E-6

/* Funktion f */ 
float f(float *x)
{
    return((*x)*(*x));
}

/* Berechnet Integral der Funktion f zwischen a und b */
float simpson(float (*f)(), float a, float b, int n)
{
    float x,
          integ,
          h;

    integ=0;

    /* Schrittweite */ 
    h=(b-a)/n;

    /* Integral berechnen */ 
    for(x=a+h/2;x<b-h/2+FLOATFEHLER;x+=h) 
        integ+=(*f)(&x);

    integ*=2;

    for (x=a;x<=b+FLOATFEHLER;x+=h) 
        integ+=(*f)(&x);

    return(integ*h/3);
}

/* Hauptprogramm */ 
void main(void)
{
    int n;

    /* Ausgabe bei verschiedenen Unterteilungen */ 
    for(n=5;n<=50;n+=5) 
        printf("Integral bei %2d Unterteilungen: %12.10g\n", n,simpson(f,0,1,n));

    /* Auf Taste warten... */
    getchar();
}

(************************************************) 
(* Beispielprogramm zur Gauß-Integration.       *)
(* -------------------------------------------- *)
(* Entwickelt mit ST Pascal Plus. 14.02.1989    *)
(************************************************)

(*---------- + ------------------ *)
(* Listing 2 / by D. Rabich       *)
(*  (c) MAXON Computer GmbH       *)

program gauss_intgration;

(* Typ *)
type poss_n = 0..2;

(* Variablen *) 
var i : poss_n;

(* Funktion f *)
function f (x : real) : real;

begin 
    f:=x*x 
end;

(* Berechnet Integral von Funktion f zwischen a und b *)
function gauss(function f (x : real) : real;
                                 a,b : real; 
                                 n   : poss_n)
                                     : real;

var fakt,h ; real;

begin 
    h:=(b-a)/2; 
    case n of

        (* Formel für eine Stützstelle *)
        0 : gauss:=2*h*f(a+h);

        (* Formel für zwei Stützstellen *)
        1 : begin
                fakt :=1.0/sqrt(3.0); 
                gauss:=h*(f(a+h*(1.0-fakt))+f(a+h*(1.0+fakt)))
            end;

        (* Formel für drei Stützstellen *)
        2 : begin
                fakt :=sqrt(0.6);
                gauss:=h*(5.0*(f(a+h*(1.0-fakt))+f(a+h*(1.0+fakt)))+8.0*f(a+h))/9.0
            end

        end
    end;

(* Hauptprogramm *) 
begin

    for i:=0 to 2 do 
        writeln(gauss(f,0.0,1.0,i));

    repeat
    until keypress 

end.

(************************************************************) 
(* Beispielprogramm zur numerischen Differentiation mittels *)
(* Differentiation des Interpolationspolynoms.              *)
(* -------------------------------------------------------- *)
(* Entwickelt mit ST Pascal Plus.                14.02.1989 *) 
(************************************************************) 

(* --------- + ----------------- *)
(* Listing 3 / by D. Rabich      *)
(* (c) MAXON Computer GmbH       *)

program intpol_differentiation;

(* Typ *)
type poss_n = (drei,fuenf,sieben);

(* Variable *) 
var i : poss_n;

(* Funktion f *)
function f (x : real) : real;

begin
    f:=x*x+x*x*x 
end;

(* Berechnet Ableitung der Funktion f an der Stelle x *) 
function diff_intpol(function f (x real) : real;
                                     x,h : real; 
                                       n : poss_n): real;

begin 

    case n of 

        (* 3 Stützstellen *)
        drei   : diff_intpol:=(- f(x- h)+ f(x+ h)) / (2*h);

        (* 5 Stuetzstellen *)
        fuenf  : diff_intpol:=(  f(x-2*h)- 8*f(x- h)+ 8*f(x+ h) - f(x+2*h)) /(12*h);

        (* 7 Stützstellen *)
        sieben : diff_intpol:=(- f(x-3*h)+ 9*f(x-2*h)-45*f(x- h)+45*f(x+ h)- 9*f(x+2*h)+ f(x+3*h)) /(60*h)
    end

end;

(* Hauptprogramm *) 
begin

    for i:=drei to sieben do 
        writeln(diff_intpol(f,0.5,1,i));

    repeat
    until keypress 

end.

(*************************************************************) 
(* Beispielprogramm zur numerischen Differentiation nach dem *)
(* Verfahren von Romberg.                                    *)
(* --------------------------------------------------------- *)
(* Entwickelt mit ST Pascal Plus.                 14.02.1989 *) 
(*************************************************************) 

(* --------- + -------------- *)
(* Listing 4 / by D. Rabich   *)
(*  (c) MAXON Computer GmbH   *)

program romberg_differentiation;

(* Konstante *) 
const real_fehler = 1.0E-8; 
      max_ber     = 50;

(* Funktion f *)
function f (x : real) : real;

begin 
    f:=x*x+x*x*x 
end;

(* Berechnet Ableitung der Funktion f an der Stelle x *) 
function romberg (function f (x : real) : real;
                                    x,h : real): real;

var d     : array [0..max_ber] of real;
    i,j,p : short_integer;

begin

    (* Romberg-Verfahren *) 
    d[0]:=(f(x+h)-f(x-h))/(2*h); 
    i:=1; 
    repeat
        h:=h/2;
        d[i]:=(f(x+h)-f(x-h))/(2*h); 
        p:=4;
        for j:=i-1 downto 0 do 
        begin
            d[j]:=(p*d[j+1]-d[j])/(p-1); 
            p:=p*4 
        end; 
        i:=i+1
    until (abs(d[1]-d[0])<real_fehler) or (i>max_ber);

    (* Fehlermeldung, falls Genauigkeit nicht erreicht *) 
    if abs(d[1]-d[0])>=real_fehler then 
        writeln('*** Genauigkeit nicht erreicht ***');

    romberg:=d[0] 
end;

(* Hauptprogramm *) 
begin

    writeln(romberg(f,0.5,1)); 

    repeat
    until keypress 
end.

Dietmar Rabich

---

---