Lektion 12


In [1]:
from sympy import *
init_printing()
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

Matrixexponentiale


In [2]:
x = Symbol('x', real=True)

In [3]:
A = Matrix(3,3, [x,x,0,0,x,x,0,0,x])
A


Out[3]:
$$\left[\begin{matrix}x & x & 0\\0 & x & x\\0 & 0 & x\end{matrix}\right]$$

In [4]:
A.exp()


Out[4]:
$$\left[\begin{matrix}e^{x} & x e^{x} & \frac{x^{2} e^{x}}{2}\\0 & e^{x} & x e^{x}\\0 & 0 & e^{x}\end{matrix}\right]$$

Koeffizienten


In [5]:
x = Symbol('x')
p = series(exp(x+2*x**2)**2, x, 0, 10)
p


Out[5]:
$$1 + 2 x + 6 x^{2} + \frac{28 x^{3}}{3} + \frac{50 x^{4}}{3} + \frac{108 x^{5}}{5} + \frac{1324 x^{6}}{45} + \frac{10424 x^{7}}{315} + \frac{3958 x^{8}}{105} + \frac{21428 x^{9}}{567} + \mathcal{O}\left(x^{10}\right)$$

In [9]:
p.coeff(x, 0)


Out[9]:
$$1$$

In [10]:
a = ((1-exp(x)+exp(2*x))**5).expand()
a


Out[10]:
$$e^{10 x} - 5 e^{9 x} + 15 e^{8 x} - 30 e^{7 x} + 45 e^{6 x} - 51 e^{5 x} + 45 e^{4 x} - 30 e^{3 x} + 15 e^{2 x} - 5 e^{x} + 1$$

In [11]:
a.coeff(exp(6*x))


Out[11]:
$$45$$

In [12]:
b = (1+x**2)**3
b


Out[12]:
$$\left(x^{2} + 1\right)^{3}$$

In [13]:
b.coeff(x**2)


Out[13]:
$$0$$

In [14]:
b.expand().coeff(x**2)


Out[14]:
$$3$$

Böse Falle!

Gekoppelte Pendel

\begin{align} y'' &= w - y + \cos(2t)\\ w'' &= y - 3w \end{align}

Übersetzt sich in \begin{align*} y_0' &= y_1 \\ y_1' &= y_2 - y_0 + \cos(2t) \\ y_2' &= y_3 \\ y_3' &= y_0 - 3 y_2 \end{align*}


In [15]:
y = Function('y')
w = Function('w')
t = Symbol('t', real=True)

In [16]:
dgl = {}
dgl[0] = Eq(y(t).diff(t,2), w(t) - y(t) + cos(2*t))
dgl[1] = Eq(w(t).diff(t,2), y(t) - 3*w(t))

In [17]:
A = Matrix(4,4,[0,1,0,0,-1,0,1,0,0,0,0,1,1,0,-3, 0])
A


Out[17]:
$$\left[\begin{matrix}0 & 1 & 0 & 0\\-1 & 0 & 1 & 0\\0 & 0 & 0 & 1\\1 & 0 & -3 & 0\end{matrix}\right]$$

In [18]:
A.eigenvals()


Out[18]:
$$\left \{ - i \sqrt{- \sqrt{2} + 2} : 1, \quad i \sqrt{- \sqrt{2} + 2} : 1, \quad - i \sqrt{\sqrt{2} + 2} : 1, \quad i \sqrt{\sqrt{2} + 2} : 1\right \}$$

Fundamentalsystem


In [19]:
%time Phi = (x*A).exp()  # Fundamentalsystem für das System


CPU times: user 2.22 s, sys: 0 ns, total: 2.22 s
Wall time: 2.22 s

Das Fundamentalsystem wird leider zu kompliziert


In [20]:
% time len(latex(Phi))


CPU times: user 30.3 s, sys: 8 ms, total: 30.3 s
Wall time: 30.3 s
Out[20]:
$$123417$$

In [21]:
a = [Symbol('a_'+str(j), real=True) for j in range(2)]
b = [Symbol('b_'+str(j), real=True) for j in range(2)]

Ansatz für eine spezielle Lösung des inhomogenen Systems


In [22]:
phi = a[0]*cos(2*t) + a[1]*sin(2*t)
psi = b[0]*cos(2*t) + b[1]*sin(2*t)

In [23]:
ers = {}
ers[y(t)] = phi
ers[w(t)] = psi

In [26]:
glg = [dgl[j].subs(ers).doit() for j in range(2)]
glg


Out[26]:
$$\left [ - 4 \left(a_{0} \cos{\left (2 t \right )} + a_{1} \sin{\left (2 t \right )}\right) = - a_{0} \cos{\left (2 t \right )} - a_{1} \sin{\left (2 t \right )} + b_{0} \cos{\left (2 t \right )} + b_{1} \sin{\left (2 t \right )} + \cos{\left (2 t \right )}, \quad - 4 \left(b_{0} \cos{\left (2 t \right )} + b_{1} \sin{\left (2 t \right )}\right) = a_{0} \cos{\left (2 t \right )} + a_{1} \sin{\left (2 t \right )} - 3 b_{0} \cos{\left (2 t \right )} - 3 b_{1} \sin{\left (2 t \right )}\right ]$$

In [27]:
gls = []
for gl in glg:
    diff = (gl.rhs - gl.lhs).expand()
    gls.append(diff.coeff(cos(2*t)))
    gls.append(diff.coeff(sin(2*t)))

gls


Out[27]:
$$\left [ 3 a_{0} + b_{0} + 1, \quad 3 a_{1} + b_{1}, \quad a_{0} + b_{0}, \quad a_{1} + b_{1}\right ]$$

In [28]:
Lsg = solve(gls)
Lsg


Out[28]:
$$\left \{ a_{0} : - \frac{1}{2}, \quad a_{1} : 0, \quad b_{0} : \frac{1}{2}, \quad b_{1} : 0\right \}$$

Probe:


In [29]:
yp = phi.subs(Lsg)
yp


Out[29]:
$$- \frac{1}{2} \cos{\left (2 t \right )}$$

In [30]:
wp = psi.subs(Lsg)
wp


Out[30]:
$$\frac{1}{2} \cos{\left (2 t \right )}$$

In [31]:
ers_probe = {}
ers_probe[y(t)] = yp
ers_probe[w(t)] = wp

In [32]:
[dgl[j].subs(ers_probe).doit() for j in range(2)]


Out[32]:
$$\left [ \mathrm{True}, \quad \mathrm{True}\right ]$$

Numerische Lösungen


In [33]:
x = Symbol('x')
y = Function('y')

In [34]:
dgl = Eq(y(x).diff(x,2), -sin(y(x)))
dgl


Out[34]:
$$\frac{d^{2}}{d x^{2}} y{\left (x \right )} = - \sin{\left (y{\left (x \right )} \right )}$$

In [36]:
#dsolve(dgl)  # NotImplementedError

die Funktion mpmath.odefun löst die Differentialgleichung $[y_0', \dots, y_n'] = F(x, [y_0, \dots, y_n])$.


In [37]:
def F(x, y):
    y0, y1 = y
    w0 = y1
    w1 = -mpmath.sin(y0)
    return [w0, w1]

In [38]:
F(0,[0,1])


Out[38]:
[1, mpf('0.0')]

In [39]:
ab = [mpmath.pi/2, 0]
x0 = 0

In [40]:
phi = mpmath.odefun(F, x0, ab)
phi(1)


Out[40]:
[mpf('1.0749116843722417'), mpf('-0.97551004396953367')]

In [41]:
xn = np.linspace(0, 25, 200)
wn = [phi(xx)[0] for xx in xn]
dwn = [phi(xx)[1] for xx in xn]

In [42]:
plt.plot(xn, wn, label="$y$")
plt.plot(xn, dwn, label="$y'$")
plt.legend();


Ergebnisse werden intern gespeichert (Cache)


In [43]:
%time phi(50)


CPU times: user 1.98 s, sys: 0 ns, total: 1.98 s
Wall time: 1.98 s
Out[43]:
[mpf('-0.084824921581113169'), mpf('1.4116688868027041')]

In [44]:
%time phi(60)


CPU times: user 772 ms, sys: 0 ns, total: 772 ms
Wall time: 770 ms
Out[44]:
[mpf('1.3469819877501543'), mpf('-0.6662588540047103')]

In [45]:
%time phi(40)


CPU times: user 0 ns, sys: 0 ns, total: 0 ns
Wall time: 428 µs
Out[45]:
[mpf('-1.2600336115087054'), mpf('-0.78202929230655394')]

Die Pendelgleichung


In [46]:
dgl


Out[46]:
$$\frac{d^{2}}{d x^{2}} y{\left (x \right )} = - \sin{\left (y{\left (x \right )} \right )}$$

In [47]:
eta = Symbol('eta')
y0 = Symbol('y0')

Wir lösen die AWA $y'' = -\sin(y)$, $y(0) = y_0$, $y'(0) = 0$.


In [48]:
H = Integral(-sin(eta), eta).doit()
H


Out[48]:
$$\cos{\left (\eta \right )}$$

In [52]:
E = y(x).diff(x)**2/2 - H.subs(eta, y(x))  # Energie
E


Out[52]:
$$- \cos{\left (y{\left (x \right )} \right )} + \frac{1}{2} \frac{d}{d x} y{\left (x \right )}^{2}$$

In [53]:
E.diff(x)


Out[53]:
$$\sin{\left (y{\left (x \right )} \right )} \frac{d}{d x} y{\left (x \right )} + \frac{d}{d x} y{\left (x \right )} \frac{d^{2}}{d x^{2}} y{\left (x \right )}$$

In [54]:
E.diff(x).subs(dgl.lhs, dgl.rhs)


Out[54]:
$$0$$

Die Energie ist eine Erhaltungsgröße.


In [55]:
E0 = E.subs({y(x): y0, y(x).diff(x): 0})
E0


Out[55]:
$$- \cos{\left (y_{0} \right )}$$

In [56]:
dgl_E = Eq(E, E0)
dgl_E


Out[56]:
$$- \cos{\left (y{\left (x \right )} \right )} + \frac{1}{2} \frac{d}{d x} y{\left (x \right )}^{2} = - \cos{\left (y_{0} \right )}$$

In [ ]:
# dsolve(dgl_E)  # abgebrochen

Lösen wir mit der Methode der Trennung der Variablen.


In [57]:
Lsg = solve(dgl_E, y(x).diff(x))
Lsg


Out[57]:
$$\left [ - \sqrt{- 2 \cos{\left (y_{0} \right )} + 2 \cos{\left (y{\left (x \right )} \right )}}, \quad \sqrt{- 2 \cos{\left (y_{0} \right )} + 2 \cos{\left (y{\left (x \right )} \right )}}\right ]$$

In [59]:
print(dgl_E)


Eq(-cos(y(x)) + Derivative(y(x), x)**2/2, -cos(y0))

In [60]:
h = Lsg[0].subs(y(x), eta)
h


Out[60]:
$$- \sqrt{2 \cos{\left (\eta \right )} - 2 \cos{\left (y_{0} \right )}}$$

In [62]:
I1 = Integral(1/h, eta).doit()
I1


Out[62]:
$$- \frac{\sqrt{2}}{2} \int \frac{1}{\sqrt{\cos{\left (\eta \right )} - \cos{\left (y_{0} \right )}}}\, d\eta$$

In der Tat nicht elementar integrierbar.

Trennung der Variablem führt zu $$ -\frac{\sqrt2}2 \int_{y_0}^{y(x)} \frac{d\eta}{\sqrt{\cos(\eta)-\cos(y_0)}} = x. $$ Insbesondere ist $$ -\frac{\sqrt2}2 \int_{y_0}^{-y_0} \frac{d\eta}{\sqrt{\cos(\eta)-\cos(y_0)}} $$ gleich der halben Schwingungsperiode.


In [63]:
I2 = Integral(1/h, (eta, y0, -y0))
I2


Out[63]:
$$\int_{y_{0}}^{- y_{0}} - \frac{1}{\sqrt{2 \cos{\left (\eta \right )} - 2 \cos{\left (y_{0} \right )}}}\, d\eta$$

In [64]:
def T(ypsilon0):
    return 2*re(I2.subs(y0, ypsilon0).n())

In [65]:
T(pi/2)


Out[65]:
$$7.41629870920549$$

In [66]:
phi(T(pi/2)), mpmath.pi/2


Out[66]:
([mpf('1.5707963267948966'), mpf('6.5086722915372182e-17')],
 mpf('1.5707963267948966'))

In [67]:
xn = np.linspace(0.1, .95*np.pi, 5)
wn = [T(yy) for yy in xn]

In [68]:
plt.plot(xn, wn);



In [ ]: