Computeralgebra mit SymPy - Mathematical Programming

In diesem Abschnitt wollen wir das Python-Modul SymPy vorstellen, welches uns ermöglicht, symbolische Berechnungen durchzuführen. Insbesondere werden wir lernen,

mit mathematischen Ausdrücken zu arbeiten,
Funktionen zu differenzieren und zu integrieren,
(nichtlineare) Gleichungssysteme zu lösen,
sowie einfache Differentialgleichungen zu behandeln.

Zunächst müssen wir das Paket installieren, beispielsweise über unsere conda-Umgebung:

conda install -c conda-forge sympy

Anschließend binden wir das Modul in unser Python-Skript ein:

import sympy as sy

Arbeiten mit Funktionen¶

Beginnen wir zunächst damit eine mathematische Funktion zu definieren. Dazu müssen wir zunächst eine symbolische Variable deklarieren und anschließend die Funktionsvorschrift angeben:

x = sy.symbols('x')
f = x**3/3 - 2*x**2 + 3

Die Variable x ist vom Typ

type(x)

sympy.core.symbol.Symbol

und der Ausdruck f vom Typ

type(f)

sympy.core.add.Add

Objekte vom Typ sympy.core.add.Add entstehen durch die Addition mehrerer Ausdrücke. Die einzelnen Summanden können wir über das Attribut args abfragen:

f.args

(3, -2*x**2, x**3/3)

Die drei Argumente selbst sind wiederum symbolische Ausdrücke, welche aus weiteren Rechenoperationen hervorgegangen sind:

g1 = f.args[0]
print("Erster Summand", g1, "vom Typ", type(g1))

g2 = f.args[1]
print("Zweiter Summand", g2, "vom Typ", type(g2))

g3 = f.args[2]
print("Dritter Summand", g3, "vom Typ", type(g3))

Erster Summand 3 vom Typ <class 'sympy.core.numbers.Integer'>
Zweiter Summand -2*x**2 vom Typ <class 'sympy.core.mul.Mul'>
Dritter Summand x**3/3 vom Typ <class 'sympy.core.mul.Mul'>

Der Ausdruck g1 ist lediglich die Zahl 3.0, während g2 und g3 vom Typ sympy.core.mul.Mult sind. Diese entstehen durch die Multiplikation symbolischer Ausdrücke.

Ein mathematischer Ausdruck wird intern als Baumstruktur (Ausdrucksbaum) dargestellt. Die vollständige interne Darstellung können wir uns mit folgendem Befehl anzeigen lassen:

sy.srepr(f)

"Add(Mul(Rational(1, 3), Pow(Symbol('x'), Integer(3))), Mul(Integer(-1), Integer(2), Pow(Symbol('x'), Integer(2))), Integer(3))"

Vielleicht haben wir durch diese Betrachtung bereits eine grobe Vorstellung davon bekommen, wie symbolische Ausdrücke intern aufgebaut sind. Auf die genauen Details wollen wir hier jedoch nicht weiter eingehen.

Schauen wir uns nun an, was wir mit der eben definierten Funktion alles anstellen können.

Auswerten symbolischer Ausdrücke¶

Symbolische Ausdrücke lassen sich auch für konkrete Werte der Variablen auswerten. Dazu verwenden wir die Methode subs (substitute), mit der wir eine symbolische Variable durch einen numerischen Wert ersetzen können.

f.subs(x, 2)

Hier wird die Variable x im Ausdruck f durch den Wert 2 ersetzt und das Ergebnis berechnet.

Wir können auch mehrere Variablen gleichzeitig ersetzen. Dazu übergeben wir ein Dictionary mit den gewünschten Ersetzungen:

x, y = sy.symbols("x y")
g = x**2 + y**2

r = g.subs({x:1, y:2})

print(f"r = {r} ({type(r)})")

r = 5 (<class 'sympy.core.numbers.Integer'>)

Das Ergebnis dieser Operation ist weiterhin ein symbolischer Ausdruck. Falls wir eine numerische Auswertung als Fließkommazahl wünschen, können wir zusätzlich die Methode evalf() verwenden:

r.evalf()

Vereinfachen von Ausdrücken¶

SymPy kann mathematische Ausdrücke auch vereinfachen und umformen. Dies ist besonders hilfreich, wenn symbolische Rechnungen zu komplizierten Zwischenergebnissen führen.

Betrachten wir zunächst folgenden Ausdruck:

f = (x**2 - 1)/(x-1)

Mathematisch wissen wir, dass

\frac{x^2-1}{x-1} = x+1

(1)

für alle $x\ne 1$ gilt. Mit der Funktion simplify kann SymPy versuchen, solche Vereinfachungen automatisch durchzuführen:

sy.simplify(f)

Mit expand können Produkte ausmultipliziert werden:

f = (x+1)**3
sy.expand(f)

Die Funktion factor versucht dagegen einen Ausdruck zu faktorisieren:

f = x**2 - 1
sy.factor(f)

Diese Funktionen sind besonders hilfreich, wenn wir symbolische Rechnungen weiterverarbeiten oder mathematische Strukturen besser sichtbar machen wollen.

Differenzieren¶

Um die Ableitung eines symbolischen Ausdrucks bezüglich einer symbolischen Variable zu berechnen, können wir entweder die Methode diff verwenden, die allen SymPy-Ausdrücken zur Verfügung steht,

f = x**3/3 - 2*x**2 + 3
Df = f.diff(x)
print("Ableitung:")
Df

Ableitung:

oder die freie Funktion sympy.diff verwenden:

Df = sy.diff(f,x)
print("Ableitung:")
Df

Ableitung:

Beide Varianten liefern dasselbe Ergebnis.

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)

Höhere Ableitungen und partielle Ableitungen¶

Mit SymPy lassen sich auch höhere Ableitungen eines symbolischen Ausdrucks berechnen.
Dazu kann man entweder die Methode .diff mehrfach aufrufen oder direkt einen zweiten Parameter angeben, der angibt, wie oft abgeleitet werden soll.

Beispiel:

x = sy.symbols('x')

# Erste Ableitung
Df1 = f.diff(x)
print("1. Ableitung:", Df1)

# Zweite Ableitung
Df2 = f.diff(x, 2)
print("2. Ableitung:", Df2)

# Dritte Ableitung
Df3 = f.diff(x, 3)
print("3. Ableitung:", Df3)

1. Ableitung: x**2 - 4*x
2. Ableitung: 2*(x - 2)
3. Ableitung: 2

SymPy erlaubt es, auch Funktionen mit mehreren Variablen symbolisch abzuleiten. Wir können partielle Ableitungen nach jeder Variablen berechnen, sowie gemischte Ableitungen.

x, y = sy.symbols('x y')

# Funktion definieren
h = x**2 * y + sy.sin(x*y)

# Partielle Ableitung nach x
df_dx = h.diff(x)
print("Partielle Ableitung nach x:", df_dx)

# Partielle Ableitung nach y
df_dy = h.diff(y)
print("Partielle Ableitung nach y:", df_dy)

# Zweite Ableitung nach x
d2f_dx2 = h.diff(x, 2)
print("Zweite Ableitung nach x:", d2f_dx2)

# Gemischte Ableitung: zuerst x, dann y
d2f_dxdy = h.diff(x).diff(y)
print("Gemischte Ableitung d²f/(dx dy):", d2f_dxdy)

# Alternativ kann man direkt beide Ableitungen in einem Schritt angeben
d2f_dxdy_alt = h.diff(x, y)
print("Gemischte Ableitung (direkt):", d2f_dxdy_alt)

Partielle Ableitung nach x: 2*x*y + y*cos(x*y)
Partielle Ableitung nach y: x**2 + x*cos(x*y)
Zweite Ableitung nach x: y*(-y*sin(x*y) + 2)
Gemischte Ableitung d²f/(dx dy): -x*y*sin(x*y) + 2*x + cos(x*y)
Gemischte Ableitung (direkt): -x*y*sin(x*y) + 2*x + cos(x*y)

Integration¶

Mit der Methode integrate können wir die Stammfunktion einer Funktion berechnen lassen:

F = f.integrate(x)
print("Stammfunktion:")
F

Stammfunktion:

Manchmal ist SymPy bei zu komplizierten Funktionen überfordert. Dann gibt die integrate-Methode einfach den symbolischen Ausdruck zurück:

G = g.integrate(x)
print("Stammfunktion:")
G

Stammfunktion:

zurück.

Bestimmte Integrale:

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

Integral über f mit Grenzen 0 und 1:

Uneigentliche Integrale:

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 Integral:")
v

Uneigentliches Integral:

Potenzreihen¶

SymPy bietet viele nützliche Funktionen für die Arbeit mit Reihen, insbesondere sympy.series(...). Mit dieser Funktion können wir die Taylor- oder Potenzreihe einer Funktion beliebiger Ordnung berechnen:

G = g.series(x, 0, 4) # Taylorapproximation im Entwicklungsunkt 0, bis Ordnung 4
print("Potenzreihe:")
G

Potenzreihe:

Optional können wir das Restglied $O(x^n)$ entfernen, um nur die exakten Terme der Reihe zu erhalten:

G = G.removeO()
print("Potenzreihe ohne Restglied:")
G

Potenzreihe ohne Restglied:

Plotten¶

SymPy stellt auch ein direktes Interface zu MatplotLib bereit. Der Unterschied zu den in Visualisierung mit Matplotlib 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()

Nichtlineare Gleichungen¶

Mit SymPy können wir auch nichtlineare Gleichungen bzw. Gleichungssysteme algebraisch lösen, sofern möglich.

Die zentrale Funktion hierfür ist sympy.solve(...). Wir betrachten ein Beispiel, in dem wir die stationären Punkte der Funktion

f(x) = \frac{x}{x^2+1}

(3)

berechnen.

f = x/(x**2+1)
Df = f.diff()
p3 = sy.plot(f,(x,-4,4))  # Funktion plotten

eq = sy.Eq(Df, 0)         # Gleichung mit linker Seite DF und rechter Seite 0 erstellen
sols = sy.solve(eq, x)    # Gleichung lösen
print("Extrema bei", sols) # Lösung ausgeben

Extrema bei [-1, 1]

Die stationären Punkte werden als Liste zurückgegeben. Wir können nun über diese iterieren, prüfen, um welche Werte es sich handelt, und die Ableitung an diesen Stellen evaluieren:

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))

Stationärer Punkt x = -1
Ist vom Typ <class 'sympy.core.numbers.NegativeOne'>
Erfüllt f'(x) = 0
Stationärer Punkt x = 1
Ist vom Typ <class 'sympy.core.numbers.One'>
Erfüllt f'(x) = 0

Auch Systeme von Gleichungen können mit sympy.solve gelöst werden. Ein Beispiel:

# Symbole definieren
x1, x2 = sy.symbols("x1 x2")

# Gleichungssystem als Liste
equations = [
    sy.Eq(2*x1 + 1*x2, 10),
    sy.Eq(1*x1 - 2*x2, 11) 
]

# System lösen
sols = sy.solve(equations)
sols

{x1: 31/5, x2: -12/5}

Hierbei wird das Ergebnis als Dictionary

type(sols)

dict

zurückgegeben. In einem Dictionary werden Key-Value-Paare gespeichert: Der Key ist der Name der Variable, der Value deren Wert.

Wir können über die Dictionary-Einträge iterieren, um die Lösungen geordnet auszugeben:

print("Die Lösung unseres Gleichungssystems lautet")
for a in sols:
    print(a, "=", sols[a])

Die Lösung unseres Gleichungssystems lautet
x1 = 31/5
x2 = -12/5

Alternativ lassen sich lineare Gleichungssysteme auch mit sympy.linsolve lösen. Beispiel:

# Symbole definieren
x1, x2, x3 = sy.symbols('x1, x2, x3')

# Lineares Gleichungssystem lösen
sols = sy.linsolve([2*x1+3*x2-6*x3 - 3, 1*x1-3*x2+3*x3 - 2], (x1,x2,x3))
sols

Da wir hier 2 Gleichungen und 3 Unbekannte haben, erwartet man unendlich viele Lösungen. linsolve gibt diese allgemein in Abhängigkeit der freien Variablen zurück. Mit der Methode subs kann man eine freie Variable durch einen bestimmten Wert ersetzen:

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))

Allgemeine Lösung : (x3 + 5/3, 4*x3/3 - 1/9, x3)
Spezielle Lösung  : (5/3, -1/9, 0)
Spezielle Lösung  : (11/3, 23/9, 2)

Gewöhnliche Differentialgleichungen¶

Auch einfache Differentialgleichungen (ODEs) können mit SymPy gelöst werden. Zunächst müssen wir dafür die gesuchte Funktion als Function definieren und die Differentialgleichung selbst als symbolischen Ausdruck konstruieren. Die Methode Derivative erlaubt es uns, Ableitungen einzubauen:

x = sy.symbols('x')
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

Die Differentialgleichung lösen wir mit dsolve:

sol = sy.dsolve(ode, y(x))
print(type(sol))
sol

<class 'sympy.core.relational.Equality'>

Das Ergebnis ist ein Equality-Objekt mit linker und rechter Seite. Die rechte Seite können wir extrahieren und daraus eine normale SymPy-Funktion erstellen:

f = sol.rhs
type(f)

sympy.core.add.Add

Unsere eben berechnete Funktion hängt von folgenden Variablen ab:

f.free_symbols

{C1, C2, x}

Da wir keine Anfangsbedingungen definiert haben, enthält die Lösung zwei Konstanten C1 und C2. Diese können wir mit einem Dictionary und subs durch konkrete Werte ersetzen:

constants = {'C1': 1, 'C2': 3} # Dictionary anlegen
my_f = f.subs(constants)       # Konstanten durch Werte ersetzen
my_f

Nun hängt my_f nur noch von x ab, und wir können sie wie gewohnt plotten:

sy.plot(my_f, (x,0,2))

<sympy.plotting.backends.matplotlibbackend.matplotlib.MatplotlibBackend at 0x79ffad86e710>

Schauen wir uns nun noch an, wie wir Anfangswertprobleme, also Differentialgleichungen zusammen mit entsprechenden Anfangsbedingungen, lösen können. Ruft man mit

sy.dsolve?

den Hilfetext zur Funktion dsolve auf, stößt man schnell auf den Parameter ics. Dieser erwartet ein Dictionary, welches die Anfangs- oder Randbedingungen enthält. Dabei wird folgende Syntax verwendet:

initial_conditions = {y(0)                    : 0, # y(0) = 0
                      y(x).diff(x).subs(x, 0) : 1} # y'(0)= 1

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

Offensichtlich hängt die Lösung nur noch von der Variablen x ab. Das ist erwartbar, da die Konstanten der allgemeinen Lösung durch die Anfangsbedingungen bestimmt wurden.

f aus dem obigen Block ist wieder ein mathematischer Ausdruck, mit dem wir wie gewohnt weiterarbeiten können:

print("f ist vom Typ", type(f))
sy.plot(f, (x,0,1))

f ist vom Typ <class 'sympy.core.add.Add'>

<sympy.plotting.backends.matplotlibbackend.matplotlib.MatplotlibBackend at 0x79ffad86dd10>

In folgendem Beispiel wird die Riccati-Differentialgleichung

y'(x) + y(x)^2+\frac{2}{x}\,y(x) = -\frac{2}{x^2}

(4)

gelöst. Interessant ist auch die Ausgabe der Funktion sympy.classify_ode(...). SymPy analysiert dabei die Differentialgleichung und wendet - abhängig von der Klassifikation - eine passende Lösungsstrategie an:

ode = sy.Derivative(y(x),x) + y(x)**2 + 2/x*y(x) + 2/(x**2)
sy.classify_ode(ode)

('1st_rational_riccati',
 'Riccati_special_minus2',
 'separable_reduced',
 'lie_group',
 'separable_reduced_Integral')

sol = sy.dsolve(ode, y(x))
sol

f = sol.rhs
my_f = f.subs('C1', 1)
my_f