Franz J. Vesely > CompPhys Tutorial > Finite Differences  

< >
Ch. 1 Sec. 4

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:

(An-) Harmonic Oscillator:
Start Applet

(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$.

(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} $
\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).

(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

< >