During our high school days we are taught that a simple pendulum executes an
approximately simple harmonic motion if the angle of swing is small. However, high
school textbooks avoid discussing the general case: the motion of a
pendulum that may swing to larger angles. The main reason is that this leads
to an unmanageable differential equation that cannot be solved without a computer.
Consider the following diagram.
Simple pendulum
Taking $O$ as the origin and positive
$x$- $y$- and $\theta$-directions as shown, the position of
the bob is
$$\begin{eqnarray*}
x & = & L\sin(\theta)\\
y & = & -L\cos(\theta).
\end{eqnarray*}$$
Remember that $\theta $ is a function of time $t.$ So the above
equations actually mean
$$\begin{eqnarray*}
x(t) & = & L\sin(\theta(t))\\
y(t) & = & -L\cos(\theta(t)).
\end{eqnarray*}$$
The forces on the bob along the positive $x$- and $y$-directions
are, respectively,
$$\begin{eqnarray*}
F_x & = & -T\sin(\theta)\\
F_y & = & T\cos(\theta)-mg.
\end{eqnarray*}$$
Here $T$ is the tension in the rod. It is also a function of
$t.$
To derive the equations of motion we shall use Newton's second law of
motion, which says
$$\begin{eqnarray*}
F_x & = & m x''\\
F_y & = & m y'',
\end{eqnarray*}$$
where $x''$, $y''$ denote the second derivatives of $x(t)$, $y(t)$ with respect
to $t.$
Differentiating $x(t)$ and $y(t)$ twice we get
$$\begin{eqnarray*}
x'' & = & -L \sin(\theta)(\theta')^2 + L \cos(\theta)\theta ''\\
y'' & = & L \cos(\theta)(\theta')^2 + L \sin(\theta)\theta ''.
\end{eqnarray*}$$
Putting these in Newton's second law, and simplifying, we get
$$
\theta '' = -\frac gL \sin(\theta).
$$
At this point most textbooks use
the ``$\sin(\theta)\approx \theta $'' approximation for "small"
$\theta$ to reduce the above differential equation to
$$
\theta '' = -\frac gL \theta,
$$
which can be solved easily by hand to produce simple harmonic
motion. The approximation is pretty good if the pendulum swings
within $4^\circ$. But not all pendulums swing within that
range. What if you have a pendulum that swings $30^\circ?$
That's what we are going to explore now.
We first reduce the second
order differential equation $\theta '' = -\frac gL\sin\theta$
to a system of first order equations.
$$\begin{eqnarray*}
\theta ' & = & \omega\\
\omega' & = & -\frac gL \sin\theta.
\end{eqnarray*}$$
Notice that $(\theta',\omega')$ is given as a function
of $(\theta,\omega).$ The entire motion of the pendulum is
determined if we know $(\theta,\omega)$ at some instant. So
we call $(\theta,\omega)$ the phase of the system. We
are given the initial phase of the system, i.e., we know from which initial angle we have released the
pendulum, and with what angular velocity. Our aim is to know the phase at
all time points during the swing.
Thus, at $t=t_0$ (specified number), we know
$$\begin{eqnarray*}
\theta & = & \theta_0\mbox{ (specified number)},\\
\omega & = & \omega_0\mbox{ (specified number)}.
\end{eqnarray*}$$
We want to know the values $\theta(t)$ and $\omega(t)$ at any given $t > t_0.$
Thanks to the differential equations, we also know
the rate at which they are increasing at $t=t_0:$
$$\begin{eqnarray*}
\theta'(t_0) & = & \omega_0,\\
\omega'(t_0) & = & -\frac gL\sin \theta_0.
\end{eqnarray*}$$
Now advance time a little to $t_1=t_0+\delta t$, say. By this
time $\theta$ and $\omega$ will roughly change to
$$\begin{eqnarray*}
\theta_1 & = & \theta_0 + \theta'(t_0)\delta t = \theta_0+\omega_0\delta t\\
\omega_1 & = & \omega_0 +\omega'(t_0)\delta t =
\omega_0-\frac gL\sin \theta_0 \delta t.
\end{eqnarray*}$$
So we get the phase (approximately) at $t_1= t_0+\delta t.$
Now we keep on advancing time by $\delta t$ increments. The same
logic may be used repeatedly to give, at $t_k = t_0+k\cdot\delta t,$
$$\begin{eqnarray*}
\theta_k & = & \theta_{k-1} + \omega_{k-1}\delta t\\
\omega_k & = & \omega_{k-1} -\frac gL\sin \theta_{k-1} \delta t.
\end{eqnarray*}$$
Admittedly, this is a rather crude approximation. However, if $\delta t$ is pretty small, the accuracy increases.
Let's explore this numerically using R. First, we decide about a
time resolution, and the number of steps to run the process:
of time points:
dt = 0.1
n = 100
Next,
we create three arrays, one for $t$, one for $\theta$,
and the other
for $\omega$:
Yikes! What nonsense!
We want to play with the value of $\delta t$ (and the number
of steps). It is difficult to re-run all these lines of R
everytime we change the values. So at this point it is a good
idea to turn the above R code into a function:
to produce the same plot as before. Also the
variable result1 will store the numerical values in
case we need them.
Now let's increase the time resolution by setting $\delta
t=0.01$. To maintain the same time range don't forget to
change $n$ accordingly:
pendulum(0, 30* (pi/180), 1000, 0.01)
The result is
Much better
But still do you see something fishy?
EXERCISE:
Execute the above code with different initial values, and see if
the output changes as it should. Make a plot of the velocity over
time. Draw the phase diagram, i.e., a parametric plot
of $(\theta(t), \omega(t))$ with $t$ as the parameter.
EXERCISE:
Modify the code to allow the user to specify a non-zero initial velocity.
EXERCISE:
Plot the potential energy, kinetic energy and total mechanical
energy of the system as functions of time. Check if the total
mechanical energy curve is indeed a horizontal straight line, as
it should be.
EXERCISE:
Try to produce an animation of the pendulum in R. Hint: the function
The rod in the above pendulum is an inextensible one. So we could treat
$L$ as a constant. What if we replace it by a spring with constant
$\gamma?$ Then the tension in the spring will be
$$
T = \gamma\cdot \left(\sqrt{x^2+y^2}- \ell\right),
$$
where $\ell$ is the unstretched length of the spring.
Numerically solve these assuming that
$$
\ell=4,~~x(0) = 1,~~y(0)= -2,~~x'(0)=y'(0)=0.
$$
Animate to see if the solution looks natural. You may need to tweak $\delta t$ to make it look more natural (i.e.,
to make it more accurate).
Our differential equation was of the form
$$y'(t) = f(y),$$
where $y(t_0) = y_0.$ In our pendulum example we had $y
= (\theta,\omega).$ Also
$$
f(y) \equiv f(\theta, \omega) = \left(\omega, -\frac gL\sin \theta\right).
$$
The crude approximation that we used is called Euler's
method. It works with the more general form:
$$
y'(t) = f(t,y),\quad y(t_0) = y_0.
$$
We shall be given the function
$f(\cdot,\cdot)$ and the initial values $t_0,y_0.$ We are also given a
positive integer $n$ and a step size $\delta t.$ We have to find out the
function $y(t)$ at the points $t_1,...,t_n,$ where
$$
t_i = t_0+i \delta t.
$$
Euler's method works by making local linear
approximations to the unknown $y(t).$
For this we need to know the derivative of $y(t).$ If
at some instant $t$ we can guess the value of $y(t),$ then the value of $y'(t)$ may be obtained from
differential equation: $y'(t)=f(t,y(t)).$
Here we are starting from $y_0 =
y(t_0).$ So we know that $y'(t_0) = f(t_0,y_0).$ This is the slope
of the tangent to $y(t)$ at $t=t_0.$ Follow this tangent for
a little time $\delta t$ to arrive at $y_1 =
y_0+f(t_0,y_0)\delta t.$ The point
$(t_1,y_1)$ may not lie exactly on the curve of $y(t).$ But if $\delta t$
is small, then this should lie close to it. So we take $y_1$ as an
approximation to $y(t_1).$ Now we repeat the process again to get
$y_2 = y_1+\delta t\,f(t_1,y_1).$ In general, we define
$$
y_k = y_{k-1}+\delta t\, f(t_{k-1},y_{k-1}) \quad\mbox{ for } k=1,...,n.
$$
This is Euler's method.
How good is this method? To explore this we shall try out a simple example, where
the solution is known.
EXAMPLE:
Suppose we are working with $f(t,y) = -\sin(t+y).$ In other
words, we are solving $\frac{dy}{dt} = -\sin(t+y).$ We shall
start from the point $\left(0,-\frac \pi2\right),$ i.e., $y = -\frac \pi2$
when $t=0.$
Of course,
we can solve it analytically, by taking $v = t+y.$
Then $\frac{dv}{dt} = 1+\frac{dy}{dt}= 1-\sin(v),$ which may be solved by direct integration. The
answer is (check!)
$$
y = -\sin ^{-1} \left( \frac{1-t^2}{1+t^2} \right) - t.
$$
We
shall compare our approximation with this to
see how Euler's method performs. We shall take $n$ steps
in $[0,2].$ Then the Euler iterations are
$$
y_i = y_{i-1} - \frac 2n\times \sin(t_{i-1}+y_{i-1})
$$
for $i=1,...,n$ starting with $t_0 = 0$ and $y_0 =-\frac \pi2.$
The result (with $n=10$) is shown in the following graph. The continuous curve is the true solution. The dashed
polyline (with 10 segments) is Euler's approximation:
Euler with 10 steps
If we increase the number of steps to 20 then the approximation
is somewhat better:
Euler with 20 steps
If you use 100 steps the accuracy is pretty good:
Euler with 100 steps
EXERCISE:
Write an R function of the following form
euler = function(x0, y0, dx, n) {
...
}
to produce these plots. To draw multiple lines in the same plot
use the lines function:
Note that the range of $x$ and $y$-values are
determined by the first plot command. You might use
its xlim and ylim parameters to change
these ranges. See the help of plot function that:
?plot
Now try your hand at the following problem.
EXERCISE:
$\frac{dy}{dx} = e^{-xy^2}$ starting from $(0,0).$
One problem with Euler's method is that unless $\delta t$ is very small the
$y_i$'s may move away from the curve of $y(t).$ Taylor's
method generalises Euler's method to improve the accuracy. In Euler's
method we used local linear approximations (tangents)
to $y(t)$ at each $t_i:$
$$
y(t) \approx y(t_i) + y'(t_i)(t-t_i).
$$
These are just the first two terms of the Taylor expansion
of $y(t):$
$$
y(t_i) + y'(t_i)(t-t_i) + \frac{y''(t_i)}{2}(t-t_i)^2
+\cdots \frac{y^{(k)}(t_i)}{k!} (t-t_i)^k + \cdots
$$
In Taylor's method we take more terms from this series. Thus,
1-st order Taylor's method is the same as Euler's method, while
the $k$-th order Taylor's method uses
$$
y(t)\approx y(t_i) + y'(t_i)(t-t_i) + \frac{y''(t_i)}{2}(t-t_i)^2
+\cdots \frac{y^{(k)}(t_i)}{k!} (t-t_i)^k.
$$
EXERCISE: Check that it is the unique $\leq k$ degree polynomial that has the same
derivatives as $y(t)$ up to order $k$ at $t_i.$
EXAMPLE: Solve the same differential equation
$y'(t) = -\sin(t+y)$ with $y(0) = -\frac \pi2$
over $[0,2]$
using 2-nd order Taylor method.
SOLUTION:
The 2-nd degree Taylor polynomial for $y(t)$ at
any $t_{k-1}$ is
$$
y(t_{k-1}) + y'(t_{k-1})(t-t_{k-1}) + \frac{y''(t_{k-1})}{2} (t-t_{k-1})^2.
$$
For this we need $y'(t_{k-1})$ and $y''(t_{k-1}).$ These may be obtained approximately as follows.
We can write the differential equation as
$$
y'(t) = -\sin(t+y(t)).
$$
Differentiating w.r.t. $t,$
$$
y''(t) = -\cos(t+y(t))(1+y'(t)).
$$
Hence
$$\begin{eqnarray*}
y'(t_{k-1}) & = & -\sin(t_{k-1}+y(t_{k-1}))\\
y''(t_{k-1}) & = & -\cos(t_{k-1}+y(t_{k-1}))(1+y'(t_{k-1})).
\end{eqnarray*}$$
As before we take a grid of values in $[0,2],$ say $10$ steps. So
there are 11 points starting with $t_0=1$ and ending
with $t_{10}=2.$ The general formula for $t_k$ is
$$
t_k = 1+k\cdot\delta t,
$$
where $\delta t = \frac 2n.$
So the 2nd order Taylor iterations are
$$y_{k} = y_{k-1} +\frac 2ny'(t_{k-1}) + \frac{2}{n^2}y''(t_{k-1}).$$
The result is shown below. The red curve is the 2nd order Taylor approximation with $n=10.$ The dashed polyline is
the Euler's approximation (i.e., 1st order Taylor) with the same $n.$
Euler and 2nd order Taylor
($n=10$)
EXERCISE:
Implement the code in R to produce the above plot. Also try other
values of $n.$
This is a more complicated example.
Many science museums (including the Birla Industrial and
Technological Museum here in Kolkata) has a model to demonstrate
Einstein's theory of gravitation. The model consists of some
balls rolling on a
large curved plastic funnel. See
this Youtube video here. The funnel
represents the space-time warped by a heavy star (yellow ball) sitting at the
center. The smaller balls tend to roll into
the cavity, but owing to their initial tangential velocities end up
orbiting the star.
Consider the following funnel-like
surface. It is obtained by rotating the curve $z = f(y)$ around the
the $z$-axis. For instance, $f(y)=\sqrt{y-1}$ would produce a
surface like the following.
Ball in funnel
A ball is moving along the inner surface of the funnel. We shall ignore
the radius of the ball and the friction of the surface. (Thus the ball is a point mass slipping,
not rolling, on the funnel.) We know the
initial position and velocity of the ball. We want to find out the path
that the ball will follow.
There are two forces acting on the ball: its weight and the reaction of
the surface. The first works downwards, and so is
$$
\left[\begin{array}{ccccccccccc}0\\0\\-mg
\end{array}\right].
$$
The reaction acts inwards along the normal to the surface at the current
position of the ball. Let the current position of the ball be
$$
\left[\begin{array}{ccccccccccc}t\\y\\f(u)
\end{array}\right],
$$
where $u = \sqrt{x^2+y^2}.$ A little coordinate geometry shows that
the normal lies along
$$
\left[\begin{array}{ccccccccccc}-x\\-y\\u/f'(u)
\end{array}\right].
$$
If this is not clear, then think like this:
First consider the curve $z = f(y)$ and the point on the
graph at $y = u.$ The slope of the tangent there
is $f'(u)$ and hence the slope of the normal
is $-\frac{1}{f'(u)}.$ Understand how we get from this the normal arrow
shown below:
The arrow has coordinates $\left( 1,-\frac{1}{f'(u)} \right)$
Now when we rotate the curve around the $z$-axis to generate
the surface, the normal arrow also rotates with the curve. The
following diagram shows a snapshot during the rotation:
The arrow has coordinates $\left( \frac xu,
\frac yu, ,-\frac{1}{f'(u)} \right)$
The pink plane is the rotated version of the graph paper from the
earlier 2D graph of $z = f(y).$ If the new position
of $P$ is at $(x,y,0),$ then we must
have $\sqrt{x^2+y^2} = u.$ From this you should be able to
see that a vector along the (outward) normal direction is $\left( \frac xu,
\frac yu, ,-\frac{1}{f'(u)} \right)$, or equivalently, $\left( x,
y,-\frac{u}{f'(u)} \right)$. So the inward normal direction is given
by its negative: $\left( -x, -y, \frac{u}{f'(u)} \right)$
So the reaction force is
$$
R\left[\begin{array}{ccccccccccc}-x\\-y\\u/f'(u)
\end{array}\right],
$$
for some unknown function $R$ of $t$ (time). So we have the equation of
motion:
$$
m\left[\begin{array}{ccccccccccc}x''\\y''\\z''
\end{array}\right] = R\left[\begin{array}{ccccccccccc}-x\\-y\\u/f'(u)
\end{array}\right]+\left[\begin{array}{ccccccccccc}0\\0\\-mg
\end{array}\right].
$$
or, dividing by $m,$
$$
\left[\begin{array}{ccccccccccc}x''\\y''\\z''
\end{array}\right] = \tilde R\left[\begin{array}{ccccccccccc}-x\\-y\\u/f'(u)
\end{array}\right]+\left[\begin{array}{ccccccccccc}0\\0\\-g
\end{array}\right],
$$
where $\tilde R = \frac Rm.$
Notice
that $z$ is a known function of $x$ and $y.$
Find
$z''$ in terms of $x,y$ and their derivatives. Then use
the 3rd equation to get
$$
\tilde R = \frac{f'(u)\left(x'^2+y'^2-u'^2\right)/u+u'^2f''(u)+g}
{u\left(f'(u)+\frac{1}{f'(u)}\right)}.
$$
Thus, now $\tilde R$ is expressed in terms of $x$
and $y$ (and their derivatives).
With this $\tilde R$ we now have just two equations in two unknowns:
$$\begin{eqnarray*}
x'' & = & -x\tilde R\\
y'' & = & -y\tilde R.
\end{eqnarray*}$$
Use 2nd order Taylor method to solve this for the initial condition
$$
t(0)=10,~~y(0)=0,~~x'(0)=0,~~y'(0)=5.
$$
Take $g = 9.8.$
EXERCISE:
You can make the above problem more realistic by taking friction into
account. Remember that kinetic frictional force has magnitude proportional
to the normal reaction and acts opposite to the velocity vector.
Many other physics examples are
discussed in the web page www.myphysicslab.com. That page
has many interactive animations. However, they use a method more
sophisticated than what we have used. We
shall learn that method
in the second half of this course.
Comments
To post an anonymous comment, click on the "Name" field. This
will bring up an option saying "I'd rather post as a guest."
