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.
Eq(y.subs(x,0),3)
Eq(y.diff(x).subs(x, 0),2)
Loes = dsolve(DGL,ics={y.subs(x,0):3,y.diff(x).subs(x, 0): 2})
Loes
Die Lösung können wir dann wie folgt visualisieren:
plot(Loes.rhs)
Die Lösung können wir uns auch mittels "matplotlib" visualisieren:
import matplotlib.pyplot as plt
import numpy as np
xvals = np.linspace(-10, 10, 300)
func = lambdify(x, Loes.rhs)
plt.ylabel('y(x)')
plt.xlabel('x')
plt.plot(xvals,func(xvals));
func(xvals[:10])
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')
sum(f.diff(x, k).subs(x,x_0)/factorial(k)*(x - x_0)**k for k in range(5))
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))
P_3=sum(f.diff(x, k).subs(x,1)/factorial(k)*(x - 1)**k for k in range(4))
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
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 der 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} $$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
Mit den oben angegebenen zwei Punkten erhält man das folgende Lagrange Polynom $P_1(x)$:
List_ranges=[]
for i in range(len(points)):
List_range=list(range(len(points)))
List_range.pop(i)
List_ranges.append(List_range)
List_ranges
P_1=sum(f_points[k] * prod((x - points[i])/(points[k] - points[i]) for i in List_ranges[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.
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)$:
List_ranges=[]
for i in range(len(points)):
List_range=list(range(len(points)))
List_range.pop(i)
List_ranges.append(List_range)
List_ranges
P_2=sum(f_points[k] * prod((x - points[i])/(points[k] - points[i]) for i in List_ranges[k]) for k in range(len(points)))
P_2.expand()
Leitet man diesen Ausdruck nach $x$ ab, so erhält man $f^\prime(x) \approx P_2^\prime(x)$
P_2.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-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
List_ranges=[]
for i in range(len(points)):
List_range=list(range(len(points)))
List_range.pop(i)
List_ranges.append(List_range)
List_ranges
P_2_M=sum(f_points[k] * prod((x - points[i])/(points[k] - points[i]) for i in List_ranges[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
List_ranges=[]
for i in range(len(points)):
List_range=list(range(len(points)))
List_range.pop(i)
List_ranges.append(List_range)
List_ranges
P_4_M=sum(f_points[k] * prod((x - points[i])/(points[k] - points[i]) for i in List_ranges[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()