Das Computeralgebrasystem: Python Jupyter + SymPy¶
In diesem Jupyter Notebook wird die Programmiersprache Python, in Verbindung mit der Python-Bibliothek SymPy benutzt, um symbolische Berechnungen durchzuführen (Matrix-Berechnungen, Differentiation, Integration, analytisches Lösen von Differentialgleichungen). In diesem Sinne agiert das Jupyter Notebook somit als ein Computeralgebrasystem - ein Computerprogramm, das der Bearbeitung algebraischer Ausdrücke dient. Es löst nicht nur mathematische Aufgaben mit festen Zahlenwerten, sondern auch solche mit symbolischen Ausdrücken (Gleichungen von Variablen, Funktionen, Polynomen und Matrizen) und ist dadurch im Erscheinungsbild den kommerziellen Software-Produkten Mathematica oder Maple sehr ähnlich.
Zunächst wird das Python Modul "sympy" eingebunden. Falls Sie das "sympy" Modul das erste Mal verwenden, müssen Sie es zunächst in Ihrer Python 3 Umgebung installieren (z.B. in einem Linux Terminal mit "pip3 install sympy").
from sympy import *
init_printing()
Berechnungen mit Matrizen¶
Wir definieren uns eine (2x3)-Matrix $\hat{\bf {A}}$ und eine (3x3)-Matrix $\hat{\bf {B}}$:
A = Matrix([[1, 2, 3, 4, 5],[21,22,23,24,25],[51,52,53,54,55]])
A
B = Matrix([[1, 2, 3],[4, 5, 6],[7, 8, 9],[10,11,12],[13,14,15]])
B
Das Produkt der Matrizen $\hat{\bf {A}}$ und $\hat{\bf {B}}$ berechnet sich zu
C = A*B
C
, die transponierte Matrix $(\hat{\bf {A}})^T$
transpose(A)
Man kann die einzelnen Elemente der Matrix auch mittels symbolischer Variablen benennen und dann allgemeine Ausdrücke erhalten.
a,b,c,d = symbols('a,b,c,d')
F = Matrix([[a, b],[c,d]])
F
Z.B. ergibt sich für die Determinante der Matrix $\hat{\bf {F}}$
det(F)
Berechnungen mit Vektoren können in gleicher Weise durchgeführt werden und die Operationen des Skalar- und Kreuzproduktes sind ebenfalls in SymPy schon vordefiniert. Nehmen wir z.B. die folgenden Vektoren $\vec{a}$ und $\vec{b}$ und berechnen das Skalar- und Kreuzprodukt:
$$ \begin{equation} \vec{a} =\left( \begin{array}{ccc} 4 \\ 7 \\ 9 \end{array} \right) \quad , \quad \vec{b} =\left( \begin{array}{ccc} -3 \\ 5 \\ -7 \end{array} \right) \qquad , \quad \hbox{Skalarprodukt:} \,\, s = \vec{a} \cdot \vec{b} \quad , \quad \hbox{Kreuzprodukt:} \,\, \vec{c} = \vec{a} \times \vec{b} \end{equation} $$
va = Matrix([4, 7, 9])
vb = Matrix([-3, 5, -7])
Das Skalarprodukt ($s = \vec{a} \cdot \vec{b} = \sum_{i=1}^3 a_i \, b_i $) ergibt sich mittels dem Befehl $\vec{a}$.dot( $\vec{b}$ ):
va.dot(vb)
Das Kreuzprodukt ( $\vec{c} = \vec{a} \times \vec{b} = \sum_{i,j,k=1}^3 \epsilon_{ijk} \, a_i \, b_j \, \vec{e}_k$ ) ergibt sich mittels dem Befehl $\vec{a}$.cross( $\vec{b}$ ), wobei $\epsilon_{ijk}$ der total antisymmetrischer Tensor (auch Epsilon-Tensor bzw. Levi-Civita-Symbol) ist.
va.cross(vb)
Weitere Matrix Operationen finden Sie unter http://docs.sympy.org/latest/tutorial/matrices.html
Symbolische Berechnungen von Gleichungen¶
Das Modul sympy ist auch dafür gemacht, symbolische Berechnungen durchzuführen. Betrachten wir z.B. die folgende Gleichung $0=x^2 + p \, x + q$. Zunächst definieren wir die symbolischen Variablen der Gleichung und geben der Gleichung den Variablennamen "Gl".
x = symbols('x')
p,q = symbols('p,q')
Gl = Eq(x**2 + p*x + q, 0)
Gl
Die Lösungen erhält man mittels des SymPy-Befehls "solve(...)"
Loes_Gl=solve(Gl,x)
Loes_Gl
Wir legen nun die symbolischen Variablen $p$ und $q$ auf die folgenden Werte fest (SymPy-Befehl "subs(...)", $0=x^2-3x+9$)
Gl_1=Gl.subs(p,-3).subs(q,9)
Gl_1
und lösen die Gleichung erneut:
Loes_Gl_1=solve(Gl_1)
Loes_Gl_1
Die Gleichung besitzt keine reellwertige Lösung, wie man auch anhand der folgenden Abbildung erkennt.
plot(Gl_1.lhs, xlim=(-5,7), ylim=(-2,50));
Differentation und Integration¶
Die erste Ableitung einer Funktion $f(x)$:
hier speziell $f(x)=10 \cdot e^{-x/5} \cdot \hbox{sin}(3 x)$
f = 10 * exp(-x/5) * sin(3*x)
f
Erste Ableitung: $f^\prime(x)$=
f.diff(x)
Zweite Ableitung: $f^{\prime \prime}(x)$=
f.diff(x,x)
Dritte Ableitung: $f^{\prime \prime \prime}(x)$=
f.diff(x,3)
Das unbestimmte Integral einer Funktion $f(x)$: $\int f(x) \, dx$
f.integrate(x)
Das bestimmte Integral einer Funktion $f(x)$: $\int_a^b f(x) \, dx$
f.integrate((x, a, b))
Das bestimmte Integral einer Funktion $f(x)$: $\int_1^2 f(x) \, dx$
f.integrate((x, 1, 2))
f.integrate((x, 1, 2)).evalf()
Analytisches Lösen von Differentialgleichungen¶
Zusätzlich ist es in sympy auch möglich Differentialgleichungen (DGLs) zu lösen (falls eine analytische Lösung existiert). Betrachten wir z.B. die folgende DGL zweiter Ordnung: $\frac{d^2y(x)}{dx^2}=-3y(x)$
y = Function('y')(x)
DGL = Eq(y.diff(x,x), -3*y)
DGL
Die allgemeine Lösung lautet:
dsolve(DGL)
Die allgemeine Lösung hängt von den zwei Parametern $C_1$ und $C_2$ ab, die durch eine Wahl der Anfangs/Randbedingungen festgelegt werden können. Wir wählen z.B. $\quad y(0)=1 \quad $ und $ \quad \left. \frac{dy(x)}{dx}\right|_{x=0}=2\quad $.
Eq(y.subs(x,0),1)
Eq(y.diff(x).subs(x, 0),2)
Die spezielle Lösung lautet dann:
Loes = dsolve(DGL,ics={y.subs(x,0):1, y.diff(x).subs(x, 0): 2})
Loes
Die Lösung können wir wie folgt, direkt mittels SymPy visualisieren:
plot(Loes.rhs, (x, 0, 15));
Die Lösung können wir uns auch mittels "matplotlib" visualisieren, wobei man dazu den analytischen Ausdruck mittels "lambdify(...)" in eine numerische Funktion umwandelt.
import matplotlib.pyplot as plt
import numpy as np
xvals = np.linspace(0, 15, 300)
func = lambdify(x, Loes.rhs)
plt.ylabel('y(x)')
plt.xlabel('x')
plt.plot(xvals,func(xvals));
Die oben definierte numerische Funktion mit dem Namen "func", gibt dabei die Funktionswerte der jeweiligen Gitterpunkte aus. In der folgenden Ausgabe z.B. für die ersten 10 Gitterpunkte.
func(xvals[:10])
array([1. , 1.09643549, 1.18459782, 1.26382176, 1.33350954,
1.39313532, 1.44224919, 1.48048057, 1.50754098, 1.52322624])
Legen wir hingegen nur eine der Randbedingungen fest (z.B. $y(0)=1$), so können wir in einer Animation die Auswirkungen der anderen Randbedingung uns visualisieren. Hierzu führt man einen Animationsparameter ein, z.B. den Parameter $\alpha = \left. \frac{dy(x)}{dx}\right|_{x=0}$, und berechnet sich den analytischen Ausdruck der Lösung der DGL.
alpha = symbols("alpha")
Loes_alpha = dsolve(DGL,ics={y.subs(x,0):1, y.diff(x).subs(x, 0): alpha})
Loes_alpha
Nun kann man z.B. direkt die analytische Lösung mittels der spb-Python-Bibliothek (das SymPy Plotting Backends-Modul) als eine interaktive Animation darstellen (siehe SPB Modules Reference) und speichern. Wir animieren im Folgenden die Lösung im Bereich $\alpha \in [-2,6]$.
from spb import plot as spb_plot
from IPython.display import IFrame
anim = spb_plot( Loes_alpha.rhs, (x, 0, 15),
params={alpha: (-2, 6)},
animation={"fps": 20, "time": 5},
title = (r"$\alpha$ = {:.2f}", alpha),
xlabel = ('x'), ylabel= ('y(x)'),
ylim = (-4, 4),
imodule = "panel", show = False
)
#anim.show().save("Animation.html", embed=True)
IFrame(src="Animation.html", width="100%", height="600px")
Alternativ kann man auch in Matplotlib die animierte Grafik erzeugen, speichern und im Jupyter Notebook darstellen (näheres siehe matplotlib.animation und Beispiele von Animationen).
func_alpha = lambdify((x,alpha), Loes_alpha.rhs)
from matplotlib.animation import FuncAnimation
from IPython.display import HTML
frames_N = 20
fig, ax = plt.subplots(figsize=(7, 5))
def animate(i):
ax.cla()
alpha_s = 8*i/frames_N - 2
ax.set_title(r"$\alpha$ = {:.2f}".format(alpha_s))
ax.set(xlabel='x', ylabel='y(x)', ylim=(-4, 4))
ax.plot(xvals, func_alpha(xvals, alpha_s))
return fig,
ani = FuncAnimation(fig, animate, frames=frames_N, interval=250)
ani.save("./Animation.gif", writer='pillow')
plt.close(fig)
HTML(ani.to_html5_video())
Anwendungsbeispiel: Taylorreihen Entwicklung¶
Wir möchten nun die Funktion $f(x)=\frac{1}{x}$ um die Stelle $x_0=1$ in eine Taylorreihe entwickeln. Hierzu definieren wir zunächst die Funktion:
f = 1/x
Das N-te Taylorpolynom $P_N(x)$ der Taylorreihe einer Funktion um eine Stelle $x_0$ ist wie folgt definiert:
$$ \begin{equation} P_N(x) = \sum_{k=0}^N \frac{f^{(k)}(x_0)}{k!}\cdot(x-x_0)^k \quad, \end{equation} $$
wobei $f^{(k)}(x_0)$ die k-te Ableitung der Funktion $f(x)$ an der Stelle $x=x_0$ darstellt: $f^{(k)}(x_0):=\left. \frac{d^k f(x)}{dx^k} \right|_{x=x_0}$
Das 4-te Taylorpolynom $P_4(x)$ der Taylorreihe der Funktion $f(x)=\frac{1}{x}$ um eine Stelle $x_0$ lautet z.B.:
x_0 = Symbol('x_0')
P_4 = sum(f.diff(x, k).subs(x,x_0)/factorial(k)*(x - x_0)**k for k in range(5))
P_4
Durch das Ausmultiplizieren der Klammern (Befehl "expand()") vereinfacht sich der obere Ausdruck:
P_4.expand()
Wir stellen uns im Folgenden die ersten drei Taylorpolynome $P_1(x), P_2(x)$ und $P_3(x)$ der Funktion $f(x)=\frac{1}{x}$ mit $x_0=1$ grafisch dar und vergleichen diese mit der ursprünglichen Funktion $f(x)$.
P_1 = sum(f.diff(x, k).subs(x,1)/factorial(k)*(x - 1)**k for k in range(2))
P_2 = sum(f.diff(x, k).subs(x,1)/factorial(k)*(x - 1)**k for k in range(3)).expand()
P_3 = sum(f.diff(x, k).subs(x,1)/factorial(k)*(x - 1)**k for k in range(4)).expand()
x_min = 0.4
x_max = 2.25
pf = plot(f, P_1, P_2, P_3, (x,x_min,x_max), line_color='black', legend = True, show=False)
pf[1].line_color='red'
pf[2].line_color='blue'
pf[3].line_color='green'
pf.show()
Wir stellen uns nun das 30-te Taylorpolynom $P_{30}(x)$ der Funktion $f(x)=\frac{1}{x}$ mit $x_0=1$ grafisch dar und vergleichen es wieder mit der ursprünglichen Funktion $f(x)$.
P_30 = sum(f.diff(x, k).subs(x,1)/factorial(k)*(x - 1)**k for k in range(31))
P_30.expand()
x_min = 0.05
x_max = 2.1
pf = plot(f, P_30, (x,x_min,x_max), line_color='black', show=False)
pf[1].line_color='red'
pf.show()
Anwendungsbeispiel: Herleitung der Differentiationregeln mittels der Methode der Lagrange Polynome¶
Die Ableitung (Differentiation) einer Funktion $f(x)$ stellt ein fundamentales Konzept der Mathematik dar, welches bei der Formulierung von physikalischen Bewegungsgleichungen essenziell ist. Die mathematische Definition der Ableitung einer Funktion, der sogenannte Differentialquotient, benutzt dabei eine Grenzwert-Formulierung ($h \to 0$), wobei in der numerischen Mathematik hingegen mehrere Differentiationregeln formuliert werden, die den wirklichen Wert der Ableitung approximieren. Die analytische Herleitung dieser Regeln benutzt die Methode der Lagrange Polynome (siehe Anwendungsbeispiel: Interpolation und Polynomapproximation). In diesem Anwendungsbeispiel wird die Herleitung dieser Differentiationregeln behandelt.
Die mathematische Definition der Ableitung einer Funktion $f(x)$ an der Stelle $x_0$ lautet
$$ \begin{equation} f^\prime(x_0) = \lim \limits_{h \to 0} \frac{f(x_0+h) - f(x_0)}{h} \quad. \end{equation} $$
Zur Approximation dieser Zahl werden wir nun die Funktion $f(x)$ in ein Lagrange Polynom entwickeln und diesen Ausdruck dann Differenzieren. Wir möchten jetzt die Funktion $f(x)$ mittels ($n+1$) vorgegebener Punkte ($(x_0,f(x_0)), (x_1,f(x_1)), ..., (x_n,f(x_n)) $) in ein Lagrange Polynom vom Grade $n$ entwickeln, wobei wir annehmen, dass $x_0 \in [a,b]$ und $x_n = x_0 + n \cdot h$. Der Wert von $h$ sollte hierbei hinreichend klein sein um sicherzustellen, dass $x_n \in [a,b]$ gilt.
Das n-te Lagrange Polynom $P_n(x)$ einer Funktion ist wie folgt definiert:
$$ \begin{equation} P_n(x) = \sum_{k=0}^n f(x_k) \cdot L_{n,k}(x) \quad, \end{equation} $$
mit
$$ \begin{equation} L_{n,k}(x) = \prod_{i=0 , \, i \neq k}^n \frac{(x-x_i)}{(x_k-x_i)} \quad. \end{equation} $$
Wir approximieren nun die Funktion $f(x)$ durch das n-te Lagrange Polynom $P_n(x)$
$$ \begin{equation} f(x) = P_n(x) + \underbrace{ \frac{f^{(n+1)}(\xi(x))}{(n+1)!} \prod_{i=0}^n (x-x_i) }_{\hbox{Fehler}} \,\, , \,\, \forall \,\, \xi(x) \in [x_0,x_n] \quad, \end{equation} $$ wobei $f^{(n+1)(x)}$ die $(n+1)$-te Ableitung der Funktion darstellt und $\xi(x)$ ein x-Wert im Bereich der Stützstellen ist.
Vernachlässigt man den Fehlerterm und leitet den Ausdruck ab, erhält man:
$$ \begin{equation} f^\prime(x) \approx \frac{d}{dx} P_n(x) = \frac{d}{dx} \left( \sum_{k=0}^n f(x_k) \cdot L_{n,k}(x) \right) = \sum_{k=0}^n f(x_k) \cdot \frac{d}{dx} L_{n,k}(x) = \sum_{k=0}^n f(x_k) \cdot L^\prime_{n,k}(x) \end{equation} $$
Die unterschiedlichen numerischen Differentiationregeln für $f^\prime(x)$ ergeben sich nun durch die Wahl wieviele Punkte ($n+1$) man bei der Berechnung berücksichtigt.
Wir definieren zunächst eine allgemeine SymPy-Funktion $f(x)$.
f = Function('f')(x)
Zweipunkteformel (n=1)¶
Wir betrachten nun den Fall $n=1$ und approximieren die Funktion mittels des ersten Lagrange Polynoms $P_1(x)$: $f(x) \approx P_1(x)$
x0,h = symbols('x_0,h')
points = [x0, x0+h]
f_points=[f.subs(x,points[0]), f.subs(x,points[1])]
points, f_points
In dem Unterpunkt Anwendungsbeispiel: Interpolation und Polynomapproximation der Vorlesung 5 hatten wir in dem Jupyter Notebook 'LagrangePolynome.ipynb' (View Notebook, Download Notebook) den analytischen Ausdruck des Lagrange Polynoms $P_n(x)$ mittels der folgenden Anweisungen berechnet (hier speziell $P_1(x)$):
P_1 = sum(f_points[k] * prod((x - points[i]) / (points[k] - points[i])
for i in range(len(points)) if i != k)
for k in range(len(points))
)
P_1.simplify()
Leitet man diesen Ausdruck nach $x$ ab, so erhält man $f^\prime(x_0) \approx P_1^\prime(x_0)$ .
P_1.diff(x).simplify()
Der obere Ausdruck stellt die sogenannte Zweipunkte-Differentiationsregel dar.
Zweipunkteformel (n=1): $ f^\prime(x_0) \approx \frac{1}{h} \left[ f(x_0 + h) - f(x_0) \right]$
Dreipunkte-Endpunkt-Formel (n=2)¶
Wir betrachten nun den Fall $n=2$ und approximieren die Funktion mittels des zweiten Lagrange Polynoms $P_2(x)$: $f(x) \approx P_2(x)$, wobei wir hierbei die folgenden drei Punkte verwenden: ($(x_0,f(x_0)), (x_0 + h,f(x_0 + h)), (x_0 + 2h,f(x_0 + 2h)) $).
points = [x0, x0+h, x0+2*h]
f_points=[f.subs(x,points[0]), f.subs(x,points[1]), f.subs(x,points[2])]
points, f_points
Mit den oben angegebenen zwei Punkten erhält man das folgende Lagrange Polynom $P_2(x)$:
P_2 = sum(f_points[k] * prod((x - points[i]) / (points[k] - points[i])
for i in range(len(points)) if i != k)
for k in range(len(points))
)
P_2.simplify()
Leitet man diesen Ausdruck nach $x$ ab, so erhält man $f^\prime(x) \approx P_2^\prime(x)$
P_2.diff(x)
Betrachtet man nun speziell die Ableitung bei $f^\prime(x_0) \approx P_2^\prime(x_0)$, so erhält man die Dreipunkte-Endpunkt-Formel:
$$ \begin{equation} f^\prime(x_0) \approx \frac{1}{2h} \left[ - 3\,f(x_0) + 4 \, f(x_0 + h) - f(x_0 + 2 \, h) \right] \end{equation} $$
P_2.diff(x).subs(x,x_0).simplify()
Dreipunkte-Mittelpunkt-Formel (n=2)¶
Wir betrachten nun abermals den Fall $n=2$ und approximieren die Funktion mittels des zweiten Lagrange Polynoms $P_2(x)$: $f(x) \approx P_2(x)$, wobei wir nun die folgenden drei Punkte verwenden: ($(x_0 -h,f(x_0 -h)), (x_0 ,f(x_0)), (x_0 + h,f(x_0 + h)) $).
points = [x0 -h,x0,x0+h]
f_points = [f.subs(x,points[0]), f.subs(x,points[1]), f.subs(x,points[2])]
points, f_points
P_2_M = sum(f_points[k] * prod((x - points[i]) / (points[k] - points[i])
for i in range(len(points)) if i != k)
for k in range(len(points))
)
P_2_M.simplify()
Leitet man diesen Ausdruck nach $x$ ab, so erhält man $f^\prime(x) \approx P_2^\prime(x)$
P_2_M.diff(x).expand()
Betrachtet man nun speziell die Ableitung bei $f^\prime(x_0) \approx P_2^\prime(x_0)$, so erhält man die Dreipunkte-Mittelpunkt-Formel:
$$ \begin{equation} f^\prime(x_0) \approx \frac{1}{2h} \left[ f(x_0 + h) - f(x_0 - h) \right] \end{equation} $$
P_2_M.diff(x).subs(x,x_0).simplify()
Fünfpunkte-Mittelpunkt-Formel (n=4)¶
Wir betrachten nun den Fall $n=4$ und approximieren die Funktion mittels des vierten Lagrange Polynoms $P_4(x)$: $f(x) \approx P_4(x)$, wobei wir die folgenden fünf Punkte verwenden: ($(x_0 - 2 \, h,f(x_0 - 2 \, h)), (x_0 - h,f(x_0 - h)), (x_0 ,f(x_0)), (x_0 + h,f(x_0 + h)), (x_0 + 2 \, h,f(x_0 + 2 \, h)) $).
points = [x0-2*h, x0-h, x0, x0+h, x0+2*h]
f_points = [f.subs(x,points[0]), f.subs(x,points[1]), f.subs(x,points[2]), f.subs(x,points[3]), f.subs(x,points[4])]
points, f_points
P_4_M = sum(f_points[k] * prod((x - points[i]) / (points[k] - points[i])
for i in range(len(points)) if i != k)
for k in range(len(points))
)
P_4_M.simplify()
Leitet man diesen Ausdruck nach $x$ ab, so erhält man $f^\prime(x) \approx P_4^\prime(x)$
P_4_M.diff(x).expand()
Betrachtet man nun speziell die Ableitung bei $f^\prime(x_0) \approx P_4^\prime(x_0)$, so erhält man die Fünfpunkte-Mittelpunkt-Formel:
$$ \begin{equation} f^\prime(x_0) \approx \frac{1}{12 \, h} \left[ f(x_0 - 2 \, h) - 8 \, f(x_0 - h) + 8 \, f(x_0 + h) - f(x_0 + 2 \, h) \right] \end{equation} $$
P_4_M.diff(x).subs(x,x_0).simplify()