Computeralgebra mit SymPy#
In diesem Abschnitt wollen wir das Python-Modul SymPy vorstellen, welches uns ermöglicht symbolische Berechnungen durchzuführen. Insbesondere werden wir lernen, wie wir mit mathematischen Ausdrücken arbeiten können, wie wir Funktionen differenzieren und integrieren, und wie wir (nicht)lineare Gleichungssysteme und Differentialgleichungen lösen.
Wir müssen zunächst das Paket installieren, beispielsweise über unsere conda-Umgebung
conda install -c conda-forge sympy
und anschließend in unser Python-Skript einbinden:
import sympy as sy
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
Cell In[1], line 1
----> 1 import sympy as sy
ModuleNotFoundError: No module named 'sympy'
Arbeiten mit Funktionen#
Beginnen wir zunächst damit eine mathematische Funktion zu definieren. Dazu müssen wir die veränderliche Variable kennzeichnen und anschließend die Funktionsvorschrift eingeben:
x = sy.symbols('x')
f = 1./3*x**3 - 2.*x**2 + 3.
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[2], line 1
----> 1 x = sy.symbols('x')
2 f = 1./3*x**3 - 2.*x**2 + 3.
NameError: name 'sy' is not defined
Die Variable x ist vom Typ
type(x)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[3], line 1
----> 1 type(x)
NameError: name 'x' is not defined
und der Ausdruck f vom Typ
type(f)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[4], line 1
----> 1 type(f)
NameError: name 'f' is not defined
Objekte vom Typ sympy.core.add.Add entstehen bei der Summation von, in unserem Fall,
f.args
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[5], line 1
----> 1 f.args
NameError: name 'f' is not defined
Die 3 Argumente selbst sind symbolische Ausdrücke, welche wiederrum aus anderen Rechenoperationen hervorgegangen sind:
g1 = f.args[0]
print("Erster Summand", g1, "vom Typ", type(g1))
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[6], line 1
----> 1 g1 = f.args[0]
2 print("Erster Summand", g1, "vom Typ", type(g1))
NameError: name 'f' is not defined
g2 = f.args[1]
print("Zweiter Summand", g2, "vom Typ", type(g2))
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[7], line 1
----> 1 g2 = f.args[1]
2 print("Zweiter Summand", g2, "vom Typ", type(g2))
NameError: name 'f' is not defined
g3 = f.args[2]
print("Dritter Summand", g3, "vom Typ", type(g3))
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[8], line 1
----> 1 g3 = f.args[2]
2 print("Dritter Summand", g3, "vom Typ", type(g3))
NameError: name 'f' is not defined
Der Ausdruck g1 ist nur die Zahl 3.0 und g2 und g3 sind vom Typ sympy.core.mul.Mult, sind also Ausdrücke die aus der Multiplikation weiterer Ausdrücke hervorgegangen sind. Ein mathematischer Ausdruck wird schließlich in Form eines Binärbaums dargestellt. Die genaue Darstellung bekommen wir wie folgt:
sy.srepr(f)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[9], line 1
----> 1 sy.srepr(f)
NameError: name 'sy' is not defined
Vielleicht haben wir durch diese Betrachtung schon eine grobe Vorstellung wie symbolische Ausdrücke funktionieren. Über die genauen Details wollen wir uns aber hier keine weiteren Gedanken machen.
Schauen wir uns an, was wir mit der eben definierten Funktion anstellen können
Differenzieren#
Df = f.diff(x)
print("Ableitung:")
Df
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[10], line 1
----> 1 Df = f.diff(x)
2 print("Ableitung:")
3 Df
NameError: name 'f' is not defined
bzw. mit der freien Funktion sympy.diff
Df = sy.diff(f,x)
print("Ableitung:")
Df
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[11], line 1
----> 1 Df = sy.diff(f,x)
2 print("Ableitung:")
3 Df
NameError: name 'sy' is not defined
Natürlich kann SymPy auch kompliziertere Funktionen ableiten. Zu beachten ist hier, dass wir für alle mathematischen Funktionen die entsprechenden Funktionen aus dem SymPy-Modul verwenden:
g = sy.exp(-x**2)*sy.sin(x)
g.diff(x)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[12], line 1
----> 1 g = sy.exp(-x**2)*sy.sin(x)
2 g.diff(x)
NameError: name 'sy' is not defined
Integration#
Mit der Methode integrate können wir die Stammfunktion einer Funktion berechnen lassen
F = f.integrate(x)
print("Stammfunktion:")
F
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[13], line 1
----> 1 F = f.integrate(x)
2 print("Stammfunktion:")
3 F
NameError: name 'f' is not defined
Manchmal ist aber auch SymPy überfordert. Bei zu komplizierten Funktionen gibt die integrate-Methode lediglich
G = g.integrate(x)
print("Stammfunktion:")
G
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[14], line 1
----> 1 G = g.integrate(x)
2 print("Stammfunktion:")
3 G
NameError: name 'g' is not defined
zurück. Wir können auch bestimmte Integrale berechnen, indem wir als zweiten Parameter ein Tupel bestehend aus der Integrationsvariablen und den beiden Grenzen angeben:
v = sy.integrate(f, (x,0,1))
print("Integral über f mit Grenzen 0 und 1:")
v
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[15], line 1
----> 1 v = sy.integrate(f, (x,0,1))
2 print("Integral über f mit Grenzen 0 und 1:")
3 v
NameError: name 'sy' is not defined
Um uneigentliche Integrale, beispielsweise \(\int_0^\infty g(x)\,\mathrm dx\) zu berechnen nutzt man als obere Grenze einfach sy.oo (unendlich):
g = x**2*sy.exp(-x)
v = sy.integrate(g,(x,0,sy.oo))
print("Uneigentliches Inegral:")
v
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[16], line 1
----> 1 g = x**2*sy.exp(-x)
2 v = sy.integrate(g,(x,0,sy.oo))
3 print("Uneigentliches Inegral:")
NameError: name 'x' is not defined
Potenzreihen#
SymPy bietet viele weitere nützliche Funktionen, beispielsweise sympy.series(...). Diese Funktion berechnet eine Potenzreihe beliebiger Ordnung für eine gegebene Funktion:
G = g.series(x, 0, 4) # Taylorapproximation im Entwicklungsunkt 0, bis Ordnung 4
print("Potenzreihe:")
G
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[17], line 1
----> 1 G = g.series(x, 0, 4) # Taylorapproximation im Entwicklungsunkt 0, bis Ordnung 4
2 print("Potenzreihe:")
3 G
NameError: name 'g' is not defined
G = G.removeO()
print("Potenzreihe ohne Restglied:")
G
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[18], line 1
----> 1 G = G.removeO()
2 print("Potenzreihe ohne Restglied:")
3 G
NameError: name 'G' is not defined
Plotten#
SymPy stellt auch ein direktes Interface zu MatplotLib bereit. Der Unterschied zu den in visualisation gelernten Funktionen ist, dass wir hier nicht erst ein Gitter für die x-Werte anlegen müssen. Dies geschieht in den plot-Funktionen von SymPy automatisch. Wir müssen lediglich die SymPy-Expression und gegebenenfalls zusätzliche Funktionsparameter der sympy.plot-Funktion übergeben:
p1 = sy.plot(g, (x,-1,1), show=False, legend=True, title="Potenzreihenapproximation")
p2 = sy.plot(G, (x,-1,1), show=False)
p1.append(p2[0])
p1.show()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[19], line 1
----> 1 p1 = sy.plot(g, (x,-1,1), show=False, legend=True, title="Potenzreihenapproximation")
2 p2 = sy.plot(G, (x,-1,1), show=False)
4 p1.append(p2[0])
NameError: name 'sy' is not defined
Übungsaufgabe
Berechne die Ableitung, die Stammfunktion, und erstelle einen Plot der Funktionen
Nichtlineare Gleichungen#
Mit SymPy können wir auch nichtlineare Gleichungen bzw. Gleichungssysteme lösen. Dies geschieht algebraisch, sofern möglich.
Die Funktionen, welche wir dafür benötigen, ist sympy.solve(...). Wir konstruieren ein einfaches Beispiel in dem wir die stationären Punkte der Funktion \(f(x)=x\,\ln(x)\) berechnen:
f = x/(x**2+1)
Df = f.diff()
p3 = sy.plot(f,(x,-4,4)) # Plot function
eq = sy.Eq(Df, 0) # Create equation with lhs DF and rhs 0
sols = sy.solve(eq, x) # Solve the equation
print("Extrema at", sols) # Print solution
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[20], line 1
----> 1 f = x/(x**2+1)
2 Df = f.diff()
3 p3 = sy.plot(f,(x,-4,4)) # Plot function
NameError: name 'x' is not defined
Die stationären Punkte werden in einer Liste abgelegt. Über diese können wir anschließend iterieren. So wie in folgender Probe:
for sol in sols:
print("Stationärer Punkt x =", sol)
print("Ist vom Typ", type(sol))
print("Erfüllt f'(x) =", Df.subs(x, sol))
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[21], line 1
----> 1 for sol in sols:
2 print("Stationärer Punkt x =", sol)
3 print("Ist vom Typ", type(sol))
NameError: name 'sols' is not defined
Auch Systeme von Gleichungen können mit diesem Befehl gelöst werden wie folgendes Beispiel zeigt:
x1, x2 = sy.symbols("x1 x2")
equations = [
sy.Eq(2*x1 + 1*x2, 10),
sy.Eq(1*x1 - 2*x2, 11)
]
sols = sy.solve(equations)
sols
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[22], line 1
----> 1 x1, x2 = sy.symbols("x1 x2")
3 equations = [
4 sy.Eq(2*x1 + 1*x2, 10),
5 sy.Eq(1*x1 - 2*x2, 11)
6 ]
8 sols = sy.solve(equations)
NameError: name 'sy' is not defined
Die zu lösenden Gleichungen wurden hier einfach in eine Liste gepackt und dem solve-Befehl übergeben. Das Ergebnis ist hier ein Objekt vom Typ
type(sols)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[23], line 1
----> 1 type(sols)
NameError: name 'sols' is not defined
Dictionary ist auch eine Containerklasse, welche Key-Value-Paare abspeichert. Der Key (Schlüssel) ist hier der Name der Variable und der Value deren Wert. Wir können auch mit einer for-Schleife über ein Dictionary iterieren. Die lokale Variable a in folgendem Beispiel ist dabei lediglich der Schlüssel. Den zugehörigen Wert bekommt man mit sols[a]:
print("Die Lösung unseres Gleichungssystems lautet")
for a in sols:
print(a, "=", sols[a])
Die Lösung unseres Gleichungssystems lautet
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[24], line 2
1 print("Die Lösung unseres Gleichungssystems lautet")
----> 2 for a in sols:
3 print(a, "=", sols[a])
NameError: name 'sols' is not defined
Alternativ lassen sich lineare Gleichungssysteme auch mit den Befehl sympy.linsolve(...) lösen:
x1, x2, x3 = sy.symbols('x1, x2, x3')
sols = sy.linsolve([2*x1+3*x2-6*x3 - 3, 1*x1-3*x2+3*x3 - 2], (x1,x2,x3))
sols
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[25], line 1
----> 1 x1, x2, x3 = sy.symbols('x1, x2, x3')
2 sols = sy.linsolve([2*x1+3*x2-6*x3 - 3, 1*x1-3*x2+3*x3 - 2], (x1,x2,x3))
3 sols
NameError: name 'sy' is not defined
Für dieses System aus 2 Gleichungen und 3 Unbekannten hat der linsolve-Befehl eine Lösung ermittelt, die noch von \(x_3\) abhängt. Natürlich erwarten wir in diesem Beispiel unendlich viele Lösungen. Mit der subs-Methode können wir die freie Variable mit einem speziellen Wert substituieren:
for sol in sols:
print("Allgemeine Lösung :", sol)
print("Spezielle Lösung :", sol.subs('x3', 0))
print("Spezielle Lösung :", sol.subs('x3', 2))
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[26], line 1
----> 1 for sol in sols:
2 print("Allgemeine Lösung :", sol)
3 print("Spezielle Lösung :", sol.subs('x3', 0))
NameError: name 'sols' is not defined
Gewöhnliche Differentialgleichungen#
Auch einfache Differentialgleichungen können mit SymPy gelöst werden. Zunächst müssen wir dafür eine Differentialgleichung definieren. Dazu legen wir zunächst ein Objekt vom Typ Function an und definieren die Differentialgleichung als mathematischen Ausdruck. Die Funktion Derivative erlaubt es die Ableitungen der gesuchten Funktion einzuarbeiten:
y = sy.Function('y')
ode = sy.Derivative(y(x),x,x) + 2*sy.Derivative(y(x),x) - 6*y(x) - sy.exp(3*x)
ode
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[27], line 1
----> 1 y = sy.Function('y')
2 ode = sy.Derivative(y(x),x,x) + 2*sy.Derivative(y(x),x) - 6*y(x) - sy.exp(3*x)
3 ode
NameError: name 'sy' is not defined
Die Funktion zum Lösen der Differentialgleichung ist sympy.dsolve (das “d” steht für differential equation). Wir übergeben den eben für die ODE definierten Ausdruck sowie die Funktion nach der wir diese lösen sollen:
sol = sy.dsolve(ode, y(x))
sol
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[28], line 1
----> 1 sol = sy.dsolve(ode, y(x))
2 sol
NameError: name 'sy' is not defined
type(sol)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[29], line 1
----> 1 type(sol)
NameError: name 'sol' is not defined
Das Ergebnis ist ein Objekt vom Tyo Equality, bestehend aus einer linken und einer rechten Seite einer Gleichung. Um die rechte Seite zu extrahieren und daraus eine Funktion zu definieren nutzen wir
f = sol.rhs
type(f)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[30], line 1
----> 1 f = sol.rhs
2 type(f)
NameError: name 'sol' is not defined
Unsere eben berechnete Funktion hängt von folgenden Variablen ab:
f.free_symbols
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[31], line 1
----> 1 f.free_symbols
NameError: name 'f' is not defined
Dass unsere Lösung von 2 beliebigen Konstanten C1 und C2 abhängt ist nicht verwunderlich, da wir keine Anfangsbedingungen für unsere Differentialgleichung zweiter Art definiert haben. Wir können nun auch feste Werte für die Konstanten wählen. Dazu definieren wir ein Dictionary und nutzen die Funktion subs (Substituiere):
constants = {'C1': 1, 'C2': 3} # Dictionary anlegen
my_f = f.subs(constants) # Konstanten durch Werte ersetzen
my_f
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[32], line 2
1 constants = {'C1': 1, 'C2': 3} # Dictionary anlegen
----> 2 my_f = f.subs(constants) # Konstanten durch Werte ersetzen
3 my_f
NameError: name 'f' is not defined
Unsere Funktion my_f hängt nun nur noch von x ab und wir können wie gewohnt weiter arbeiten:
sy.plot(my_f, (x,0,2))
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[33], line 1
----> 1 sy.plot(my_f, (x,0,2))
NameError: name 'sy' is not defined
Schauen wir uns nun noch an, wie wir Anfangswertprobleme, also Differentialgleichungen zusammen mit entsprechenden Anfangsbedingungen, lösen können. Ruft man mit sp.dsolve? den Hilfetext zur Funktion dsolve auf, stößt man schnell auf den Parameter ics. Dieser soll ein Dictionary sein, welches alle Anfangs- und Randbedingungen beinhaltet. Dafür soll folgende Syntax genutzt werden:
initial_conditions = {y(0) : 0, # y(0) = 0
y(x).diff(x).subs(x, 0) : 1} # y'(0)= 1
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[34], line 1
----> 1 initial_conditions = {y(0) : 0, # y(0) = 0
2 y(x).diff(x).subs(x, 0) : 1} # y'(0)= 1
NameError: name 'y' is not defined
Wir haben hier die Anfangsbedingungen \(y(0)=0\) und \(y'(0)=1\) vorgegeben. Der Aufruf von dsolve ergibt nun:
sol = sy.dsolve(ode, y(x), ics=initial_conditions)
f = sol.rhs
f
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[35], line 1
----> 1 sol = sy.dsolve(ode, y(x), ics=initial_conditions)
2 f = sol.rhs
3 f
NameError: name 'sy' is not defined
Offensichtlich hängt die Lösung nur noch von der Variablen \(x\) ab. Klar, da die Konstanten in der allgemeinen Lösung durch die Anfangsbedingungen eliminiert wurden.
f aus oberem Block ist wieder ein mathematischer Ausdruck, mit dem wir wie gewohnt weiter arbeiten können:
print("f ist vom Typ", type(f))
sy.plot(f, (x,0,1))
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[36], line 1
----> 1 print("f ist vom Typ", type(f))
2 sy.plot(f, (x,0,1))
NameError: name 'f' is not defined
In folgendem Beispiel wird die Riccati-Differentialgleichung
gelöst. Interessant ist auch die Ausgabe der Funktion sympy.classify_ode(...). SymPy analysiert die Differentialgleichung und wendet, abhängig von der Klassifikation der Gleichng, eine entsprechende Lösungsstrategie an:
ode = sy.Derivative(y(x),x) + y(x)**2 + 2/x*y(x) + 2/(x**2)
sy.classify_ode(ode)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[37], line 1
----> 1 ode = sy.Derivative(y(x),x) + y(x)**2 + 2/x*y(x) + 2/(x**2)
2 sy.classify_ode(ode)
NameError: name 'sy' is not defined
sol = sy.dsolve(ode, y(x))
sol
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[38], line 1
----> 1 sol = sy.dsolve(ode, y(x))
2 sol
NameError: name 'sy' is not defined
f = sol.rhs
my_f = f.subs('C1', 1)
my_f
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[39], line 1
----> 1 f = sol.rhs
2 my_f = f.subs('C1', 1)
3 my_f
NameError: name 'sol' is not defined