\documentclass[11pt]{article}
%\usepackage[latin1]{inputenc}
\usepackage{a4wide}
\usepackage{german}
\usepackage{epsfig}

\topmargin=-1.5cm
\textheight=26cm

\pagestyle{empty}

\begin{document}
\centerline{\it  Physik auf dem Computer I}
\bigskip
\centerline{\bf \"Ubung 7 Runge-Kutta-Verfahren, Chaos}
\medskip

Bitte protokollieren Sie Ihre Arbeit an diesem \"Ubungsblatt durch
Beantwortung der Fragen {\it in der Aufgabe 2}. Vergessen Sie nicht, Ihre(n)
Namen auf dem Protokoll zu vermerken und {\it geben Sie dies am Ende der
\"Ubung Ihrem Assistenten.} Sollten Sie an keiner der
\"Ubungsgruppen teilnehmen k\"onnen, so senden Sie Ihre L\"osungen 
bis sp\"atestens zur n\"achsten Vorlesung per e-mail (reines ASCII Format,
bitte keine attachments).

\bigskip

{\bf Aufgabe 1.\ } {\it Runge-Kutta f\"ur Pendel} \\
Programmieren Sie das chaotische getriebene Pendel $(g/l = 1)$ 
im Gravitationsfeld mit D\"ampfung $(q=2\gamma),$ das durch 
die Differentialgleichung 2. Ordnung
\begin{equation} \label{e:theta}
\ddot{\theta} = - \sin (\theta) -
  2 \gamma \dot{\theta} + \delta \sin (\omega_d t)
\end{equation}
beschrieben wird.
Implementieren Sie diese Gleichung numerisch als gekoppeltes Gleichungssystem
1. Ordnung:
\begin{eqnarray} 
\dot{\omega} &=& - \sin (\theta) -
  2\gamma \dot{\theta} + \delta \sin (\omega_d t), \\
\dot{\theta} &=& \omega.
\end{eqnarray}
Da das Euler-Verfahren instabil ist und der Verlet-Algorithmus nur
f\"ur geschwindigkeitsunabh\"angige Kr\"afte geeignet ist,
verwenden Sie diesmal bitte das Runge-Kutta-Verfahren vierter Ordnung,
welches die Auswertung der Beschleunigung an drei St\"utzstellen
vornimmt (zweite Ordnung siehe Vorlesung). 
F\"ur eine Differentialgleichung erster Ordnung,
\begin{equation}
 \dot{x}  =  f(x,\dot{x},t),
\end{equation}
lautet der Runge-Kutta-Algorithmus mit Zeitschritt $\Delta t,$ 
\begin{eqnarray*}
  m_0 &\equiv & f(x, t),\\
  m_1 &\equiv & f(x+(1/2)\Delta t\, m_0, t+(1/2)\Delta t), \\
  m_2 &\equiv & f(x+(1/2)\Delta t\, m_1, t+(1/2)\Delta t), \\
  m_3 &\equiv & f(x+\Delta t\, m_2, t+ \Delta t), \\
  x(t+\Delta t) & = & x + (1/6)\Delta t (m_0 + 2m_1 +
                                2 m_2 + m_3). 
\end{eqnarray*}
Diese Vorschrift ist in der folgenden Abbildung noch einmal grafisch
veranschaulicht: 

\medskip
\centerline{\epsfig{file=runge-kutta.eps,height=0.3\textheight}}

\newpage

F\"ur das Problem des Pendels  
kann man jetzt $x$ als den Vektor $(\theta,\omega)$ und $f(x,t)$ als 
den Vektor $(\omega, -\sin(\theta)-2 \gamma \dot{\theta} + \delta \sin
(\omega_d t) )$ auffassen. 
Das angegebene Gleichungssystem wird dann zu einem Schema 
aus vektoriellen Gleichungen, von
denen hier nur die Gleichungen f\"ur die Komponenten der ersten beiden
aufgeschrieben seien:
\begin{eqnarray*}
m^{(\theta)}_0 & = & \omega,\\
m^{(\omega)}_0 & = & -\sin(\theta) -
  2 \gamma \omega + \delta \sin (\omega_d t); \\
%
m^{(\theta)}_1 & = & \omega + (1/2)\Delta t m^{(\omega)}_0,\\
m^{(\omega)}_1 & = & -\sin(\theta + (1/2)\Delta t m^{(\theta)}_0 ) -
  2 \gamma (\omega + (1/2)m^{(\omega)}_0)  + \delta \sin (\omega_d
  (t+(1/2)\Delta t) ); \\
  \dots & & 
\end{eqnarray*}

Ein f\"ur die Implementierung geeignetes Programmger\"ust ist in\\ 
\verb|/user0/sts/COMPUTERPHYSIK/BLATT7/pendulum.cc| \\
vorgegeben. 
Ersetzen Sie darin die Funktion \verb|accel| durch eine, die den
Beschleunigungsterm der Gleichung (\ref{e:theta}) berechnet. 
Die eigentliche Integration soll in der Funktion \verb|runge_kutta_4|
durchgef\"uhrt werden, die bisher ein einfaches Euler-Verfahren
implementiert. 

Stellen Sie sicher, da\ss{} nach der Integration die Winkelkoordinate
in das Intervall $[-\pi\dots\pi]$ verlegt wird. Hierzu kann ein Test
dienen, ob dieser Bereich verlassen worden ist, und dann kann die
Einhaltung des Wertebereiches entsprechend durch Addition bzw.\ 
Subtraktion von \verb|2*M_PI| sichergestellt werden.
  
{\it Sollten Sie diese Aufgabe bis zur H\"alfte der
verf\"ugbaren Zeit nicht geschafft haben, dann bearbeiten Sie die
folgende Aufgabe bitte mit dem in \verb|pendulum_rk4.cc| gegebenen
L\"osungsvorschlag!} 

\bigskip

{\bf Aufgabe 2.\ } {\it Chaos} \\
Setzen Sie $\gamma=0.25, \omega_d=2/3, {\rm d}t = 0.1.$
Die Anfangsbedingungen seien $\theta(0)=0, \omega(0)=1.$ 
Legen sie dreispaltige Dateien mit jeweils einer Spalte f\"ur 
$t,\ \theta,$ und $\omega$ an.\\
a) Plotten Sie die Position $\theta$ sowie die Winkelgeschwindigkeit
$\omega$ gegen die Zeit $t$ f\"ur $\alpha = 0,0.5,1.2,2.2.$
Ermitteln Sie in den F\"allen, in denen das m\"oglich ist, die 
Periode der Oszillation und vergleichen Sie mit $\omega_d.$
Entspricht das Ihrer Erwartung von den Resultaten f\"ur den angetriebenen 
{\it harmonischen} Oszillator?
\\
b) Plotten Sie das Phasenraumportr\"at des Oszillators, d.h.
$\omega$ gegen $\theta$, f\"ur $\alpha=0,0.5,1.2,2.2.$ 


\bigskip

{\bf Aufgabe 3.\ } {\it freiwillig, interaktive grafische Darstellung}\\
Wir wollen die Bewegung des Pendels interaktiv mit Hilfe der 
X11-Fensteroberfl\"ache darstellen. Ein Beispielprogramm zu einer 
Animation eines ged\"ampften Oszillators ist dazu in\\ 
\verb|/user0/sts/COMPUTERPHYSIK/BLATT7/grafik.cc| \\
gezeigt. Kopieren Sie \verb|grafik.cc| sowie 
die ``Projektdatei'' \verb|Makefile|
in Ihr lokales Verzeichnis und kompilieren Sie dort mit \verb|make|.

Zur Visualisierung des getriebenen Pendels kann etwa eine Linie vom
Ursprung $(0,0)$ bis zu einer Position auf dem Einheitskreis 
$(\sin \theta, -\cos \theta)$ gezeichnet werden, die der
momentanen Auslenkung entspricht. Die dazu n\"otigen Programmelemente
sind im Beispiel vorhanden. 
\bigskip

Viel Spa\ss{}! 

\vfill
\begin{flushright}
\begin{tabular}{lll}
Michael Brunner & Di 9:45 & {\tt michi@ica1.uni-stuttgart.de} \\
Ansgar Hennecke & Mo 11:30 & {\tt ansgar@theo2.physik.uni-stuttgart.de} \\
Stefan Schwarzer & Di 11:30 & {\tt sts@ica1.uni-stuttgart.de} \\ 
\end{tabular}
\end{flushright}

\end{document}








