Chapter 3 Problem 3.4

计算物理

---

Name: 陈锋
Student Number: 2013301020145
April 16, 2016

Abstract

This page is to solve the general form harmonic equation with $alpha = 1$ and analyze the relationship between the initial value of x and the period.

Introduction

The general form for the equation of motion for harmonic motion is:

$$\begin{equation} \frac{d^2x}{dt^2}=-kx^{\alpha} \end{equation}$$
When $\alpha = 1$, it becomes the case of simple harmonic oscillation, which has a constant period related to k. When $\alpha = 3$, the situation becomes different, especially for the period which is also related to the initial value of x (we call it amplitude sometimes). This page will focus on how the period depends on the amplitude and the constant k in the differential equation.

Numerical Approach

For the equation:

$$\begin{equation} \frac{d^2x}{dt^2}=-kx^{\alpha} \end{equation}$$
Then, we can rewrite this as two first-order differential equations:
$$\begin{equation} \frac{d\omega}{dt} = - kx^{\alpha}\\ \frac{dx}{dt} =\ \omega \end{equation}$$ where $\omega$ is similar to the angular velocity of the pendulum.
Considering the Euler_Crome method, I convert equation (2) into difference equations via using a time step $\Delta t$ so that rime is discretized with $t = i\Delta t$. Letting $x_i$ and $\omega_i$ be the numerically approximated angular displacement and velocity of the pendulum at time step i, the equations in (2) become:
$$\begin{equation} \omega_{i + 1} = \omega_i - k{x_i}^{\alpha}\Delta t\\ x_{i+1} = x_i +\omega_{i+1} \Delta t \end{equation}$$
Also, I will use 2nd-order Runge-Kutta method and 4th-order Runge-Kutta method to solve this equation and make a comparison. Their introduction have showed in ode.Instruction and the algorithms have been programmed in ode.py.

Data Analysis

The case of $\alpha = 1$:

Seting the initial value $x_0 = 1$, time step = 0.04s

Figure 1: In this case, the Euler-Cromer method is as precise as Second-order Runge-Kutta method and Forth-order Runge-Kutta method.

The case of $\alpha = 3$:

Setting time step = 0.04s, and the initial value change from 0.2 to 1.0:

Figure 2: As the figure showed above, the period becomes longer as the amplitude is decreased and I will give a detail discussion in the next section.

Find the period

This section is about how to find the period in the numerical results. Firstly, we have set the inital phase as zero, so the initial value became the ampilitude. Secondly, we should make a pragramme which can read all the datum and find out the first point except the initial point where the value is near the ampilitude. The corresponding time is the period we need to find. Generally, we might find two or even more points satisfying the former requirements, then we should make an average over these points in order to get a precise period. As a matter of fact, when the precision is high enough, which means the time step and the difference between the displacement and the ampilitude is sufficiently small, there are a little difference among these points. Thus, I only select the first two points and make an average over them.
The code are as follows:

    def find_period(self):
        t_period = []
        for i in range(len(self.x)):
            if abs(self.x[i] - self.x_0) < 0.000001 and self.t[i] > 20*self.dt:
                t_period.append(self.t[i])
            else:
                pass
        s = 0
        for j in range (2):
            s = s + t_period[j]
        print "the period is:%f" %(s/2)
        return s/2

Result:
Setting the time step = 0.001s

It seems that there is a relationship between the amplitude and the period and I try to figure out such relation via fitting in the next subsection.

Fitting result:

General model:
f(x) = a/x
Coefficients (with 95% confidence bounds):
a = 7.415 (7.413, 7.416)

Intuitive explantion

For the equation of angular velocity (x is smaller than 1):

$$\begin{equation} \frac{d\omega}{dt} = - k x^3 < -k x \end{equation}$$
We should notice that the initial value of angular veloctiy is zero. Also, according to conservation of energy, the displacement must smaller or equal the initial value (for the initial phase is equal to zero.). Therefore, the change of angular velocity in this case is much smaller than that of simple harmonic. As a direct result, the angular velocity will also smaller than that case.
Since:
$$\begin{equation} T = 4 \int_0^{x_0}\frac{dx}{\omega} \end{equation}$$
A smaller angular velocity means a longer period, so the period becomes longer as the amplitude is decreased.

Analytical method for the period

For the equation:

$$\begin{equation} \frac{d^2x}{dt^2}=-kx^{\alpha} \end{equation}$$
Setting $\omega = \frac{dx}{dt}$, then we have:
$$\begin{equation} \omega\frac{d\omega}{dx}=-kx^{\alpha} \end{equation}$$
which means:
$$\begin{equation} {\omega}d{\omega}=-kx^{\alpha}dx \end{equation}$$
Integrating both sides of it:
$$\begin{equation} \int_0^\omega \omega d\omega = -k \int_{x_0}^x x^{\alpha}dx \end{equation}$$
Then we have:
$$\begin{equation} \omega^2=\frac{2k}{\alpha+1}(x_0^{\alpha+1}-x^{\alpha+1}) \end{equation}$$
Therefore, we have:
$$\begin{equation} T=-4\sqrt{\frac{\alpha+1}{2k}}\int_{x_0}^0\frac{dx}{\sqrt{x_0^{\alpha+1}-x^{\alpha+1}}} \end{equation}$$ where T is the period.
when $\alpha = 1$
$$\begin{equation} T=-4\sqrt{\frac{1}{k}}\int_{x_0}^0\frac{dx}{\sqrt{x_0^{2}-x^{2}}}=-4\sqrt{\frac{1}{k}}arcsin\frac{x}{x_0}|^0_{x_0}=2\pi\sqrt{\frac{1}{k}} \end{equation}$$
when $\alpha = 3$
$$\begin{equation} T=-4\sqrt{\frac{2}{k}}\int_{x_0}^0\frac{dx}{\sqrt{x_0^{4}-x^{4}}}=4\sqrt{\frac{2}{k}}\frac{\sqrt{\pi}\Gamma(\frac{5}{4})}{x_0\Gamma(\frac{3}{4})} \approx 7.4163\frac{1}{x_0\sqrt{k}} \end{equation}$$
when setting $k = 1$, the coefficient is nearly the same as the numerical fitting result.

Reference

1.Giodano, N.J., Nakanishi, H. Computational Physics. Tsinghua University Press, December 2007.
2.Wikipedia contributors. "Pendulum (mathematics)." Wikipedia, The Free Encyclopedia. Wikipedia, The Free Encyclopedia, 7 Apr. 2016. Web. 16 Apr. 2016.

Programme Code

• Chapter 3 Problem 3.4