1.4 Sample Applications
A. Newton's equation for an oscillator
$
\frac{\textstyle d^{2}x}{\textstyle dt^{2}} = - \omega_{0}^{2} x
$
(1)
may be approximated as
$
\begin{eqnarray}
\frac{\textstyle \delta^{2} x_{n}}{\textstyle (\Delta t)^{2}} & = &
- \omega_{0}^{2} x_{n} + O[(\Delta t)^{2}]
\end{eqnarray}
$
which leads us to
$
\begin{eqnarray}
x_{n+1} = 2x_{n} - x_{n-1} + (- \omega_{0}^{2} x_{n}) (\Delta t)^{2}
+ O[(\Delta t)^{4}]
\end{eqnarray}
$
This is the famous Stoermer-Verlet algorithm for the step-by-step integration
of an equation of motion such as (1). The procedure is
applicable for any right hand side, and for any number of coupled equations of motion!
Here is a realization of this algorithm, along with several others that will be
explained later:
(1) To try out various harmonic and anharmonic oscillators, you may change the parameters
$\beta$ and $\omega_{0}^{2}$ in the text panels.
(2) To get a feeling for the limitations of the Verlet algorithm, play around with
$\Delta t$.
EXERCISE:
(see also
here)
a) Write a program to tabulate and/or draw the analytical
solution to the HO equation. (You may achieve a very concise visualization
by displaying the trajectory in phase space, i.e. in the
coordinate system
$\{ x; \dot{x} \}$; where for
$\dot{x}$ the approximation
$\dot{x} \approx (x_{n+1}-x_{n-1})/2 \Delta t$
may be used.) Choose
values of $\omega_{0}^{2}$, $\Delta t$ and
$x_{0}, \dot{x}_{0}$,
and use these to determine the exact value of
$x_{1}$. Starting
with $x_{0}$ and $x_{1}$, employ the above algorithm to compute
the further path $\{x_{n} ; n=2,3,\dots \}$. Test the performance of
your program by varying $\Delta t$ and
$\omega_{0}^{2}$.
b) Now apply your code to the anharmonic oscillator
\frac{\textstyle d^{2}x}{\textstyle dt^{2}} = - \omega_{0}^{2} x - \beta x^{3}
To start the algorithm you may use the approximate value given by
x_{1} \approx x_{0} + \dot{x}_{0} \Delta t + \ddot{x}_{0}
\frac{\textstyle (\Delta t)^{2}}{\textstyle 2}
c) The planar pendulum is described by the equation of motion
\frac{\textstyle d^{2}\phi}{\textstyle dt^{2}} =
- \frac{\textstyle g}{\textstyle l} sin \phi
Solution strategies vary between the three regions:
$\phi$ very small; $\phi$ small; $\phi$ arbitrary.
Very small $\phi$: Here $\sin \phi \approx \phi$,
and the e.o.m. is that
of a harmonic oscillator, with the usual analytic solution.
Small $\phi$: We may put $\sin \phi \approx \phi - \phi^{3}/6$,
adding an anharmonic term to the e.o.m. Again, an analytical solution
may be found,
but it is more involved than in the harmonic case; see
Landau-Lifshitz, Mechanics.
Any $\phi$:
The exact solution is known in implicit form:
t(\phi) =
\int \limits_{0}^{\phi}
\frac{\textstyle d \phi'}{\textstyle
\sqrt{\textstyle (2/ml^{2}) \left[E+mgl \cos \phi' \right]}}
One may clumsily invert this equation for regular times $t_{k}$, using
Newton-Raphson.
Instead, one may choose to simulate the pendulum, using the Stoermer-Verlet
algorithm:
\phi(t_{n+1}) =
2 \phi(t_{n})-\phi(t_{n-1})-\frac{\textstyle g}{\textstyle l} \sin \phi(t_{n})
(\Delta t)^{2}
As a starting value for $\phi(t_{1})$ one may use the Taylor approximation
\begin{eqnarray}
\phi_{1} &\approx& \phi_{0} + \dot{\phi}_{0} \Delta t + \ddot{\phi}_{0}
\frac{\textstyle (\Delta t)^{2}}{\textstyle 2}
\end{eqnarray}
B. Thermal Conduction
\begin{eqnarray}
\frac{\textstyle \partial T(x,t)}{\textstyle \partial t} & = &
\lambda \frac{\textstyle \partial^{2} T(x,t)}{\textstyle \partial x^{2}}
\end{eqnarray}
Writing $T_{i}^{n} \equiv T(x_{i},t_{n}) $,
we arrive at the "FTCS scheme" ("forward-time, centered-space")
$
\begin{eqnarray}
\frac{\textstyle 1}{\textstyle \Delta t}[T_{i}^{n+1}-T_{i}^{n}] & \approx &
\frac{\textstyle \lambda}{\textstyle (\Delta x)^{2}}
[T_{i+1}^{n}-2 T_{i}^{n}+T_{i-1}^{n}]
\end{eqnarray}
$
or
\begin{eqnarray}
T_{i}^{n+1}&=&(1-2a) T_{i}^{n} + a(T_{i-1}^{n}+T_{i+1}^{n}) \; ,
\;\;\;
with
\;\;\;
a &\equiv& \lambda \Delta t/(\Delta x)^{2}
\end{eqnarray}
for $i=1, \dots N-1$
(and $T_{0}^{n+1}$, $T_{N}^{n+1}$ given as
boundary conditions).
EXERCISE:
(see also
here)
Let us divide the rod into $N=10$ pieces of equal length, with node points
$i=0,\dots N$, and assume the boundary conditions $T(0,t) \equiv T_{0}^{n}=1.0$
and $T(L,t) \equiv T_{10}^{n}=0.5$. The values for the temperature
at time $t=0$ (the initial values) are
$T_{1}^{0}=T_{2}^{0}=$ $ \dots T_{5}^{0}= 1.0$ and
$T_{6}^{0}=T_{7}^{0}=$ $ \dots T_{9}^{0}= 0.5$ (step function).
vesely
2005-10-10