from IPython.display import HTML
HTML('''<script>
code_show=true;
function code_toggle() {
if (code_show){
$('div.input').hide();
} else {
$('div.input').show();
}
code_show = !code_show
}
$( document ).ready(code_toggle);
</script>
<form action="javascript:code_toggle()"><input type="submit" value="Click here to toggle on/off the raw code."></form>''')
Bestimmen Sie für die Differentialgleichung
\begin{equation*} \ddot{q} + \omega_0^2 (q+\beta^2q^2) = P \cos(\Omega t) \end{equation*}1) eine harmonische Lösung unter Nutzung der Störungsrechnung. Wie verhalten sich die Amplituden im Resonanzfall bei Abbruch nach dem ersten Glied bzw. nach dem zweiten Glied?
2) eine subharmonische Lösung mit der Periode $T = 4\pi/\Omega$
Überführen der DGL in die Form
\begin{align*} \eta^2 x^{\prime \prime} + x + \mu (x^2 - P_0 \cos(\tau)) = 0 \end{align*}Potenzreihen:
\begin{align*} x(\tau) &= x_0(\tau) + \mu x_1(\tau) + \mu^2 x_2(\tau) \\ \Omega(\tau) &= \Omega_0(\tau) + \mu \Omega_1(\tau) + \mu^2 \Omega_2(\tau) \\ \text{bzw. }\; \eta(\tau) &= 1 + \mu \eta_1(\tau) + \mu^2 \eta_2(\tau) \\ \end{align*}Aus den Lösungen der Koeffizientenvergleiche in $\mu^0, \mu^1, \mu^2$ unter der Forderung nach periodischen Lösungen und einer Anfangsbedingung mit $A_1 = 0$ folgt
\begin{align*} \eta = 1 - \mu \frac{P_0}{2 A_0} - \mu^2 \left( \frac{5}{12} A_0^2 + \frac{P_0^2}{8A_0^2} \right) + \mathcal{O}(\mu^3) \end{align*}import matplotlib.pyplot as plt
import numpy as np
A = np.hstack((np.linspace(-1.0, -0.001, 100), np.NaN, np.linspace(0.001, 1.0, 100)))
A_abs = np.abs(A)
mu = 0.2
P0 = 2E-4
# Abbruch nach dem 1. Glied
eta1 = 1.0-mu*P0/2/A
# Abbruch nach dem 2. Glied
eta2 = 1.0-mu*P0/2/A-mu**2*(5/12*A**2+P0**2/8/A**2)
plt.plot(eta1, A_abs)
plt.plot(eta2, A_abs)
plt.xlabel("η")
plt.ylabel("|A|")
plt.legend(('Abbruch nach dem 1. Glied', 'Abbruch nach dem 2. Glied'))
plt.xlim(0.98, 1.02)
plt.show()
from sympy import *
init_printing()
## variables
t = symbols('t', real=True)
## functions
q = Function('q')(t)
## constants
om, Om, beta, P = symbols('omega_0 Omega beta P', real=True)
q0, A, B = symbols('q_0 A B', real=True)
dgl = Eq(Derivative(q, t, t) + om**2*(q+beta**2*q**2) - P*cos(Om*t),0)
dgl
einsetzen:
dgl = dgl.replace(q,q0+A*cos(Om/2*t)+B*sin(Om/2*t))
dgl
Ableitung ausführen und Terme zusammenfassen:
dgl = dgl.expand().doit().simplify().trigsimp()
# sin^2(x) bzw. cos^2(x) ersetzen ...
dgl = dgl.rewrite(sin, exp).rewrite(cos, exp).expand().rewrite(exp, sin).simplify()
dgl
KV_skal = Eq(dgl.lhs.as_independent(t)[0],dgl.rhs.as_independent(t)[0])
KV_skal
KV_cos1_2 = Eq(dgl.lhs.coeff(cos(Om/2*t)),dgl.rhs.coeff(cos(Om/2*t)))
KV_cos1_2
KV_cos = Eq(dgl.lhs.coeff(cos(Om*t)),dgl.rhs.coeff(cos(Om*t)))
KV_cos
KV_sin1_2 = Eq(dgl.lhs.coeff(sin(Om/2*t)),dgl.rhs.coeff(sin(Om/2*t)))
KV_sin1_2
KV_sin = Eq(dgl.lhs.coeff(sin(Om*t)),dgl.rhs.coeff(sin(Om*t)))
KV_sin
Aus dem Koeffizientenvergleich in $sin(\Omega t)$ folgt $A = 0 \lor B=0$. In Kombination mit dem Koeffizientenvergleich in $\cos(\Omega t)$ mit $\beta^2 > 0$, $\omega_0^2>0$, $A^2>0$, $B^2>0$ folgt $B = 0$.
Die gesuchte Lösung $A(P,\Omega)$ aus dem Koeffizientenvergleich in $\cos(\Omega t)$ :
Ares = solve(KV_cos.replace(B,0),A)
Ares[0] = Ares[0].replace(om,Om/2)
Ares[1] = Ares[1].replace(om,Om/2)
Eq(A,Ares[1].replace(om,Om/2))
$q_0$ kann anschließend aus dem skalaren Koeffizientenvergleich ermittelt werden:
q0res = solve(KV_skal.replace(B,0).replace(A,Ares[1]),q0)
q0res
$q_0 = -\frac{1}{2\beta^2} \pm \sqrt{\frac{1}{4\beta^4}-\frac{4P}{\Omega^2\beta^4}}$
Zuletzt folgt aus dem Koeffizientenvergleich in $\cos(\Omega/2 t)$:
KV_cos1_2_1 = KV_cos1_2.replace(A,Ares[1].replace(om,Om/2)).replace(q0,q0res[0]).replace(om,Om/2).expand().doit().simplify()
KV_cos1_2_1
Mit $P \neq 0$, $\beta \neq 0$ folgt eine implizite Darstellung für $\Omega$.
KV_cos1_2_2 = Eq(KV_cos1_2_1.lhs*2*beta/sqrt(2)/sqrt(P)+Om,KV_cos1_2_1.rhs*2*beta/sqrt(2)/sqrt(P)+Om)
KV_cos1_2_2
Die Lösung existiert nur für $16 P \beta^2 < \Omega^2$ bzw. $4 P \beta^2 < \omega_0^2$