An implementation of John Conway's Game of Life in Python.
Damped Harmonic Oscillator - Derivation and solution of the differential equations
Damped Harmonic Oscillator
Problem Statement
The damped harmonic oscillator is a classic problem in mechanics. It describes the movement of a mechanical oscillator (eg spring pendulum) under the influence of a restoring force and friction. This article deals with the derivation of the oscillation equation for the damped oscillator.
Although a basic understanding of differential calculus is assumed, the aim of this article is to provide the derivation with as many details as possible. Unfortunately many other sources available on the Internet omit important secondary calculations or only present them in abbreviated form.
The equation of motion
An equation of motion is a mathematical equation that completely describes the spatial and temporal development of a dynamic system under the influence of external forces. As a rule, these are second-order differential equations. In this section, we first consider the forces that affect the movement of the pendulum.
The pendulum motion depends on the balance of three forces: Inertia, the restoring force of the spring and friction.
Inertia: (Newtons 2nd law) $$F = m \ddot {x}$$
Restoring force of the spring: $$F_f = - k x$$
Friction: $$F_r = - \mu \dot{x}$$
Newtons second law states that any change in the motion of a body is proportional to the force acting on it and that this change in movement always takes place in the same direction in which the force acts.
On a pendulum we have a restoring force \(F_f\) that is acting on the oscillator. It pulls the pendulum mass back towards its resting position. This force is directly proportional to the deflection of the pendulum. It's direction is opposite to the direction of pendulum deflection. It's strength is computed by multiplying the distance of deflection with a material dependent spring constant \(k\).
Lastly the pendulum shall experience friction. The friction force \(F_r\) is proportional to the velocity of the pendulum. Its direction is opposite to the the direction of movement of the pendulum. The parameter \(\mu\) is a material and shape dependent coefficient of friction. (see also Stokes law ). The inertial force is therefore opposed by two forces:
\begin{equation}
\nonumber m \ddot {x} = - \mu \dot{x} - k x \label{eq:eqn_of_motion1}
\end{equation}
which can be transformed into:
\begin{equation}
\ddot {x} + {\mu \over m} \dot{x} + {k \over m} x = 0 \label{eq:eqn_of_motion2}
\end{equation}
This equation is a linear, homogeneous, second-order differential equation with constant coefficients. Determining the type of differential equation is important because the approach necessary for solving it depends on it. The exponential ansatz is usually chosen to solve differential equations of this type.
The Exponential Ansatz
Equation \eqref{eq:eqn_of_motion2} is a homogeneous linear differential equation of second order with constant coefficients. The assumption is that the solution to such a differential equation is an exponential function. This method was first described by the german mathematician Leonard Euler. It is therefore called "exponential ansatz".
The exponential ansatz states that a special solution to the differential equation takes the form:
\begin{equation}
x(t) = Ce^{\lambda t} \label{eq:exp_sol1}
\end{equation}
By forming the first and second derivations of \eqref{eq:exp_sol1} we obtain:
\begin{equation}
\dot{x}(t) = \lambda Ce^{\lambda t}
\\
\ddot{x}(t) = \lambda^2 Ce^{\lambda t} \label{eq:exp_sol2}
\end{equation}
By inserting the equations \eqref{eq:exp_sol1} and \eqref{eq:exp_sol2} into the differential equation \eqref{eq:eqn_of_motion2} we obtain:
\begin{equation}
\nonumber \lambda^2 Ce^{\lambda t} + {\mu \over m} \lambda Ce^{\lambda t} + {k \over m} Ce^{\lambda t} = 0
\end{equation}
Division by \(Ce^{\lambda t}\) will simplify this to:
\begin{equation}
\lambda^2 + {\mu \over m} \lambda + {k \over m} = 0
\end{equation}
This equation is also called the characteristic equation. It is a quadratic equation in normal form. Therefore the following two solutions exist for \(\lambda\):
\begin{equation}
\lambda_{1,2} = - {\mu \over {2m} } \pm \sqrt { \Big({\mu \over {2m} }\Big)^2 - {k \over m} } \label{eq:lambda1}
\end{equation}
For further simplification we introduce the new constante \(\delta\) and \(\omega_0\):
\begin{equation}
\nonumber \delta = {\mu \over {2m} } \quad,\quad \omega_0 = \sqrt {k \over m}
\end{equation}
Equation \eqref{eq:lambda1} can now be written as:
\begin{equation}
\lambda_{1,2} = - \delta \pm \sqrt { \delta^2 - \omega_0^2 } \label{eq:lambda2}
\end{equation}
The further solution depends on the values of \(\lambda_1\) and \(\lambda_2\). We need to consider several different possibilities. The term inside the square root is called the discriminant. Depending on the choice of the constants \(\delta\) and \(\omega_0\) the discriminant can be greater than 0, less than 0 or equal to 0. Therefore \(\lambda_1\) and \(\lambda_2\) can be:
- Two real valued solutions that differ from each other.
- Two conjugate complex solutions.
- Two identical real valued solutions.
Each of the three options requires a different approach to the solution. The general solution of the homogeneous differential equation has the form:
\begin{equation}
x(t) = C_1 x_1(t) + C_2 x_2(t) \label{eq:general_sol}
\end{equation}
The functions \(x_1(t)\) and \(x_2(t)\) are determined by the value of the discriminant in equation \eqref{eq:lambda2}. Lets now look at the different possibilities.
Overdamped case
If \(\delta > \omega_0\) we have a case of strong friction. Therefore the discriminant in equation \eqref{eq:lambda2} is positive and two different real valued solutions exist for \(\lambda\). The solution to the differential equation then looks like:
\begin{equation}
x_1(t) = C_1e^{\lambda_1 t}
\\
\nonumber x_2(t) = C_2e^{\lambda_2 t}
\end{equation}
By inserting this into Equation \eqref{eq:general_sol} we obtain for the general solution of the differential equation:
\begin{equation}
\nonumber x(t) = C_1e^{\lambda_1 t} + C_2e^{\lambda_2 t}
\end{equation}
Which can be further transformed be inserting the lambdas from equation \eqref{eq:lambda2}:
\begin{equation}
x(t) = C_1e^{(-\delta + \sqrt {\delta^2 - \omega_0^2 })\;t} + C_2e^{(-\delta -\sqrt { \delta^2 - \omega_0^2 })\;t} \label{eq:general_sol2}
\end{equation}
Lets now have a look at the exponent in equation \eqref{eq:general_sol2}. The exponents must be negative values because for overdamping case the following must be true:
\begin{equation}
\delta > \sqrt {\delta^2 - \omega_0^2 }
\end{equation}
Therefore equation \eqref{eq:general_sol2} is an addition of two exponentially decaying functions. We will see those again in the diagram below. To simplify the equation we will substitute the square root with a new constant:
\begin{equation}
\alpha = \sqrt {\delta^2 - \omega_0^2 }
\end{equation}
After inserting the new constant in equation \eqref{eq:general_sol2} we can rewrite it as:
\begin{equation}
x(t) = e^{-\delta t} \Big(C_1e^{\alpha t} + C_2e^{-\alpha t}\Big) \label{eq:general_sol3}
\end{equation}
This is the solution for the overdamped case. The integration constants \(C_1\) and \(C_2\) for a special problem can be determined from given initial conditions. For instance the start position and the initial velocity of the pendulum could be given.
\begin{equation}
x(0) = x_0 \quad,\quad \dot{x}(0) = v_0
\end{equation}
The parameters do not represent overdamping! \( \delta \le \omega_0 \) !
Critically damped case
We speak of critical damping when \(\delta = \omega_0\). This is the transition from overdamping to the oscillation. In this case equation \eqref{eq:lambda2} has a single solution for \(\lambda\):
\begin{equation}
\lambda_1 = \lambda_2 = \lambda = -\delta \label{eq:lambda_aperiod}
\end{equation}
Using the exponential ansatz in this case means using the following approach to obtain the two partial solutions to the differential equation:
\begin{equation}
x_1(t) = C_1 e^{\lambda t}
\\
x_2(t) = t \; C_2 e^{\lambda t} \label{eq:ansatz_aperiod}
\end{equation}
Inserting equation \eqref{eq:lambda_aperiod} and \eqref{eq:ansatz_aperiod} in the general solution \eqref{eq:general_sol} yields the solution for the critically damped case:
\begin{equation}
x(t) = e^{-\delta t} (C_1 + t \; C_2)
\end{equation}
The integration constants \(C_1\) and \(C_2\) must be obtained from the initial conditions of a specific problem:
\begin{equation}
x(0) = x_0 \quad,\quad \dot{x}(0) = v_0
\end{equation}
The details for the calculation of the constants can be found at Wolfram Alpha.
Underdamped case
Damped oscillation occurs for \(\delta < \omega_0 \). In this case, the discriminant in equation \eqref{eq:lambda2} is negative. Therefore \(\lambda_1\) and \(\lambda_2\) are complex numbers. The exponential ansatz \(x(t)=C e^{\lambda t}\) is again used to solve the differential equation.
\begin{equation}
x_1(t) = C_1 e^{\lambda_1 t}
\\
\nonumber x_2(t) = C_2 e^{\lambda_2 t}
\end{equation}
By inserting the exponential approach into equation \eqref{eq:general_sol} and replacing \(\lambda\) with equation \eqref{eq:lambda2} we obtain:
\begin{equation}
x(t) = e^{-\delta t} \Big( C_1 e^{\sqrt{\delta^2 - \omega_0^2}\; t} + C_2 e^{-\sqrt{\delta^2 - \omega_0^2}\; t}\Big) \label{eq:solution_schwingfall}
\end{equation}
We have already seen this equation in a slightly modified form as equation \eqref{eq:general_sol3}. It is the general solution to the differential equation of the harmonic oscillator. We are now dealing with a complex valued solution because in the underdamped case the term under the square root is negative. Therefore the constants \(C_1\) and \(C_2\) are complex values. For dealing with complex numbers we are using Euler's formula. This equation establishes a connection between complex exponential functions and trigonometric functions:
\begin{equation}
\nonumber e^{i \phi} = cos \phi + i sin \phi
\end{equation}
It is helpfull to rearrange equation \eqref{eq:solution_schwingfall} in a way that its imaginary part is isolated:
\begin{equation}
\sqrt{\delta^2 - \omega_0^2} = \sqrt{-1 * (\omega_0^2 - \delta) } = i \sqrt{\omega_0^2 - \delta^2} \label{eq:iomega}
\end{equation}
We obtain a term consisting of the imaginary unit \(i\) multiplied with a real valued square root. To simplify further calculations we substitute the square root with a new constant named \(\omega\):
\begin{equation}
\omega = \sqrt{\omega_0^2 - \delta^2} \label{eq:omega}
\end{equation}
As we will see shortly, this constant represents the natural frequency of the harmonic oscillator. With the relationships from equations \eqref{eq:iomega} and \eqref{eq:omega} equation \eqref{eq:solution_schwingfall} can now be transformed into:
\begin{equation}
x(t) = e^{-\delta t} \Big( \underbrace{C_1 e^{i \omega t} + C_2 e^{- i \omega t}}_{\hat{x}(t)}\Big) \label{eq:sol_complex}
\end{equation}
From a physical point of view only purely real valued solution are of interest. To find these we have to separate the real and imaginary parts from one another. As already mentioned \(C_1\) and \(C_2\) are complex-valued constants. We are now transforming them into their polar form.
\begin{equation}
C_1 = \hat{C_1} e^{i \phi_1}
\\
C_2 = \hat{C_2} e^{i \phi_2} \label{eq:const_pol}
\end{equation}
The first part of equation \eqref{eq:sol_complex} is an exponentially decaying function. Because the constant \(\delta\) is not a complex value this part of the equation cannot yield complex results. We therefore focus on the second part which we label \(\hat{x}(t)\). We replace the constants with the ones written in polar form \eqref{eq:const_pol} and apply Euler's formula.
\begin{align}
\hat{x}(t) & = \hat{C_1} e^{i \phi_1} e^{i \omega t} + \hat{C_2} e^{i \phi_2} e^{- i \omega t}
\nonumber\\
& = \hat{C_1} e^{i (\phi_1 + \omega t)} + \hat{C_2} e^{i (\phi_2 - \omega t)}
\nonumber\\
& = \hat{C_1} cos(\phi_1 + \omega t) + i \hat{C_1} sin(\phi_1 + \omega t) + \hat{C_2} cos(\phi_2 - \omega t) + i \hat{C_2} sin(\phi_2 - \omega t)
\nonumber\\
& = \hat{C_1} cos(\phi_1 + \omega t) + \hat{C_2} cos(\phi_2 - \omega t) + i \Big(\hat{C_1} sin(\phi_1 + \omega t) + \hat{C_2} sin(\phi_2 - \omega t) \Big) \label {eq:sep_imag}
\end{align}
In equation \eqref{eq:sep_imag} real and imaginary parts of the equation are separated. We are only interested in solutions for which the imaginary part disappears. This happens when the following is true:
\begin{align}
\hat{C_1} sin(\phi_1 + \omega t) + \hat{C_2} sin(\phi_2 - \omega t) & = 0
\nonumber\\
\hat{C_1} sin(\phi_1 + \omega t) & = - \hat{C_2} sin(\phi_2 - \omega t)
\nonumber\\
\hat{C_1} sin(\phi_1 + \omega t) & = \hat{C_2} sin( - \phi_2 + \omega t ) \label{eq:purge_imag}
\end{align}
Thus the imaginary part of equation \eqref{eq:sol_complex} disappears if:
\begin{equation}
\hat{C_1} = \hat{C_2} \; , \; \phi_1 = - \phi_2 \label{eq:complex_const2}
\end{equation}
Or in other words: the constants \(C_1\) and \(C_2\) must be complex conjugate:
\begin{equation}
\nonumber C_1 = C_2^*
\end{equation}
By inserting equation \eqref{eq:complex_const2} in the remaining real part of equation \eqref{eq:sep_imag} we obtain for \(\hat{x}(t)\):
\begin{align}
\hat{x}(t) & = \hat{C_1} cos(\phi_1 + \omega t) + \hat{C_1} cos(-\phi_1 - \omega t)
\nonumber \\
& = \hat{C_1} cos(\phi_1 + \omega t) + \hat{C_1} cos(\phi_1 + \omega t)
\nonumber \\
& = 2 \hat{C_1} cos(\phi_1 + \omega t) \label{eq:xhat_final}
\end{align}
By using equation \eqref{eq:xhat_final} in equation \eqref{eq:sol_complex} we obtain the solution of the differential equation for the underdamped case case:
\begin{equation}
x(t) = e^{-\delta t} \Big( 2 A cos(\phi + \omega t) \Big) \label{eq:sol_complex_final}
\end{equation}
The constants \(\phi_1\) and \(\hat{C_1}\) have been renamed to \(A\) und \(\phi\). \(A\) is an amplitude and \(\phi\) the phase angle. The constants can be computed from a given set of initial conditions, such as:
\begin{equation}
x(0) = x_0 \quad,\quad \dot{x}(0) = v_0
\end{equation}
The entire calculations for computing the constants for a given set of \(v_0\) and \(x_0\) can be found here. The following graphic shows the amplitude curve calculated over time.
The parameters do not represent an underdamped case \( \omega_0 \le \delta \)!
Graphical representation of the general analytical solution of the underdamped case. The diagram shows the pendulum deflection over time for the initial conditions set in the side panel.
References
- "Der Exponentialansatz zur Lösung der linearen, homogenen Dgl. 2. Ordnung." Lukas Geiling, Martin Wagener (Dr. W. Seifert), Lehrdokument, Institute of Physics, University Halle-Wittenberg, D-06099 Halle, Germany
- "Gedämpfte harmonische Schwingungen." Barbara Herzog, Kim Holm, Lehrdokument, Institute of Physics, University Halle-Wittenberg, D-06099 Halle, Germany
- "Freie gedämpfte Schwingung", Online-Angebot der Georg-August-Universität Göttingen, Georg-August-Universität Göttingen
![[1]](http://commons.wikimedia.org/w/index.php?title=File:Damped_oscillation_graph2.svg&oldid=17548606){kind=link}
![[2]](http://commons.wikimedia.org/w/index.php?title=File:Simple_harmonic_oscillator.gif&oldid=20055037){kind=link}