Preface

This section presenmts explicit solution of simple pendulum equation through Jacobi function.

Return to computing page for the first course APMA0330
Return to computing page for the second course APMA0340
Return to Mathematica tutorial for the second course APMA0340
Return to the main page for the course APMA0330
Return to the main page for the course APMA0340
Return to Part IV of the course APMA0330

Glossary

Solution of pendulum equation

Simple pendulum equation \( \ddot{\theta} + \omega_0^2 \sin \theta =0 , \) although straightforward in appearance, is in fact rather difficult to solve because of the non-linearity of the term \( \sin \theta . \) It turns out that the general initial value problem

\begin{equation} \ddot{\theta} + \omega_0^2 \sin\theta =0 , \qquad \theta (0) =a, \quad \dot{\theta} (0) =b , \label{Eq.sol.1} \end{equation}
where 𝑎 is the initial displacement and b is the initial velocity of the bob (its mass does not matter)), could not be solved using elementary functions but is expressed through elliptic functions. Here &omega0² = ℓ/g, with ℓ being pendulum length and g being acceleration due to gravity.Accordingly, θ denotes the angle between the rod and the downward vertical direction with the counterclockwise direction taken as positive.

As usual in this tutorial, we demonstrate different approaches to obtain the solution. In order to derive the exact solution, this equation is multiplied by the integrating factor \( \dot{\theta} = {\text d}\theta / {\text d} t, \) so that it becomes

\[ \dot{\theta} \,\ddot{\theta} + \dot{\theta}\, \omega_0^2 \sin\theta =0 , \]
which can be rewritten as
\[ \frac{{\text d}}{{\text d}t} \left[ \frac{1}{2} \left( \frac{{\text d} \theta}{{\text d}t} \right)^2 - \omega_0^2\, \cos\theta \right] =0 . \]
This equation, which corresponds to the conservation of the mechanical energy, is immediately integrable, taking into account the initial conditions:
\[ \left( \frac{{\text d} \theta}{{\text d}t} \right)^2 = b^2 + 2\,\omega_0^2\left( \cos\theta - \cos a \right) = b^2 - 2\,\omega_0^2 \cos a + 2\,\omega_0^2 \cos\theta , \]
which leads to two equations for each branch of square root:
\[ \frac{{\text d} \theta}{{\text d}t} = \pm \sqrt{b^2 + 2\,\omega_0^2\left( \cos\theta - \cos a \right) } . \]
From this equation we immediately can determine the maximum amplitude of oscillations by equating the first derivative to zero:
\[ \cos\theta = \cos a - \frac{b^2}{2\,\omega_0^2} . \qquad\blacktriangleleft \]

We are going to derive the above first order differential equation for angle θ using another approach. Let \( \alpha = \dot{\theta} = {\text d}\theta / {\text d}t \) be the derivative of θ with respect to time variable t. Then from chain rule, we get

\[ \ddot{\theta} = \frac{\text d}{{\text d}t} \left( \frac{{\text d}\theta}{{\text d}t} \right) = \frac{\text d}{{\text d}t} \left( \dot{\theta} \right) = \frac{\text d}{{\text d}t} \left( \alpha \right) = \frac{{\text d}\alpha}{{\text d}\theta} \, \frac{{\text d}\theta}{{\text d}t} = \alpha\,\frac{{\text d}\alpha}{{\text d}\theta} . \]
Then the pendulum equation becomes
\[ \alpha\,\frac{{\text d}\alpha}{{\text d}\theta} + \omega_0^2 \sin\theta \qquad \Longrightarrow \qquad \alpha \, {\text d}\alpha = - \omega_0^2 \sin\theta \,{\text d} \theta . \]
Since the latter is a separable differential equation of first order, we integrate both sides to obtain
\[ \frac{\alpha^2}{2} - \frac{b^2}{2} = - \omega_0^2 \int_a^{\theta} \sin\theta \,{\text d} \theta = \omega_0^2 \left( \cos\theta - \cos a \right) . \]
Now we return to our original variable \( \dot{\theta} = \alpha \) and extract square root:
\[ \dot{\theta} = \pm \sqrt{b^2 + 2\omega_0^2 \cos\theta - 2 \omega_0^2 \cos a} , \]
which is exactly the same as we obtained previously using consevation of energy approach. Since the above first order differential equation is separable
\[ \pm\frac{{\text d}\theta}{\sqrt{A + B\,\cos\theta}} = {\text d}t , \qquad \mbox{where} \quad A = b^2 - 2\,\omega_0^2 \cos a , \quad B = 2\,\omega_0^2 , \]
we ask Mathematica to integrate
Integrate[(A + B*Cos[x])^(-1/2), x]
(2 Sqrt[(A + B Cos[x])/(A + B)] EllipticF[x/2, (2 B)/(A + B)])/Sqrt[A + B Cos[x]]
Here EllipticF denotes the elliptic integral of the first kind
\[ F\left( \phi , k \right) = \int_0^{\phi} \frac{{\text d}x}{\sqrt{1 - k^2 \sin^2 x}} . \]
The inverse function of F(φ,k) is given by the Jacobi amplitude
\[ \mbox{am} (u,k) = \phi = F^{-1} \left( u , k \right) . \]
This allows us to express the solution of the pendulum equation only implicitly:
\[ \frac{2}{\sqrt{b^2 - 2\omega_0^2 \cos a + 2\omega_0^2}} \, F \left( \frac{\theta}{2} , \frac{4\omega_0^2}{b^2 - 2\omega_0^2 \cos a + 2\omega_0^2} \right) = \frac{2}{\sqrt{b^2 - 2\omega_0^2 \cos a + 2\omega_0^2}} \, F \left( \frac{a}{2} , \frac{4\omega_0^2}{b^2 - 2\omega_0^2 \cos a + 2\omega_0^2} \right) = \pm t . \]

Even with the aid of special function (the elliptic integral of the first kind), the obtained implicit solution is not pleasant to deal with. Therefore, we consider two auxiliary initial value problems:

Problem A:

\[ \ddot{\theta} + \omega_0^2 \sin\theta =0, \qquad \theta (0) =a , \quad \dot{\theta} (0) =0 ; \]

Problem B:

\[ \ddot{\theta} + \omega_0^2 \sin\theta =0, \qquad \theta (0) =0 , \quad \dot{\theta} (0) = b , \]
It should be noted that the original initial value problem is not a sum of solutions to these two problems because the pendulum equation is nonlinear.

Upon considering problem A, we reduce it to the first order differential equation:

\[ \left( \frac{{\text d} \theta}{{\text d}t} \right)^2 = 4\,\omega_0^2\left( \sin^2 \left( \frac{\theta_0}{2} \right) - \sin^2 \left( \frac{\theta}{2} \right) \right] , \]
where we used the trigonometric identity \( \cos \theta = 1 - 2\,\sin^2 \frac{\theta}{2} . \) Now let
\[ y = \sin \left( \frac{\theta}{2} \right) \qquad\mbox{and} \qquad k= \sin \left( \frac{\theta_0}{2} \right) . \]
Then \( y(0) = \sqrt{k} . \) It is easy to obtain \( {\text d}\theta / {\text d}t \) as a function of \( {\text d}y / {\text d}t \) as follows:
\[ \frac{{\text d}y}{{\text d}t} = \frac{{\text d}y}{{\text d}\theta} \, \frac{{\text d}\theta}{{\text d}t} = \frac{1}{2}\, \frac{{\text d}\theta}{{\text d}t} \,\cos \left( \frac{\theta}{2} \right) . \]
Then square it, we have
\[ \left( \frac{{\text d}y}{{\text d}t} \right)^2 = \frac{1}{4} \, \cos^2 \left( \frac{\theta}{2} \right) \left( \frac{{\text d}\theta}{{\text d}t} \right)^2 = \frac{1}{4} \left[ 1 - \sin^2 \left( \frac{\theta}{2} \right) \right] \left( \frac{{\text d}\theta}{{\text d}t} \right)^2 = \frac{1}{4} \left( 1 - y^2 \right) \left( \frac{{\text d}\theta}{{\text d}t} \right)^2 . \]
Its square derivative will be
\[ \left( \frac{{\text d}\theta}{{\text d}t} \right)^2 = \frac{4}{1-y^2} \left( \frac{{\text d}y}{{\text d}t} \right)^2 . \]
Substituting this relation into the integrated pendulum equation, we obtain
\[ \frac{4}{1-y^2} \left( \frac{{\text d}y}{{\text d}t} \right)^2 = 4\,\omega_0^2 \left( k - y^2 \right) , \]
which can be rewritten as
\[ \left( \frac{{\text d}y}{{\text d}t} \right)^2 = \omega_0^2 k \left( 1 - y^2 \right) \left( 1 - \frac{y^2}{k} \right) . \]
To simplify the above expression, we define new variables
\[ \tau = \omega_0 t \qquad\mbox{and}\qquad z = \frac{y}{\sqrt{k}} . \]
This leads to
\[ \left( \frac{{\text d}z}{{\text d}\tau} \right)^2 = \left( 1 - z^2 \right) \left( 1 - k\,z^2 \right) , \qquad z(0) = 1, \quad \left. \frac{{\text d}z}{{\text d}\tau} \right\vert_{\tau =0} =0. \]
Separation of variables yields
\[ {\text d}\tau = \pm \frac{{\text d}z}{\sqrt{\left( 1 - z^2 \right) \left( 1 - k\,z^2 \right)}} . \]
The time τ to go from point (1,0) to point \( \left( z, {\text d}z / {\text d} \tau \right) \) in the lower half-plane of the graph of \( {\text d}z / {\text d} \tau \) as a function of z is
\[ \tau = -\int_1^z \frac{{\text d}\zeta}{\sqrt{\left( 1 - \zeta^2 \right) \left( 1 - k\,\zeta^2 \right)}} . \]
We rewrite the integral in the right hand side as
\[ \tau = \int_0^1 \frac{{\text d}\zeta}{\sqrt{\left( 1 - \zeta^2 \right) \left( 1 - k\,\zeta^2 \right)}} - \int_0^z \frac{{\text d}\zeta}{\sqrt{\left( 1 - \zeta^2 \right) \left( 1 - k\,\zeta^2 \right)}} , \]
which allows us to obtain τ as a function of z and k as
\[ \tau (z) = K(k) - F(\arcsin z, k ) , \]
where K(m) and F(φ , m) are the complete and the incomplete elliptic integral of the first kind defined as follows
\begin{eqnarray*} K(m) &=& \int_0^1 \frac{{\text d}\zeta}{\sqrt{\left( 1 - \zeta^2 \right) \left( 1 - m\,\zeta^2 \right)}} \\ F(\phi , m) &=& \int_0^{\phi} \frac{{\text d}\zeta}{\sqrt{\left( 1 - \zeta^2 \right) \left( 1 - m\,\zeta^2 \right)}} , \qquad z= \sin \phi . \end{eqnarray*}
The period of oscillation T is four times the time taken by the pendulum to swing from \( \theta =0 \ (z=0) \) to \( \theta = \theta_0 \ (z=1) . \) Therefore,
\[ T = 4\,t(0) = \frac{4\,\tau (0)}{\omega_0} = \frac{4}{\omega_0} \,K(k) = \frac{2}{\pi}\, T_0 \,K(k) = \sqrt{\frac{\ell}{g}} \,\int_0^{\pi /2} \frac{{\text d} \beta}{\sqrt{1- k^2 \sin^2 \beta}} , \]
where
\[ T_0 = \frac{2\pi}{\omega_0} = 2\pi \,\sqrt{\frac{\ell}{g}} \]
is the period of the pendulum for small oscillations. The solution formula can be rewritten as
\[ F \left( \arcsin z , k \right) = K(k) - \tau , \]
which can be written in terms of the Jacobi elliptic function sn(u;m):
\[ z = \mbox{sn} \left( K(k) - \tau ; k \right) , \]
where Mathematica has special commands for their evaluations: EllipticK[m] for complete elliptic integral of the first kind, and EllipticF[φ,m] for incomplete elliptic integral of the first kind. Mathematica also uses a special command JacobiSN[u,m] for the Jacobi elliptic function sn(u;m). These functions are discussed in the second part of the course.

Now we return to the original variables
\[ \sin \left( \frac{\theta}{2} \right) = \sin \frac{\theta_0}{2} \, \mbox{sn} \left[ K\left( \sin^2 \frac{\theta_0}{2} \right) - \omega_0 t ; \sin^2 \frac{\theta_0}{2} \right] , \]
which allows us to express θ as a function of t:
\begin{equation} \theta (t) = 2\,\arcsin \left\{ \sin\frac{\theta_0}{2} \, \mbox{sn} \left[ K\left( \sin^2 \frac{\theta_0}{2} \right) - \omega_0 t ; \sin^2 \frac{\theta_0}{2} \right] \right\} . \label{Eq.sol.2} \end{equation}
Using the relations \( \dot{\bf e}_{r} = \dot{\theta} {\bf e}_{\theta} \) and \( \dot{\bf e}_{\theta} = -\dot{\theta} {\bf e}_{r} , \) we obtain the velocity, acceleration, and jerk:
\begin{eqnarray*} {\bf v} (t) &=& \ell \dot{\theta} {\bf e}_{\theta} = \sqrt{2\ell g} \left( \cos \theta - \cos \theta_0 \right)^{1/2} {\bf e}_{\theta} , \\ {\bf a} (t) &=& \dot{\bf v} (t) = -\ell \dot{\theta}^2 {\bf e}_{r} + \ell \ddot{\theta} {\bf e}_{\theta} = - 2g \left( \cos \theta - \cos \theta_0 \right) {\bf e}_{r} - g \,\sin \theta \, {\bf e}_{\theta} , \\ {\bf j} (t) &=& \dot{\bf a} (t) = 3g\,\sin\theta \, \dot{\theta} \, {\bf e}_r - \dot{\theta} \left[ 2\ell \dot{\theta}^2 + g\,\sin \theta \right] {\bf e}_{\theta} . \end{eqnarray*}
Manipulate[
Module[{run, sol, z, app, a},
If[run == 10, run = 0];
sol = DSolve[{y''[x] + (9.8/l)*y[x] == 0, y[0] == int, y'[0] == 0}, y, x];
z[x_] := y[x] /. sol;
app = NDSolve[{\[Theta]1''[t] + (9.8/l)*Sin[\[Theta]1[t]] == 0, \[Theta]1[0] == int, \[Theta]1'[0] == 0}, \[Theta]1, {t, 0, 10}];
a[t_] := \[Theta]1[t] /. First[app];
Grid[{{
Graphics[{
{Thick, Line[{{0, 0}, {l*Sin[a[r]], -l*Cos[a[r]]}}]}, {Darker[Green], Disk[{0, 0}, .03]}, {Purple,
Disk[{l*Sin[a[r]], -l*Cos[a[r]]}, 0.1 Power[(30 m)/(4 \[Pi]), (3)^-1]]}},
PlotRange -> {{-3.4, 3.4}, {.5, -3.4}},
ImageSize -> {250, 250}],
Plot[{a[t], z[t]}, {t, 0, 10},
PlotRange -> {{0, 10}, {-1.6, 1.6}},
PlotStyle -> {{Thick, Blue}, {Red}}, ImageSize -> {250, 250},
AspectRatio -> 1.2, Frame -> True,
PlotLegends`PlotLegend -> {"numerical", "sin(\[Theta]) = \[Theta]"},
FrameLabel -> {"time", "displacement"},
Epilog -> {Blue, PointSize[.06], Point[{r, a[r]}]}]}}]],
{{int, \[Pi]/3., "initial angle \[Theta](0)"}, \[Pi]/16., \[Pi]/2., Appearance -> "Labeled"},
Delimiter,
{{m, 4, "mass m"}, 1, 10, .01, Appearance -> "Labeled"},
{{l, 2.25, "length L"}, 1, 3, .01,
Appearance -> "Labeled"}, {{r, 0, "release system"}, 0, 10, .01,
ControlType -> Trigger, AnimationRate -> .5},
TrackedSymbols :> {int, m, l, r},
Initialization :> Get["PlotLegends`"]]

Now we consider problem B. Its solution follows from the general formula:

\[ \frac{2}{b} \, F \left( \frac{\theta}{2} , \frac{4\omega_0^2}{b^2} \right) = \pm t . \]
Expressing θ from it, we obtain
\begin{equation} \theta (t) = 2\,\arcsin \left\{ k\,\mbox{sn} \left( \omega_0 t, k^2 \right) \right\} , \qquad k = \frac{b}{2\,\omega_0} . \label{Eq.sol.3} \end{equation}

 

  1. A. Belendez, C. Pascual, D.I. Mendez, T. Belendez and C. Neipp, Exact solution for the nonlinear pendulum, Revista Brasileira de Ensino de Fisica, v. 29, n. 4, p. 645--648, 2007.
  2. Alvaro H. Salas and Jairo E. Castillo H., Exact Solution to Duffing Equation and the Pendulum Equation, Applied Mathematical Sciences, Vol. 8, 2014, no. 176, 8781--8789.

 

Return to Mathematica page
Return to the main page (APMA0330)
Return to the Part 1 (Plotting)
Return to the Part 2 (First Order ODEs)
Return to the Part 3 (Numerical Methods)
Return to the Part 4 (Second and Higher Order ODEs)
Return to the Part 5 (Series and Recurrences)
Return to the Part 6 (Laplace Transform)
Return to the Part 7 (Boundary Value Problems)