Octave – Teil 3: Numerische Integration und Matrizen

Nachdem im zweiten Teil Funktionen und das Plotten von Graphen behandelt wurden, geht es nun an eher komplexe Themen wie die Programmierung und das Entwerfen von eigenen Funktionen. Auch werden Vektoren und Matrizen behandelt.

Von Jens Dörpinghaus

Hinweis: Alle Teile des Workshops findet man auf der Übersichtsseite.

Im vorherigen Abschnitt des letzten Artikels aus der Octave-Reihe wurde bereits eine Funktion f definiert. Aus der Schule kennt man die Integration, die man dort im Allgemeinen durch das Berechnen einer Stammfunktion löst. Anschaulich geht es um die Berechnung des Flächeninhaltes zwischen dem Graphen der Funktion f und der x-Achse. Es ist aber zum einen nicht immer möglich, eine Stammfunktion anzugeben und oft ist es einfach zweckmäßig eine numerische Näherung dieses Wertes zu bestimmen. Das soll an dieser Stelle die Motivation sein, ein erstes numerisches Problem zu lösen und Schleifen zu benutzen.

Numerische Integration mit der Trapezformel, for-Schleifen

Als Näherungsformel soll die Trapezformel angewandt werden.

Jens Dörpinghaus

Skizze zur Trapezformel

Dazu wählt man zwei x-Werte a und b, zwischen denen man den Flächeninhalt wie folgt annähert:

T(f) : = 1/2 ( f(a) + f(b) )&#183;(b& - a)

Diese Formel kann entweder über die Formeln von Newton-Cotes mit äquidistanten Stützstellen hergeleitet werden, oder es kann ganz anschaulich argumentiert werden. Der Flächeninhalt zwischen den Punkten a und b auf der x-Achse wird dabei durch eine Gerade approximiert, also angenähert. Der Flächeninhalt des Quadrates kann durch (b - a)f(a), der des Dreieckes durch 1/2 (b - a)(f(b) - f(a) angegeben werden. Addiert man beide Werte, so kommt man auf die obige Formel. Hier wurde stillschweigend angenommen, dass f(a) <= f(b) ist. Gilt dies nicht, kommt man aber trotzdem auf dieselbe Formel. Diese kann man direkt eingeben, wenn man die Funktionsdefinitionen aus dem zweiten Teil der Octave-Reihe im Pfad liegen hat:

>> a=2;
>> b=4;
>> 1/2* ( f(a) + f(b)) * (b-a)
ans =  20

Anschaulich kann man eine Grafik erstellen mit

>> x=a:0.1:b;
>> plot (x,f(x),'-',[a,b],[f(a),f(b)],'-r');

Die Funktion f und die Näherungsfunktion durch die Trapezformel (rot)

Jens Dörpinghaus

Die Funktion f und die Näherungsfunktion durch die Trapezformel (rot)

Zunächst wird also die Funktion f geplottet und dann die Näherungsfunktion in rot dargestellt. Es ist sinnvoll, die Trapezformel in eine separate Funktion auszugliedern:

function I = trapez(a,b,fa,fb)
  I = (1/2)* ( fa + fb) * (b-a);
end

Diese wird wieder unter einem Dateinamen, der dem Funktionsnamen entspricht, gespeichert. Ein kurzer Test mit

>> trapez(a,b,f(a),f(b))
ans =  20

liefert dasselbe Ergebnis. Man kann nun auf einem bestimmtes Intervall mit einer bestimmten Schrittweite die Ergebnisse der Trapezformel auf Teilintervallen zusammensetzen, um so ein genaueres Ergebnis zu erhalten. Die Idee ist, das Intervall [a,b] auf der x-Achse in gleiche Abstände s einzuteilen. Dann bekommt man Intervalle [a,a+s], [a+s, a+2s], ..., [b-s,b] auf die wieder die Trapezformel angewandt werden kann. Den Abstand s nennt man Schrittweite. Eine Funktion mit kleinerer Schrittweite kann »genauer« approximiert werden kann. Dazu nutzt man folgende Funktion:

function I = zusammengesetztTrapez(a,b,schrittweite)
  % Liste von Stuetzpunkten
  x = a:schrittweite:(b-schrittweite);
  integral = 0;
  % Trapezformel anwenden
  for stuetzpunkt = x,
    stuetzpunkt2 = stuetzpunkt + schrittweite;  
    integral = integral + trapez(stuetzpunkt, stuetzpunkt2, f(stuetzpunkt), f(stuetzpunkt2));
  end;
  I = integral;
end