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 [ ]: