icancode.de

Python

Approximation einer Differentialgleichung plotten

Veröffentlicht am .

Approximation einer Differentialgleichung plotten

Einleitung

Wer sich schon immer mal gefragt hat, wie man einfach mit Python plotten kann, dem wird heute hier geholfen ;)

Im Folgenden zeige ich, wie man die Differentialgleichung

u'(t) = sin(u(t)) + cos(u(t) + t)

numerisch mithilfe des expliziten Euler-Verfahrens approximiert, also eine Näherungslösung dazu erstellt, und diese Näherungen dann anschließend plottet.

Die mathematischen Winkelfunktionen sin()  und cos()  liefert uns die Python-Bibliothek math  oder auch numpy , wobei Letzterem wegen seiner besseren Optimierung klar der Vorzug gegeben werden sollte. Außerdem benötigen wir noch die Funktion arange()  aus dem Paket numpy, mit welchem wir uns einen zusätzlichen Vektor (in Python arbeiten wir auf Listen!) eine zusätzliche Liste mit den Werten für die x-Koordinaten bilden. Aber dazu gleich mehr.

Die numerische Approximation mit dem expliziten Euler-Verfahren bekommt als Eingabe, neben der Funktion, eine Schrittweite h. Diese Schrittweite ist für die Genauigkeit der Approximation wichtig, daher testen wir in unserem Beispiel die Schrittweiten h = [2, 1, 0.5, 0.25, 0.125].

Damit wir auf die Plotfunktion in Python zurückgreifen können, muss die Bibliothek pylab  geladen werden: Sie beinhaltet matplotlib , wie der Name schon andeutet, eine Bibliothek zum Plotten von mathematischen Funktionen und Werten.

Nachfolgend ersteinmal der Code, inkl. Erläuterungen:


# Importieren der nötigen Bibliotheken:
import numpy as m, pylab as p

# Definition einer neuen Funktion, in Abhängigkeit von u und t:
def Function(t, u):
  return m.sin(u) + m.cos(t + u)

Unser Verfahren benötigt als Eingabe die Funktion, einen Anfangswert, ein geschlossenes Intervall und die Schrittweite:

def ExplEuler(func, startingValue, intervalStart, 
              intervalEnd, stepSize):

Die Werte des Anfangswertes und der erste Wert des Intervalls werden nun in neue Variablen kopiert, damit man diese auch inkrementieren kann. Hier erfolgt die Instanziierung und Initialisierung in einem Schritt für beide Variablen gleichzeitig:

u, t = startingValue, intervalStart

Außerdem benötigen wir noch eine leere Liste, in der wir die Werte der jew. Approximation zwischenspeichern können:

uList = []

Damit wir alle Werte des geschlossenen Intervalls erreichen, muss in der nun folgenden Schleife <=  (größer-gleich) verwendet werden, denn das Ende des Intervalls wird ja explizit mit eingeschlossen:

while t <= intervalEnd:

Als Erstes wird nun der aktuelle Wert von u in die Liste der u-Werte kopiert, beim ersten Schleifendurchlauf ist dies noch der Anfangswert der Differentialgleichung:

uList.append(u)

Erst wird also der neue Wert für u berechnet

u += stepSize * func(t, u)

und anschließend t mit der Schrittweite inkrementiert:

t += stepSize

Nun folgt der interessante Teil: Außerhalb der Schleife packen wir nun alle u-Werte in der Liste zusammen mit einer weiteren Liste, die allen t-Werten entspricht, in die plot -Funktion. Wir lassen die t-Liste von der Funktion arange()  erstellen, müssen jedoch beachten, dass die beiden Listen (Vektoren) die gleiche Länge haben: Deshalb muss auf das Intervallende noch eine zusätzliche Schrittweite addiert werden. Komplettiert wird dies noch durch die Angabe des Labels, welches wir hier auf die Schrittweite h setzen:

p.plot(m.arange(intervalStart, intervalEnd + stepSize,
stepSize), uList, label=stepSize)

Die im Text angegebenen Schrittweiten werden nun noch in einer Liste an die obige Funktion übergeben, beachtet hier auch den Anfangswert und das geschlossene Intervall:

for stepSize in [2, 1, 0.5, 0.25, 0.125]: 
ExplEuler(Function, 0, 0, 50, stepSize)

Jetzt fehlen nur noch eine Legende und die Axenbeschriftungen, welche man wie folgt einfügen kann:

p.legend(bbox_to_anchor=(0.775, 0.8), loc=2,
borderaxespad=0., title="h") 
p.xlabel("t"), p.ylabel("u(t)")

Ist dies erledigt, kann man sich den Plot nun anzeigen, alternativ auch in ein beliebiges Grafikformat abspeichern lassen (siehe auskommentierte letzte Zeile):

p.show() 
#p.savefig('A30.svg') # or plot into a svg file

Kompletter Code

Damit der Code funktioniert, muss numpy  und die matplotlib  auf dem System vorhanden sein.

import numpy as m, pylab as p

def Function(t, u):
  return m.sin(u) + m.cos(t + u)

def ExplEuler(func, startingValue, intervalStart, 
              intervalEnd, stepSize):
  u, t = startingValue, intervalStart  
  uList = [] 
  while t <= intervalEnd: 
    uList.append(u)
    u += stepSize * func(t, u)
    t += stepSize
  p.plot(m.arange(intervalStart, intervalEnd + stepSize,
         stepSize), uList, label=stepSize)

for stepSize in [2, 1, 0.5, 0.25, 0.125]: 
  ExplEuler(Function, 0, 0, 50, stepSize)

p.legend(bbox_to_anchor=(0.775, 0.8), loc=2,
         borderaxespad=0., title="h") 
p.xlabel("t"), p.ylabel("u(t)")
p.show()                                  
#p.savefig('A30.svg')                    

num.approx. DGL

Martin Wiechmann

Martin Wiechmann

Martin studiert an der Uni Bielefeld kognitive Informatik. Er ist hier der Linux-Crack und verfügt in diesem Bereich über langjährige Erfahrungen. Hobby = Hund :-)

Navigation