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\theta$
to a system of first order equations.
$$\begin{eqnarray*}
\theta ' & = & \omega\\
\omega' & = & -\frac gL \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 the following J code:
^: (composition power): e.g.f^:3
means $f(f(f(\cdots))).$ Also, f^:(1 2 3)
means a list of functions, $f, f\circ f, f\circ f\circ f.$
i. creates a list of integers. In particular,
i. 5 is 0,1,2,3,4.
|: (transpose): matrix transposition.
plot when used as plot y makes a
line plot of $(0,y_0),...,(n-1,y_{n-1})$ if $n$ is the
length of $y.$
Here we are making 100 steps with $\delta t = 0.1.$ Try with $\delta t = 0.01$ also, and see how the plot changes.
The following R code implements the same idea:
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 phasor 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 (\sqrt{x^2+y^2}- \ell),
$$
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:
, joins two lists end to end. e.g. if $x=1, 2,
3$ then 5, x means 5,1,2,3.
_0.5p1 means $-\frac \pi2.$
(f g h) means the function $y\mapsto
g(f(y), h(y)),$ or $(x,y)\mapsto g(f(x,y), h(x,y)).$ Just as $(\sin+\cos)(x) = \sin x + \cos
x.$ In particular, 1 (-%+) x means $\frac{1-x}{1+x}.$
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.$
; makes boxes, e.g. 1; 'str' ; 3.9
will give a list of boxes.
;/ inserts a ; between consecutuve
elements of a list, e.g., ;/ 1 2 3 is the same as 1;2;3.
pd makes plots in a versatile way. Start
with pd 'reset' and end with pd
'show'. To plot $(x_i,y_i)$'s and join the points
with line segments use pd x; y
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{t^2+y^2}.$ A little coordinate geometry shows that
the normal lies along
$$
\left[\begin{array}{ccccccccccc}-t\\-y\\u/f'(u)
\end{array}\right].
$$
So the reaction force is
$$
R\left[\begin{array}{ccccccccccc}-t\\-y\\u/f'(u)
\end{array}\right],
$$
for some unknown function $R$ of $t.$ So we have the equation of
motion:
$$
m\left[\begin{array}{ccccccccccc}x''\\y''\\z''
\end{array}\right] = R\left[\begin{array}{ccccccccccc}-t\\-y\\u/f'(u)
\end{array}\right]+\left[\begin{array}{ccccccccccc}0\\0\\-mg
\end{array}\right].
$$
Thus, we have 3 equations in 3 unknowns: $t,y,$ and $R.$ Notice
that $z$ is a known function of $t$ and $y.$ To simplify
the equations first find
$z''$ in terms of $t,y$ and their derivatives. Then eliminate
$R$ to get two equations in two unknowns:
$$\begin{eqnarray*}
x'' & = & -t\tilde R\\
y'' & = & -y\tilde R,\\
\end{eqnarray*}$$
where
$$
\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)}.
$$
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."
