<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Lektion 12
Step1:  Matrixexponentiale
Step2:  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}
Step3:  Fundamentalsystem
Step4:  Das Fundamentalsystem wird leider zu kompliziert
Step5:  Numerische Lösungen
Step6:  die Funktion mpmath.odefun löst die Differentialgleichung $[y_0', \dots, y_n'] = F(x, [y_0, \dots, y_n])$.
Step7:  Ergebnisse werden intern gespeichert (Cache)
Step8:  Die Pendelgleichung
Step9:  Wir lösen die AWA $y'' = -\sin(y)$, $y(0) = y_0$, $y'(0) = 0$.
Step10:  Die Energie ist eine Erhaltungsgröße.
Step11:  Lösen wir mit der Methode der Trennung der Variablen.
Step12:  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.
<ASSISTANT_TASK:>
Python Code:
from sympy import *
init_printing()
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
Explanation: Lektion 12
End of explanation
x = Symbol('x', real=True)
A = Matrix(3,3, [x,x,0,0,x,x,0,0,x])
A
A.exp()
Explanation: Matrixexponentiale
End of explanation
A = Matrix(4,4,[0,1,0,0,-1,0,1,0,0,0,0,1,1,0,-3, 0])
A
A.eigenvals()
Explanation: 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}
End of explanation
%time Phi = (x*A).exp()  # Fundamentalsystem für das System
Explanation: Fundamentalsystem
End of explanation
% time len(latex(Phi))
t = Symbol('t', real=True)
mu = list(A.eigenvals())
mu
phi = [exp(mm*t) for mm in mu]
phi
def element(i, j):
    f = phi[j]
    return f.diff(t, i)
Phi = Matrix(4, 4, element)
Phi
P1 = Phi**(-1)
len(latex(P1))
P4 = Phi.inv()
len(latex(P4))
A3 = P1*Phi
A3[0,0].n()
A4 = simplify(A3)
A4[0,0].n()
A3[0,0].simplify()
Out[44].n()
len(latex(A3[0,0]))
A2 = simplify(P1*Phi)
A2[0,0]
A2[0,0].n()
P2 = simplify(P1.expand())
len(latex(P2))
P2
(P2*Phi).simplify()
A = Out[31]
A[0,0].n()
B = Matrix([0, cos(2*t), 0, 0])
B
P2*B
P3 = Integral(P2*B, t).doit()
P3
tmp = (Phi*P3)[0]
tmp = tmp.simplify()
expand(tmp).collect([sin(2*t), cos(2*t)])
psi2 = (Phi*P3)[2]
psi2.simplify().expand()
im(psi2.simplify()).expand()
M = Matrix([0,1,t])
Integral(M, t).doit()
Explanation: Das Fundamentalsystem wird leider zu kompliziert
End of explanation
x = Symbol('x')
y = Function('y')
dgl = Eq(y(x).diff(x,2), -sin(y(x)))
dgl
#dsolve(dgl)  # NotImplementedError
Explanation: Numerische Lösungen
End of explanation
def F(x, y):
    y0, y1 = y
    w0 = y1
    w1 = -mpmath.sin(y0)
    return [w0, w1]
F(0,[0,1])
ab = [mpmath.pi/2, 0]
x0 = 0
phi = mpmath.odefun(F, x0, ab)
phi(1)
xn = np.linspace(0, 25, 200)
wn = [phi(xx)[0] for xx in xn]
dwn = [phi(xx)[1] for xx in xn]
plt.plot(xn, wn, label="$y$")
plt.plot(xn, dwn, label="$y'$")
plt.legend();
Explanation: die Funktion mpmath.odefun löst die Differentialgleichung $[y_0', \dots, y_n'] = F(x, [y_0, \dots, y_n])$.
End of explanation
%time phi(50)
%time phi(60)
%time phi(40)
Explanation: Ergebnisse werden intern gespeichert (Cache)
End of explanation
dgl
eta = Symbol('eta')
y0 = Symbol('y0')
Explanation: Die Pendelgleichung
End of explanation
H = Integral(-sin(eta), eta).doit()
H
E = y(x).diff(x)**2/2 - H.subs(eta, y(x))  # Energie
E   
E.diff(x)
E.diff(x).subs(dgl.lhs, dgl.rhs)
Explanation: Wir lösen die AWA $y'' = -\sin(y)$, $y(0) = y_0$, $y'(0) = 0$.
End of explanation
E0 = E.subs({y(x): y0, y(x).diff(x): 0})
E0
dgl_E = Eq(E, E0)
dgl_E
# dsolve(dgl_E)  # abgebrochen
Explanation: Die Energie ist eine Erhaltungsgröße.
End of explanation
Lsg = solve(dgl_E, y(x).diff(x))
Lsg
h = Lsg[0].subs(y(x), eta)
h
I1 = Integral(1/h, eta).doit()
I1
Explanation: Lösen wir mit der Methode der Trennung der Variablen.
End of explanation
I2 = Integral(1/h, (eta, y0, -y0))
I2
def T(ypsilon0):
    return 2*re(I2.subs(y0, ypsilon0).n())
T(pi/2)
phi(T(pi/2)), mpmath.pi/2
xn = np.linspace(0.1, .95*np.pi, 5)
wn = [T(yy) for yy in xn]
plt.plot(xn, wn);
Explanation: 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.
End of explanation
<END_TASK>