12.10    Lösen von Differenzialgleichungen mit Sympy

Mit dem Modul SymPy lassen sich Differenzialgleichungen auch symbolisch lösen. Allerdings stößt SymPy schnell an seine Grenzen. Die nichtlineare DGL eines mathematischen Pendels kann SymPy nicht lösen. Deshalb wird das Federpendel aus Abbildung 12.12 gewählt.

Die DGL 2. Ordnung

können Sie direkt ohne Umformung

diff(x(t),t,2) + d/m*diff(x(t),t,1) + c/m*x(t)

in einen SymPy-Quelltext übertragen.

Listing 12.13 berechnet für die DGL des Federpendels die allgemeine und die spezielle Lösung. Gelöst wird eine DGL mit der SymPy-Methode dsolve(dgl). Die Anfangswerte x(0) und x'(0) werden in den Zeilen 10 und 11 vorgegeben. Der Verlauf der Pendelschwingungen wird in Zeile 24 als Funktionsplot dargestellt.

01  #13_sympy_dgl.py
02  from sympy import * 
03  t = symbols('t')
04  x = Function('x')
05  m=1.0   #Masse in kg
06  d=0.2   #Dämpfung 
07  c=9.87  #Federkonstante in N/m
08  #Anfangswerte
09  aw={
10      x(0):5,                     #Auslenkung
11      x(t).diff(t,1).subs(t,0):0  #Anfangsgeschwindigkeit
12      }
13  #Lösung der DGL
14  #dgl=Eq(diff(x(t),t,2) + d/m*diff(x(t),t) + c/m*x(t),0)
15  #dgl=Eq(x(t).diff(t,2) + d/m*x(t).diff(t) + c/m*x(t),0)
16  dgl=diff(x(t),t,2) + d/m*diff(x(t),t,1) + c/m*x(t)
17  aL=dsolve(dgl)        #allgemeine Lösung der DGL
18  sL=dsolve(dgl,ics=aw) #spezielle Lösung der DGL
19  x_t = sL.rhs          #rechte Seite der Funktion
20  #Ausgaben
21  print("allgemeine Lösung\n",aL)
22  print("spezielle Lösung\n",sL)
23  print("rechte Seite der Funktion\nx(t) =",x_t)
24  plot(x_t,(t,0,5),ylabel='x(t)')

Listing 12.13     Lösung einer DGL 2. Ordnung mit SymPy

Ausgabe

allgemeine Lösung
 Eq(x(t), (C1*sin(3.14006369362152*t) +  
C2*cos(3.14006369362152*t))*exp(-0.1*t))
spezielle Lösung
 Eq(x(t), (0.15923243882462*sin(3.14006369362152*t) + 5.0*cos(3.14006369362152*t))*exp(-0.1*t))
rechte Seite der Funktion
x(t) = (0.15923243882462*sin(3.14006369362152*t) + 5.0*cos(3.14006369362152*t))*exp(-0.1*t)

Abbildung 12.14     Verlauf der Federpendelschwingungen

Analyse

In den Zeilen 03 und 04 werden die unabhängige Variable t und die abhängige Variable x vereinbart. Die Variable x repräsentiert die Auslenkung der Pendelmasse in Richtung der x-Achse.

In den Zeilen 05 bis 07 können Sie die Daten des Feder-Masse-Systems ändern.

In den Zeilen 10 und 11 stehen die Anfangsbedingungen: die Auslenkung x(0) = 5 und die Anfangsgeschwindigkeit x'(t) = 0.

In Zeile 16 wird die DGL definiert. Die auskommentierten Zeilen 14 und 15 sind alternative Syntaxnotationen.

In Zeile 17 berechnet die Methode dsolve(dgl) die allgemeine Lösung aL der DGL.

In Zeile 18 bewirkt der Parameter ics=aw, dass die Methode dsolve(dgl,ics=aw) die spezielle Lösung der DGL berechnet. Die Abkürzung ics steht für initial conditions.

In Zeile 19 ermittelt die Anweisung x_t=sL.rhs die rechte Seite der Funktionsgleichung und weist sie dem Objekt x_t zu. Die Abkürzung rhs steht für right hand side. Wenn Sie mit der Anweisung print(type(sL.rhs)) den Type von rhs ermitteln, erhalten Sie die Ausgabe <class 'sympy.core.mul.Mul'>.

Die SymPy-Methode plot(x_t,(t,0,5), …) kann mithilfe des Objekts x_t den Funktionsgraphen der Lösung zeichnen (Zeile 24). Das Tupel (t,0,5) enthält die unabhängige Variable t und den Darstellungsbereich des Funktionsgraphen auf der x-Achse.