Musterlösung: Aufgabe 2, Übungsblatt 7

Einführung in die Programmierung für Studierende der Physik

(Introduction to Programming for Physicists)

Vorlesung gehalten an der J.W.Goethe-Universität in Frankfurt am Main

(Sommersemester 2022)

von Dr.phil.nat. Dr.rer.pol. Matthias Hanauske

Frankfurt am Main 28.05.2022

Analytisches Lösen von Differentialgleichungen

Der lineare harmonische Oszillator mit Dämfung

Aufgabe 2

Diese Aufgabe ist angelehnt an das Kapitel 23 "Der gedämpfte harmonische Oszillator" des Buches von Prof. Walter Greiner, Mechanik (Teil 1) [5. Auflage, 1989, siehe Seite 226- 237]. Siehe auch Vorlesungsskript von Prof. Rischke auf Seite 117- 126 http://www.th.physik.uni-frankfurt.de/~drischke/Skript_MI_WiSe2016-2017.pdf ). Wir betrachten im Folgenden den gedämpften harmonischen Oszillator am Beispiel eines reibungsfrei gelagerten Wagens (Masse=$M$) auf den eine Rückstellkraft einwirkt (die proportional zu seiner Auslenkung $x$ ist (Proportionalitätskonstante $k$)), wobei zusätzlich eine geschwindigkeitsabhängige Reibungskraft auf den Wagen einwirkt (z.B. verursacht durch den auf den Wagen einwirkende Luftwiderstand, Stokesscher Ansatz: Proportionalitätskonstante $\alpha$). Aufgrund der Rückstellkraft, besitzt das zugrundeliegende Potential $V(x)$ die Form einer Parabel $V(x)=\frac{k \, x^2}{2}$. Die Differentialgleichung des linearen harmonischen Oszillators mit Dämpfung wird mittels der folgenden Differentialgleichung zweiter Ordnung beschrieben (wir setzen $\omega_0^2=\frac{k}{M}$ und $\beta = \frac{\alpha}{2M}$): $$ \begin{equation} \ddot{x}(t) = - \omega_0^2 \, x(t) - 2 \beta \, \dot{x}(t) \end{equation} $$ Die Anfangsbedingungen seien zunächst noch allgemein gehalten: $x(0) = \alpha_1 \,\, , \,\, \dot{x}(0) = \alpha_2$. Bestimmen Sie die allgemeine Lösung der Differentialgleichung mittels eines eigenen Jupyter Notebooks. Geben Sie dann die spezielle Lösung der Differentialgleichung bei festgelegten Parameterwerten ($\omega_0^2=3$ und $\beta = 0.25$) und Anfangsbedingungen ($\alpha_1 = 0$ und $\alpha_2 = 40$) an und visualisieren Sie diese in einem x-t Diagramm. An welchem Ort befindet sich der Wagen zur Zeit $t=10$ ( $x(10)$ )?

Wir berechnen zunächst die allgemeine Lösung der Differentialgleichung (DGL):

In [1]:
from sympy import *
init_printing()

Definition der DGL:

In [2]:
t= symbols('t',real=True)
omega_0, beta = symbols('omega_0, beta', positive = True, real = True)
x = Function('x')(t)
DGL = Eq(x.diff(t).diff(t), -omega_0**2*x - 2*beta*x.diff(t))
DGL
Out[2]:
$$\frac{d^{2}}{d t^{2}} x{\left (t \right )} = - 2 \beta \frac{d}{d t} x{\left (t \right )} - \omega_{0}^{2} x{\left (t \right )}$$

Lösen der DGL (allgemein):

In [3]:
dsolve(DGL)
Out[3]:
$$x{\left (t \right )} = C_{1} e^{t \left(- \beta - \sqrt{\beta - \omega_{0}} \sqrt{\beta + \omega_{0}}\right)} + C_{2} e^{t \left(- \beta + \sqrt{\beta - \omega_{0}} \sqrt{\beta + \omega_{0}}\right)}$$

Die Konstanten $C_1$ und $C_2$ werden durch die oben angegebene Anfangsbedingung ($x(0) = \alpha_1 \,\, , \,\, \dot{x}(0) = \alpha_2$) festgelegt. Die allgemeine Lösung der DGL mit diesen Anfangsbedingung lautet dann:

In [4]:
alpha_1, alpha_2= symbols('alpha_1, alpha_2', real = True)
Loes_allg=dsolve(DGL,ics={x.subs(t,0):alpha_1,x.diff(t).subs(t, 0): alpha_2})
Loes_allg
Out[4]:
$$x{\left (t \right )} = \frac{\left(\alpha_{1} \left(\beta + \sqrt{\beta^{2} - \omega_{0}^{2}}\right) + \alpha_{2}\right) e^{t \left(- \beta + \sqrt{\beta - \omega_{0}} \sqrt{\beta + \omega_{0}}\right)}}{2 \sqrt{\beta^{2} - \omega_{0}^{2}}} - \frac{\left(\alpha_{1} \beta - \alpha_{1} \sqrt{\beta^{2} - \omega_{0}^{2}} + \alpha_{2}\right) e^{t \left(- \beta - \sqrt{\beta - \omega_{0}} \sqrt{\beta + \omega_{0}}\right)}}{2 \sqrt{\beta^{2} - \omega_{0}^{2}}}$$

Festlegung der Parameter und Anfangsbedingung auf die oben angegebenen Werte:

In [5]:
Loes_allg.rhs.subs({(omega_0,3),(beta,0.25),(alpha_1,0),(alpha_2,40)})
Out[5]:
$$6.68993608005673 i e^{t \left(-0.25 - 2.98956518577535 i\right)} - 6.68993608005673 i e^{t \left(-0.25 + 2.98956518577535 i\right)}$$

Wir sind nur an dem reelwertigen Teil der Lösung interessiert:

In [6]:
Loes_allg_FestlPara=re(Loes_allg.rhs.subs({(omega_0,3),(beta,0.25),(alpha_1,0),(alpha_2,40)}))
Loes_allg_FestlPara
Out[6]:
$$13.3798721601135 e^{- 0.25 t} \sin{\left (2.98956518577535 t \right )}$$

Wir stellen uns die berechnete analytische Lösung der DGL dar:

In [9]:
plot(Loes_allg_FestlPara,xlim=(0,10),ylim=(-12,12),ylabel="x(t)");

Zur Zeit $t=10$ befindet sich der Wagen bei $x(10)=$

In [8]:
Loes_allg_FestlPara.subs(t,10)
Out[8]:
$$-1.09688543205432$$