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