Symbolische Mathematik mit SymPy

Das Modul SymPy verfügt zwar über ähnliche Leistungsmerkmale wie Maple, Mathematica oder Maxima, jedoch erreicht es nicht deren vollen Funktionsumfang. Gegenüber diesen Computeralgebrasystemen (CAS) hat es aber den besonderen Vorteil, dass SymPy-Programme vergleichbar einfacher und flexibler erstellt werden können.

In der Tabelle sind häufig genutzte SymPy-Befehle aufgelistet.

Methode Beschreibung
apart(p) Zerlegt das Polynom p in seine Partialbrüche.
cancel(p) Erzeugt aus Partialbrüchen p eine ganzrationale Funktion.
diff(f,x,k) Berechnet die k-te Ableitung einer Funktion f.
dsolve(eq,f(x) Löst eine gewöhnliche Differenzialgleichung eq für die Funktion f(x).
epand(T) Berechnet den Term T.
integrate(f,x) Berechnet die Stammfunktion F(x) der Funktion f(x).
limit(y,x,0) Berechnet den Grenzwert einer Funktion.
simplify(term) Vereinfacht einen Term.

Differenzieren und Integrieren

Mit dem Befehl diff(y,x) kann die Ableitung einer Funktion ermittelt werden. Das Integral einer Funktion wird mit dem Befehl integrate(y,x) symbolisch berechnet. Die Befehle können direkt in die Phyton Shell eingegeben werden.

>>> from sympy import *
>>> x=symbols('x')
>>> y=x**3-2*x**2+10
>>> diff(y,x)
3*x**2 - 4*x
>>> diff(y,x,2)
2*(3*x - 2)
>>> integrate(y,x)
x**4/4 - 2*x**3/3 + 10*x

Laplace-Transformation

Die inverse Laplace-Transformation wird mit der Methode
inverse_laplace_transform(U2_s,s,t)
ermittelt. Die Übertragungsfunktion H(s) ist als algebraischer Ausdruck gegeben.

Das Beispiel zeigt, wie die inverse Laplace-Transformation die Sprungantwort eines R-L-C-Gliedes berechnet.

sprungantwort r-c-l-glied

Die Übertragungsfunktion kann aus der Schaltung entwickelt werden:

\[H\left( s\right) =\frac{U_{2}\left( s\right) }{U_{1}\left( s\right) } =\frac{\frac{1}{Cs} }{R+Ls+\frac{1}{Cs} } =\frac{1}{LCs^{2}+RCs+1} \]

Mit der Sprungfunktion aus dem Bildbereich erhält man die Ausgangsspannung für den Bildbereich:

\[U_{2}\left( s\right) =\frac{U_{1}}{s} \cdot \frac{1}{LCs^{2}+RCs+1} \]

Die Sprungantwort wird in Zeile #2 berechnet.

#Sprungantwort mit Laplace-Transformation
from sympy import *
#from sympy.integrals import inverse_laplace_transform
s,t,C,L,R = symbols("s t C L R",positive=True)
U1=10
R,L,C=2,1/2,1/4
#R,L,C=2,1,1
#R,L,C=1,4,8/9
U1_s=U1/s
Z_s=R+L*s+1/(C*s)
I_s=U1_s/Z_s
H_s=1/(R + L*s + 1/(C*s))/(C*s)         #1
#H_s=expand(H_s) #Übertragungsfunktion
U2_s=U1_s*H_s    #Sprungantwort
uc=inverse_laplace_transform(U2_s,s,t)  #2
ic=inverse_laplace_transform(I_s,s,t)   #3
#Ausgaben
print("Übertragungsfunktion")
print(H_s)
print("Spannung am Kondensator")
print("uc=",simplify(uc))
print("Kondensatorstrom")
print("ic=",ic)
plt=plot(uc,ic,(t, 0, 5),show=False)
plt[0].line_color = 'b'
plt[1].line_color = 'r'
plt.show()

Ausgabe

Übertragungsfunktion
4.0/(s*(0.5*s + 2 + 4.0/s))
Spannung am Kondensator
uc= 10.0 - 10.0*sqrt(2)*exp(-2.0*t)*sin(2.0*t + pi/4)
Kondensatorstrom
ic= 10.0*exp(-2.0*t)*sin(2.0*t)
Spannungs- und stromverlauf R-C-L-Glied

Analyse

In Zeile #1 steht die Übertragungsfunktion des R-L-C-Gliedes. In den Zeilen #2 und #3 berechnet die inverse Laplace-Transformation die Ausgangsspannung und den Kondensatorstrom des Zweitors (Übertragungssystem).

Federpendel

Eine Schraubenfeder ist mit einer Masse verbunden.

Federpendel waagerecht

Wenn die Masse ausgelenkt wird, schwingt sie hin und her.

Die Schwingungen des Federpendels können durch folgende DGL

\[\ddot{x}+\frac{d}{m} \cdot \dot{x}+\frac{c}{m} \cdot x=0\]

beschrieben werden. Dabei steht d für die Dämpfung.

Gelöst wird die DGL mit der SymPy-Methode dsolve(dgl,ics=aw). Die Abkürzung ics steht für initial/boundary conditions (dt.: Anfangs- bzw. Randbedingungen). Die Abkürzung rhs steht für right hand side also für die rechte Seite einer Gleichung.

#sympy_federpendel.py
from sympy import * 
t = symbols('t')
x = Function('x')
m=1.0   #Masse in kg
d=0.2   #Dämpfung 
c=9.87  #Federkonstante in N/m
x0=5    #Anfangswert x(0)
x_0=0   #x'(0)
dgl=diff(x(t),t,2) + d/m*diff(x(t),t,1) + c/m*x(t)
#Anfangswerte
aw={
    x(0):x0,
    x(t).diff(t,1).subs(t,0):x_0
    }
aL=dsolve(dgl)        #allgemeine Lösung der DGL
sL=dsolve(dgl,ics=aw) #spezielle Lösung der DGL
gL=sL.rhs   #rechte Seite der speziellen Lösung
print("allgemeine Lösung\n",aL)
print("spezielle Lösung\n",sL)
print("Funktionsterm\n x(t) =",gL)
p=plot(gL,(t,0,5),ylabel='x(t)',show=False)
#p.save('federpendel.svg')
p.show()

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))
Funktionsterm
x(t) = (0.15923243882462*sin(3.14006369362152*t) + 5.0*cos(3.14006369362152*t))*exp(-0.1*t)

Federpendel Schwingungsferlauf

Das Beispiel stammt aus meinem Buch "Mathematische Algorithmen mit Python".

Gekoppeltes Federpendel

Zwei Massen sind mit zwei Schraubenfedern verbunden.

Feder-Masse-System

Wenn eine oder auch beide Massen ausgelenkt werden, schwingen die Massen in x-Richtung hin und her. Die Dämpfung soll vernachlässigt werden.

Mit der Abkürzung

\[k=\frac{c}{m}\]

wird das Feder-Masse-System durch das folgende DGL-System beschrieben:

\[\ddot{x}_{1}=-2k\cdot x_{1}+k\cdot x_{2}\]

\[\ddot{x}_{2}= k\cdot x_{1}-k\cdot x_{2}\]

Im Gegensatz zur numerischen Lösung mit der SciPy-Methode solve_ivp(...) muss man dieses DGL-System nicht in ein System 1. Ordnung umwandeln, wenn man die SymPy-Methode dsolve_system(...) nutzt. Diese Methode kann nur lineare Differenzialgleichungen lösen.

#gekoppeltes_federpendel.py
from sympy import symbols,Eq,Function,plot,N
from sympy.solvers.ode.systems import dsolve_system
t = symbols("t")
x1 = Function("x1")
x2 = Function("x2")
m=2  #Masse in kg
c=1  #Federkonstante in N/m
k=c/m #Abkürzung
#DGL-System
dgl1=Eq(x1(t).diff(t,2),-2*k*x1(t) + k*x2(t)) 
dgl2=Eq(x2(t).diff(t,2), +k*x1(t) - k*x2(t))
#Anfangswerte
aw={
    x1(0): 1, 
    x2(0): -1,
    x1(t).diff(t,1).subs(t,0):0,
    x2(t).diff(t,1).subs(t,0):0
    }
#Lösung des DGL-Systems
gleichungen = [dgl1,dgl2]
aL=dsolve_system(gleichungen)        #allgemeine Lösung
sL=dsolve_system(gleichungen,ics=aw) #spezielle Lösung
gX1=sL[0][0].rhs 
gX2=sL[0][1].rhs
#Ausgaben
print("allgemeine Lösung\n",aL)
print("spezielle Lösung")
print("x1(t) =",N(gX1,3))
print("x2(t) =",N(gX2,3))
p=plot(gX1,gX2,(t,0,20),show=False,legend=True)
p.title='gekoppeltes Federpendel'
p.xlabel='t'
p.ylabel='Amplitude'
p[0].line_color='blue'
p[0].label='x1'
p[1].line_color='green'
p[1].label='x2'
#p.save('gekoppeltes_federpendel.png')
#p.save('gekoppeltes_federpendel.svg')
p.show()

Ausgabe

allgemeine Lösung
[[Eq(x1(t), 1.4142135623731*C1*sin(0.437016024448821*t) - 1.4142135623731*C2*sin(1.14412280563537*t) + 1.4142135623731*C3*cos(0.437016024448821*t) - 1.4142135623731*C4*cos(1.14412280563537*t)), Eq(x2(t), 2.28824561127074*C1*sin(0.437016024448821*t) + 0.874032048897642*C2*sin(1.14412280563537*t) + 2.28824561127074*C3*cos(0.437016024448821*t) + 0.874032048897642*C4*cos(1.14412280563537*t))]]
spezielle Lösung
x1(t) = -0.171*cos(0.437016024448821*t) + 1.17*cos(1.14412280563537*t)
x2(t) = -0.276*cos(0.437016024448821*t) - 0.724*cos(1.14412280563537*t)
Schwingungsverlauf

Dieses Beispiel stammt nicht aus dem Buch "Der Python-Kurs für Ingenieure und Naturwissenschaftler". Stattdessen wird in der 2. Auflage ein Beispiel aus der physikalischen Chemie vorgestellt. Mit der SymPy-Methode dsolve_system([dgl1,dgl2,dgl3]) wird ein DGL-System mit 3 Differenzialgleichungen für eine Folgereaktion A → B → C symbolisch und numerisch gelöst. Der Verlauf der chemischen Reaktion wird auch grafisch dargestellt.

Besuchen Sie die Seite von SymPy.

StartseiteNumPyMatplotlibSciPyVPython