Homework 08 - Harmonic and Anharmonic Motion

ComputationalPhysics_HW 岳绍圣2013301020033

---

CONTENT:

• Homework 08 - Harmonic and Anharmonic Motion
  • Abstract
  • Background
  • 1. Simple Pendulum
    • 1.1 Euler-Cromer Method
    • 1.1 Period Independence of Amplitude
    • 1.2 Conclusion
  • 2. General Oscillator
    • 2.1 Analytical Analysis
    • 2.2 Numerical Analysis
      • 2.2.1Comparasion between several methods:
      • 2.2.2 Derivation of RK4 for 2^{nd} order ODE[5]
      • 2.2.3 Solve \frac{d^2\theta}{dt^2} = - \theta^3
    • 2.3 Conslusion
  • Reference:

Abstract

In this homework several problems involving the pendulum are considered, several methods of solving ODE including
Euler method, Euler-Cromer method, Second-order Runge-Kutta method, Centered Difference Method are discussed. Also, a general problem in form $d^2x/dt^2=-kx^\alpha$ is discussed, and draw the conclusion. Both analytical and numerical solution are listed, and we finally draw the conclusion that Euler-Cromer method is the most suitable way to investigate oscillation problems and while $\alpha\neq1$ the oscillation period is related to the amplitude.

Key Words: Simple pendulum; Anharmonic Pendulum; Euler Method; Euler-Cromer Method; Runge-Kutta Method; Centered Difference Method.

Background

Many of problems encountered in physics involves ordinary differential eqautions for which some initial values are specified. For example, the motion of a projectile is described by second-order differential equations that involve the position as a function of time. There are many numerical methods to solve these ordinary equation, such as Euler method, Euler-Cromer method, Second-order Runge-Kutta method, Centered Difference Method and so on. All of these method are approximations of "the mean value theorem": there is a value $t_m$ in the interval $[t,t+\Delta t]$ such that the exact solution can be gotten while stopping at the first order in $\Delta t$

$$x(t+\Delta t)=x(t)+\frac{dx}{dt}|_{t_m}\Delta t$$

All of these method are finding a proper way to approximate the slope $\frac{dx}{dt}|_{t_m}$.

Euler method uses the slope at $t$ as an approximation of the slope at $t_m$, which as compensation wil cause asymmetric effect. Part of the problem with Euler method is the one-side approximation of the slope $dx/dt|_{t_m}$ appearing in the mean value theorem. This problem can be reduced by using a more symmetric evaluation of $dy/dt$ and $dv/dt$, namely, evaluating $dv/dt$ at $t$ while evaluating $dy/dt$ at $t+\Delta t$. And this slight change yields the Euler-Cromer method. Runge-Kutta choose a more elaborate way to balance the asymmetry, for more details please see this book.

1. Simple Pendulum

1.1 Euler-Cromer Method

For simple harmonic motion, that is, linear, undamped, no driving.

$$\frac{{{d^2}\theta }}{{d{t^2}}} = - \frac{g}{l}\theta$$

It's easy to verify that it has the general solution:

$$\theta = {\theta _0}\sin \left( {\Omega t + \phi } \right)\quad\quad \Omega=\sqrt{g/l}$$

Below is a comparison among Euler method, Euler-Cromer method and Analytical solution of angle vs time, energy vs time.
点击图片查看源代码

1.1 Period Independence of Amplitude

点击图片查看源代码

1.2 Conclusion

As shown above, the asymmetry of Euler method will case energy increasing, which is contradictory with energy conservation theorem. The more symmetric method Euler-Cromer method solved this contradictory more perfectly.

For simple pendulum, the period is independent of amplitude.

2. General Oscillator

2.1 Analytical Analysis

Since the general form for equation of motion of the anharmonic oscillator is given by:

$$\frac{d^2\theta}{dt^2} = - k\theta^\alpha$$

Although with the fact that most of ODEs have no analytical solution, we should expect explicit analytical solution form for such not so complex ODE above. And as a matter of fact, it do have analytical solution.
The result[1] given by Lev Landau is:

$$\tau=\frac{2\sqrt{2m}E^{1/n-1/2}}{A^{1/n}}\int_{0}^{1}\frac{dy}{\sqrt{1-y^n}}$$

Where $E$ is the total energy of the system, in such case:

$$\left\{\begin{matrix} A=\frac{mk}{\alpha+1}\\ n=\alpha+1\\ E=\frac{mk}{\alpha+1}\theta_0^n \end{matrix}\right.$$

Employ Mathematica [2]:

$$\tau=\theta _0^{\frac{1-\alpha}{2}} \sqrt{\frac{2 (\alpha +1)}{k}} \left(\frac{\sqrt{\pi } \Gamma \left(1+\frac{1}{\alpha +1}\right)}{\Gamma \left(\frac{1}{2}+\frac{1}{\alpha +1}\right)}+\, _2F_1\left(\frac{1}{2},\frac{1}{\alpha +1};1+\frac{1}{\alpha +1};e^{i \pi (\alpha +1)}\right)\right)$$

Condition	Period
$\alpha\leq 0$	$\infty$
$\alpha=1$	$2\pi\sqrt{\frac{1}{k}}$
$\alpha=2$	$\frac{1}{\sqrt{\theta _0}} \sqrt{\frac{6}{k}} \left(\, _2F_1\left(\frac{1}{3},\frac{1}{2};\frac{4}{3};-1\right)+\frac{\sqrt{\pi } \Gamma \left(\frac{4}{3}\right)}{\Gamma\left(\frac{5}{6}\right)}\right)$
$\alpha=3$	$\frac{1}{\theta _0}\frac{4 \sqrt{2 \pi } \sqrt{\frac{1}{k}} \Gamma \left(\frac{5}{4}\right)}{ \Gamma \left(\frac{3}{4}\right)}$
$\alpha=4$	$\frac{1}{\theta _0^{3/2}}\sqrt{\frac{10}{k}} \left(\, _2F_1\left(\frac{1}{5},\frac{1}{2};\frac{6}{5};-1\right)+\frac{\sqrt{\pi } \Gamma \left(\frac{6}{5}\right)}{\Gamma \left(\frac{7}{10}\right)}\right)$
$\alpha=5$	$\frac{1}{\theta _0^2}\frac{4 \sqrt{3 \pi } \sqrt{\frac{1}{k}} \Gamma \left(\frac{7}{6}\right)}{\Gamma \left(\frac{2}{3}\right)}$

2.2 Numerical Analysis

2.2.1Comparasion between several methods:

Local truncation errors[3] is the error caused by one iteration. Global truncation error[4] is the error at a fixed time t, after however many steps the methods needs to take to reach that time from the initial time. The global truncation error is the cumulative effect of the local truncation errors committed in each step.

Method	Global Truncation Error	Computational Cost
Euler	$\Delta t$	1
Euler-Cromer	$(\Delta t^2)^{\dagger}$	1
Second-order Runge-Kutta	$(\Delta t^2)$	2
Fourth-order Runge-Kutta	$(\Delta t^4)$	4
Verlet$^{\ddagger}$	$(\Delta t^3)$	1
Leapfrog	$(\Delta t^2)$	2

$^{\dagger}$: for general problems, the Euler-Cromer method can only promise $O(\Delta t)$, however for oscillation problems, this algorithm actually conservers energy over each complete oscillation at least $O([\Delta t]^2)$.
$^{\ddagger}$: This method requires knowledge of $\theta_{i-1}$ and is thus not self starting.

2.2.2 Derivation of RK4 for $2^{nd}$ order ODE[5]

$\ddot x + a\dot x + bx = f\left( t \right)\quad \quad x\left( {{t_0}} \right) = {x_0},v\left( {{t_0}} \right) = {v_0}$

substitue $v$ for $\dot{x}$ to obtain $\dot v + av + bx = f\left( t \right)$

by definition, $\frac{{dx}}{{dt}} = v$　　　　　　　　　　　$\therefore dx = vdt$

　　　　　　　$\frac{{dv}}{{dt}} = \dot v = {f_2}\left( {t,x,v} \right)$　　　　$\therefore dv = \dot vdt = {f_2}\left( {t,x,v} \right)$

where for the above eq'n ${f_2}\left( {t,x,v} \right)=f\left( t \right) - bx - av$

Next, apply Runge-Kutta formulas to each of the two first order equations as follows:

$\left\{ \begin{array}{l} {x_{i + 1}} = {x_i} + dx\\ {v_{i + 1}} = {v_i} + dv \end{array} \right.$

where $\left\{ \begin{array}{l} d{x_1} = {v_i}\Delta t\\ d{x_2} = \left( {{v_i} + \frac{{d{v_1}}}{2}} \right)\Delta t\\ d{x_3} = \left( {{v_i} + \frac{{d{v_2}}}{2}} \right)\Delta t\\ d{x_4} = \left( {{v_i} + d{v_3}} \right)\Delta t\\ dx = \frac{1}{6}\left( {d{x_1} + 2d{x_2} + 2d{x_3} + d{x_4}} \right) \end{array} \right.\quad \quad \left\{ \begin{array}{l} d{v_1} = {f_2}\left( {{t_i},{x_i},{v_i}} \right)\Delta t\\ d{v_2} = {f_2}\left( {{t_i} + \frac{{\Delta t}}{2},{x_i} + \frac{{d{x_1}}}{2},{v_i} + \frac{{d{v_1}}}{2}} \right)\Delta t\\ d{v_3} = {f_2}\left( {{t_i} + \frac{{\Delta t}}{2},{x_i} + \frac{{d{x_2}}}{2},{v_i} + \frac{{d{v_2}}}{2}} \right)\Delta t\\ d{v_4} = {f_2}\left( {{t_i} + \Delta t,{x_i} + d{x_3},{v_i} + d{v_3}} \right)\Delta t\\ dv = \frac{1}{6}\left( {d{v_1} + 2d{v_2} + 2d{v_3} + d{v_4}} \right) \end{array} \right.$

2.2.3 Solve $\frac{d^2\theta}{dt^2} = - \theta^3$

Here we use Euler-Cromer & Runge-Kutta 4 methods to sole this second order ordinary differential equation.

2.3 Conslusion

Up to now, we have analysed the harmonic and anharmonic oscillator analytically and numericaly, and we can draw the conslusion that the period of harmonic oscillator, that is $\alpha=1$, is independent of amplitude, as a contrast, the period of anharmonic oscillator , that is $\alpha = 3$, is dependent fo amplitude, and tin detail, the period is inversely propotional to the amplitude.

Reference: