Note: This page has been translated by MathWorks. Click here to see

To view all translated materials including this page, select Country from the country navigator on the bottom of this page.

To view all translated materials including this page, select Country from the country navigator on the bottom of this page.

This example simulates and explores the behavior of a simple pendulum by deriving its equation of motion, and solving the equation analytically for small angles and numerically for any angle.

Derive the Equation of Motion

Linearize the Equation of Motion

Solve Equation of Motion Analytically

Physical Significance of $${\omega}_{0}$$

Plot Pendulum Motion

Understanding Non-Linear Pendulum Motion Using Constant Energy Paths

Solve Non-Linear Equation of Motion

Solve Equation of Motion for Closed Energy Contours

Solve Equation of Motion for Open Energy Contours

Conclusion

The pendulum is a simple mechanical system that follows a differential equation. The pendulum is at rest in a vertical position. We displace the pendulum by an angle $$\theta $$ and release it. The force of gravity pulls it back towards its resting position, its momentum causes it to overshoot and come to an angle $$-\theta $$ (if there are no frictional forces), and so forth. The restoring force is $$-mg\mathrm{sin}\theta $$, the gravitational force along the pendulum's motion (with a minus sign to remind us that it pulls back to the vertical position). Thus, according to Newton's second law, the mass times the acceleration must equal $$-mg\mathrm{sin}\theta $$.

syms m a g theta(t) eqn = m*a == -m*g*sin(theta)

`eqn(t) = $$a\hspace{0.17em}m=-g\hspace{0.17em}m\hspace{0.17em}\mathrm{sin}\left(\theta \left(t\right)\right)$$`

The acceleration of the mass at the end of the pendulum is

$$\mathit{a}=\mathit{r}\text{\hspace{0.17em}}\frac{{\mathit{d}}^{2}\theta}{{\mathrm{dt}}^{2}}.$$

Substitute for `a`

by using `subs`

.

```
syms r
eqn = subs(eqn,a,r*diff(theta,2))
```

eqn(t) =$$m\hspace{0.17em}r\hspace{0.17em}\frac{{\partial}^{2}}{\partial {t}^{2}}\mathrm{}\theta \left(t\right)=-g\hspace{0.17em}m\hspace{0.17em}\mathrm{sin}\left(\theta \left(t\right)\right)$$

Isolate the angular acceleration in `eqn`

by using `isolate`

.

eqn = isolate(eqn,diff(theta,2))

eqn =$$\frac{{\partial}^{2}}{\partial {t}^{2}}\mathrm{}\theta \left(t\right)=-\frac{g\hspace{0.17em}\mathrm{sin}\left(\theta \left(t\right)\right)}{r}$$

Collect constants $$g$$ and $$r$$ into a single parameter, called the *natural frequency,*

$${\omega}_{0}=\sqrt{\frac{\mathit{g}}{\mathit{r}}}.$$

```
syms omega_0
eqn = subs(eqn,g/r,omega_0^2)
```

eqn =$$\frac{{\partial}^{2}}{\partial {t}^{2}}\mathrm{}\theta \left(t\right)=-{{\omega}_{0}}^{2}\hspace{0.17em}\mathrm{sin}\left(\theta \left(t\right)\right)$$

This equation is difficult to solve analytically because it is non-linear. Assuming the angles are small, we can linearize the equation by using the Taylor expansion of $$\mathrm{sin}\theta $$.

syms x approx = taylor(sin(x),x,'Order',2); approx = subs(approx,x,theta(t))

`approx = $$\theta \left(t\right)$$`

The equation of motion becomes

eqnLinear = subs(eqn,sin(theta(t)),approx)

eqnLinear =$$\frac{{\partial}^{2}}{\partial {t}^{2}}\mathrm{}\theta \left(t\right)=-{{\omega}_{0}}^{2}\hspace{0.17em}\theta \left(t\right)$$

Solve the equation `eqnLinear`

by using `dsolve`

. Specify initial conditions as the second argument.

syms theta_0 theta_t0 theta_t = diff(theta); cond = [theta(0) == theta_0, theta_t(0) == theta_t0]; thetaSol(t) = dsolve(eqnLinear,cond)

thetaSol(t) =$$\frac{0.5000\hspace{0.17em}{e}^{-1\hspace{0.17em}{\omega}_{0}\hspace{0.17em}t\hspace{0.17em}i}\hspace{0.17em}\left({\omega}_{0}\hspace{0.17em}{\theta}_{0}+1\hspace{0.17em}{\theta}_{\mathrm{t0}}\hspace{0.17em}i\right)}{{\omega}_{0}}-\frac{0.5000\hspace{0.17em}{e}^{1\hspace{0.17em}{\omega}_{0}\hspace{0.17em}t\hspace{0.17em}i}\hspace{0.17em}\left(-{\omega}_{0}\hspace{0.17em}{\theta}_{0}+1\hspace{0.17em}{\theta}_{\mathrm{t0}}\hspace{0.17em}i\right)}{{\omega}_{0}}$$

Simplify the solution by assuming $${\omega}_{0}$$ is real using `assume`

.

```
assume(omega_0,'real')
thetaSol(t) = dsolve(eqnLinear,cond)
```

thetaSol(t) =$${\theta}_{0}\hspace{0.17em}\mathrm{cos}\left({\omega}_{0}\hspace{0.17em}t\right)+\frac{{\theta}_{\mathrm{t0}}\hspace{0.17em}\mathrm{sin}\left({\omega}_{0}\hspace{0.17em}t\right)}{{\omega}_{0}}$$

$${\omega}_{0}t$$ is called the *phase*. The cosine function $\mathrm{cos}\text{\hspace{0.17em}}{\omega}_{0}\mathit{t}$ repeats every $$2\pi $$. The time needed to change $${\omega}_{0}t$$ by $$2\pi $$ is called the time period

$${\mathit{t}}_{0}=\frac{2\pi}{{\omega}_{0}}=2\pi \sqrt{\frac{\mathit{r}}{\mathit{g}}}.$$

From the equation, the time period is directly proportional to the pendulum's length. However, $${t}_{0}$$ does not depend on the mass because its moment of inertia and its weight are both directly proportional to its mass.

The time period does not depend on the initial conditions but the amplitude does. Instead, the period is governed by the equation of motion.

Plot the motion of the pendulum for the small angle approximation.

Define physical parameters.

gValue = 9.81; % m / s^2 rValue = 1; % m omega_0Value = sqrt(gValue/rValue); t_0 = 2*pi/omega_0Value;

Set initial conditions.

theta_0Value = 0.1*pi; % Solution only valid for small angles. theta_t0Value = 0; % Initially at rest.

Substitute physical parameters and initial conditions into the general solution.

vars = [omega_0 theta_0 theta_t0]; values = [omega_0Value theta_0Value theta_t0Value]; thetaSolPlot = subs(thetaSol,vars,values); figure; fplot(thetaSolPlot(t*t_0)/pi, [0 1]);

fplot(thetaSolPlot(t*t_0)); grid on; title('Harmonic Pendulum Motion'); xlabel('t/t_0'); ylabel('\theta/\pi');

Having found the solution for $$\theta (t)$$, we can visualize the motion of the pendulum below.

To understand the non-linear motion of a pendulum, visualize the pendulum path by using the equation for total energy. Since the non-linear motion must conserve total energy, paths that have constant energy describe the non-linear motion.

The total energy is

$$\mathit{E}=\frac{1}{2}\mathit{m}\text{\hspace{0.17em}}{\mathit{r}}^{2}{\left(\frac{\mathit{d}\theta}{\mathrm{dt}}\right)}^{2}+\mathit{m}\text{\hspace{0.17em}}\mathit{g}\text{\hspace{0.17em}}\mathit{r}\text{\hspace{0.17em}}\left(1-\mathrm{cos}\text{\hspace{0.17em}}\theta \right).$$

We can use the trigonometric identity

$$\left(1-\mathrm{cos}\text{\hspace{0.17em}}\theta \right)=2\text{\hspace{0.17em}}{\mathrm{sin}}^{2}\frac{\theta}{2}$$

Use the relation $${\omega}_{0}=\sqrt{g/r}$$ to write the scaled energy as

$$\frac{\mathit{E}}{\mathit{m}\text{\hspace{0.17em}}{\mathit{r}}^{2}}=\frac{1}{2}\left[{\left(\frac{\mathit{d}\theta}{\mathrm{dt}}\right)}^{2}+{\left(2{\omega}_{0}\mathrm{sin}\frac{\theta}{2}\right)}^{2}\right].$$

Since energy is conserved, the pendulum's motion can be described by constant energy paths in phase space, which is the abstract space with coordinates $$\theta $$ vs. $$d\theta /dt$$. We can visualize these paths using `fcontour`

.

syms theta theta_t omega_0 E(theta, theta_t, omega_0) = (1/2)*(theta_t^2+(2*omega_0*sin(theta/2))^2); Eplot(theta, theta_t) = subs(E,omega_0,omega_0Value); figure; fc = fcontour(Eplot(pi*theta, 2*omega_0Value*theta_t), 2*[-1 1 -1 1], ... 'LineWidth', 2, 'LevelList', 0:5:50, 'MeshDensity', 1+2^8); grid on; title('Constant Energy Contours in Phase Space ( \theta vs. \theta_t )'); xlabel('\theta/\pi'); ylabel('\theta_t/2\omega_0');

The curves are symmetric about the $$\theta $$ axis and $$d\theta /dt$$ axis, and are periodic along the $$\theta $$ axis. There are two regions of distinct behavior. The lower energies of the contour plot close upon themselves: the pendulum swings back and forth between two maximum angles and velocities.

The higher energies of the contour plot do not close upon themselves because the pendulum always moves in one angular direction.

The non-linear equations of motion are a second-order differential equation. Numerically solve these equations by using the `ode45`

solver. Because `ode45`

only accepts first-order systems, reduce the system to a first-order system. Then, generate function handles that are the input to `ode45`

.

Rewrite the second-order ODE as a system of first-order ODEs

syms theta(t) theta_t(t) omega_0 eqs = [diff(theta) == theta_t; diff(theta_t) == -omega_0^2*sin(theta)]

eqs(t) =$$\left(\begin{array}{c}\frac{\partial}{\partial t}\mathrm{}\theta \left(t\right)={\theta}_{t}\left(t\right)\\ \frac{\partial}{\partial t}\mathrm{}{\theta}_{t}\left(t\right)=-{{\omega}_{0}}^{2}\hspace{0.17em}\mathrm{sin}\left(\theta \left(t\right)\right)\end{array}\right)$$

eqs = subs(eqs,omega_0,omega_0Value); vars = [theta, theta_t];

Find the mass matrix `M`

of the system and the right sides of the equations `F`

.

[M,F] = massMatrixForm(eqs,vars)

M =$$\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)$$

F =$$\left(\begin{array}{c}{\theta}_{t}\left(t\right)\\ -9.8100\hspace{0.17em}\mathrm{sin}\left(\theta \left(t\right)\right)\end{array}\right)$$

`M`

and `F`

refer to the form

$$\mathit{M}\left(\mathit{t},\mathit{x}\left(\mathit{t}\right)\right)\frac{\mathrm{dx}}{\mathrm{dt}}=\mathit{F}\left(\mathit{t},\mathit{x}\left(\mathit{t}\right)\right).$$

To simplify further computations, rewrite the system in the form

$$\frac{\mathrm{dx}}{\mathrm{dt}}=\mathit{f}\left(\mathit{t},\mathit{x}\left(\mathit{t}\right)\right),$$

f = M\F

f =$$\left(\begin{array}{c}{\theta}_{t}\left(t\right)\\ -9.8100\hspace{0.17em}\mathrm{sin}\left(\theta \left(t\right)\right)\end{array}\right)$$

Convert `f`

to a MATLAB function handle by using `odeFunction`

. The resulting function handle is input to the MATLAB ODE solver `ode45`

.

f = odeFunction(f, vars)

`f = `*function_handle with value:*
@(t,in2)[in2(2,:);sin(in2(1,:)).*(-9.81e+2./1.0e+2)]

Solve the ODE for the closed energy contours by using `ode45`

.

From the energy contour plot, closed contours satisfy the condition $${\theta}_{0}=0$$, $${\theta}_{t0}/2{\omega}_{0}\le 1$$. Store the initial conditions of $$\theta $$ and $$d\theta /dt$$ in the variable `x0`

.

x0 = [0; 2*omega_0Value];

Choose a time interval that allows the pendulum to go through a full period. This can be found by trial and error.

tInit = 0; tFinal = 3.365*t_0;

Solve the ODE.

[t, x] = ode45(f, [tInit tFinal], x0);

$$\theta $$ is returned in the first column of `x`

, and $$d\theta /dt$$ is returned in the second column of `x`

.

Scale the time dimension by the period $${t}_{0}$$,

$${\mathit{t}}_{0}=\frac{2\pi}{{\omega}_{0}}=2\pi \text{\hspace{0.17em}}\sqrt{\frac{\mathit{r}}{\mathit{g}}},$$

t = t/t_0;

Make the values of $$\theta $$ repeat on the interval $$(-\pi ,\phantom{\rule{0.16666666666666666em}{0ex}}\pi )$$ using `wrapToPi`

, then scale $$\theta $$ by $$\pi $$.

x(:,1) = wrapToPi(x(:,1))/pi;

Scale $$d\theta /dt$$ by $$2{\omega}_{0}$$.

x(:,2) = x(:,2)/(2*omega_0Value);

Plot the closed path solution.

figure; yyaxis left; plot(t, x(:,1), '-o'); ylabel('\theta/\pi'); yyaxis right; plot(t, x(:,2), '-o'); ylabel('\theta_t/2\omega_0'); grid on; title('Closed Path in Phase Space'); xlabel('t/t_0');

This is how the pendulum motion would appear:

Solve the ODE for the open energy contours by using `ode45`

.

From the energy contour plot, open contours satisfy the condition $${\theta}_{0}=0$$, $${\theta}_{t0}/2{\omega}_{0}>1$$.

x0 = [0; 1.001*2*omega_0Value]; tFinal = 1.465*t_0; [t, x] = ode45(f, [tInit, tFinal], x0); t = t/t_0; x(:,1) = wrapToPi(x(:,1))/pi; x(:,2) = x(:,2)/(2*omega_0Value); figure; yyaxis left; plot(t, x(:,1), '-o'); ylabel('\theta/\pi'); yyaxis right; plot(t, x(:,2), '-o'); ylabel('\theta_t/2\omega_0'); grid on; title('Open Path in Phase Space'); xlabel('t/t_0');

This is how the pendulum motion would appear:

We have derived the simple pendulum's equation of motion, linearized and solved its equation of motion analytically, visualized its energy contours to understand the system qualitatively, and solved the general equation of motion numerically for particular initial conditions.

function lambda = wrapToPi(lambda) %wrapToPi Wrap angle in radians to [-pi pi] % % lambdaWrapped = wrapToPi(LAMBDA) wraps angles in LAMBDA, in radians, % to the interval [-pi pi] such that pi maps to pi and -pi maps to % -pi. (In general, odd, positive multiples of pi map to pi and odd, % negative multiples of pi map to -pi.) q = (lambda < -pi) | (pi < lambda); lambda(q) = wrapTo2Pi(lambda(q) + pi) - pi; end function lambda = wrapTo2Pi(lambda) %wrapTo2Pi Wrap angle in radians to [0 2*pi] % % lambdaWrapped = wrapTo2Pi(LAMBDA) wraps angles in LAMBDA, in radians, % to the interval [0 2*pi] such that zero maps to zero and 2*pi maps % to 2*pi. (In general, positive multiples of 2*pi map to 2*pi and % negative multiples of 2*pi map to zero.) positiveInput = (lambda > 0); lambda = mod(lambda, 2*pi); lambda((lambda == 0) & positiveInput) = 2*pi; end