Application Center - Maplesoft

App Preview:

Elliptic Functions and Simple Pendulum Animation

You can switch back to the summary page by clicking here.

Learn about Maple
Download Application


 

Image 

Animating A Simple Pendulum Using Elliptic Functions 

A.C. Baroudy
acbarou@inco.com.lb
 

Abstract 

The purpose of the present article is to built the equation of motion of a simple pendulum for large angles, then to use Maple commands: JacobiSN, to find the displacement angle for a given time.
A set of ordered pairs of coordinates for the bob is thus obtained and used for plotting.
The animation is somehow tricky since the set obtained doesn't allow for a straight forward animation and needs some careful manipulation of the data in the set.
The reward is double:
1- being able to put to good use Maple engine in finding Elliptic Functions (EF) with a simple command rather than resorting to tables for finding EF.
2- A good exercise in programming with Maple animation.
 

Equation of Motion of a Simple Pendulum with Large Angles 

From Figure-1 below we can see that
1- the
centripetal force F[n] just keeps the bob on a circular course.
   It is the result of the tension
T & the radial component `*`(mg, `*`(cos, `*`(theta))) of mg,
2- the
restoring force F[t] moves the bob along its circular path.
    It is the tangential component `+`(`-`(`*`(mg, `*`(sin, `*`(theta))))) of mg.
The
equation of motion is related to this force F[t]which is: 

                                           `and`(F[t] = ma[t], `and`(ma[t] = `*`(mr, `*`(diff(theta(t), t, t))), `*`(mr, `*`(diff(theta(t), t, t))) = `+`(`-`(`*`(mg, `*`(sin, `*`(theta)))))))rdiff(theta(t), t, t) = `+`(`-`(`*`(g, `*`(`sinθ`))))   
                                                                 
diff(theta(t), t, t) = `+`(`-`(`/`(`*`(g, `*`(`sinθ`)), `*`(r)))).                                       (1)
Note that the sign of the restoring force F[t]  
is just opposite of the sign of the angular displacement θ:
- when θ is on the positive side, F[t] is to the left
toward the equilibrium position and F[t] is negative,
- when θ is on the negative side, F[t] is to the right
toward the equilibrium position and F[t] is positive . 

Image
 

To solve the equation of motion (1) we use the following:
diff(theta(t), t, t)= `*`(`/`(`*`(d), `*`(dt)), `/`(`*`(`dθ`), `*`(dt))) = `/`(`*`(`*`(`/`(`*`(d), `*`(`dθ`)), `/`(`*`(`dθ`), `*`(dt))), `*`(`dθ`)), `*`(dt)) = `*`(`/`(`*`(d, `*`(diff(theta(t), t))), `*`(`dθ`)), diff(theta(t), t)) = `+`(`-`(`/`(`*`(g, `*`(`sinθ`)), `*`(r)))),
from which we get the differential equation for the motion of the pendulum:
                                                        `*`(`*`(d, `*`(diff(theta(t), t))), diff(theta(t), t)) = `+`(`-`(`/`(`*`(g, `*`(`sinθ`)), `*`(r))))`dθ` .                                             (2)
By integration we get : `+`(`*`(`/`(1, 2), `*`(`^`(diff(theta(t), t), 2)))) = `/`(`*`(g, `*`(cos, `*`(theta))), `*`(r)) + constant
i.e.                                                  `*`(`^`(diff(theta(t), t), 2)) = `+`(`/`(`*`(2, `*`(g, `*`(`cosθ`))), `*`(r)))  + constant.                                       (3)
If at
A (Figure-1) the initial value of θ is α then
When θ =
α , `*`(`^`(diff(theta(t), t), 2))= 0 and the constant = `+`(`-`(`/`(`*`(2, `*`(g, `*`(cos, `*`(alpha)))), `*`(r)))) ,
hence (3) becomes:                                       
                                                    `*`(`^`(diff(theta(t), t), 2)) = `+`(`/`(`*`(2, `*`(g, `*`(`+`(`cosθ`, `-`(`cosα`))))), `*`(r))).                                         (4)
Taking the square root of (4) we get :
                                              `/`(`*`(`dθ`), `*`(dt)) = sqrt(`+`(`/`(`*`(2, `*`(g, `*`(`+`(`cosθ`, `-`(`cosα`))))), `*`(r)))) .
i.e.                                          `/`(`*`(`dθ`), `*`(sqrt(`+`(`cosθ`, `-`(`cosα`))))) = `*`(sqrt(`+`(`/`(`*`(2, `*`(g)), `*`(r)))), `*`(dt)).                                  (5)
 

To make it simple we take r = g then by integrating (5) from θ = 0 to θ = α & t = 0 to t = t[alpha] we get the time from θ = α to θ = 0.
                                         = u.                                 (6)
For a complete period : from +α to -α and back, the time is 4T[alpha] .
To show that this integral is an
elliptic integral of the first kind we proceed as follows:
we use the identity :
`cosθ` = `+`(1, `-`(`*`(2, `*`(`^`(sin(`+`(`*`(`/`(1, 2), `*`(theta)))), 2))))) & `cosα` = `+`(1, `-`(`*`(2, `*`(`^`(sin(`+`(`*`(`/`(1, 2), `*`(alpha)))), 2))))) ,
and put :                                            
x = `/`(`*`(sin(`+`(`*`(`/`(1, 2), `*`(theta))))), `*`(sin(`+`(`*`(`/`(1, 2), `*`(alpha)))))) .                                                      (7)
Differentiating (7):
                                                      dx = `+`(`/`(`*`(`/`(1, 2), `*`(cos(`+`(`*`(`/`(1, 2), `*`(theta)))), `*`(`dθ`))), `*`(sin(`+`(`*`(`/`(1, 2), `*`(alpha))))))) ,
hence:
                                                       `dθ` = 2`/`(`*`(sin(`+`(`*`(`/`(1, 2), `*`(alpha)))), `*`(dx)), `*`(cos(`+`(`*`(`/`(1, 2), `*`(theta)))))) .                                                (8)
Equation (3) becomes:
= `/`(1, `*`(sqrt(2)))
`/`(1, 2) = =
replacing cos(`+`(`*`(`/`(1, 2), `*`(theta)))) with its equivalent sqrt(`+`(1, `-`(`*`(`^`(sin(`+`(`*`(`/`(1, 2), `*`(alpha)))), 2), `*`(`^`(x, 2)))))) from (7) we get:

,
                                                      =                                           (9)  
where we put `*`(`^`(k, 2)) = `*`(`^`(sin(`+`(`*`(`/`(1, 2), `*`(alpha)))), 2)).
 

In case α = π/2 then
                                                   k = `and`(sin(`+`(`*`(`/`(1, 2), `*`(alpha)))) = sin(`+`(`*`(`/`(1, 4), `*`(Pi)))), sin(`+`(`*`(`/`(1, 4), `*`(Pi)))) = `+`(`*`(`/`(1, 2), `*`(sqrt(2))))) ,                                  (10)
and
x being  `/`(`*`(sin(`+`(`*`(`/`(1, 2), `*`(theta))))), `*`(sin(`+`(`*`(`/`(1, 2), `*`(alpha)))))) while the integration at the very start was about θ and after integration we have to replace θ with α  then x = `*`(sin(`+`(`*`(`/`(1, 2), `*`(alpha)))), `/`(1, `*`(sin(`+`(`*`(`/`(1, 2), `*`(alpha))))))) = 1 in  case α = π/2 .
The integral in
x becomes:                         .                                  (11)
To put it in the
Jacobi form we set x = sin Φdx = `*`(cos, `*`(Phi, `*`(`dΦ`)))  and (11) becomes:
,
                                                   =                                                  (12)
                                                   = F(
sin Φ =1, k = `+`(`*`(`/`(1, 2), `*`(sqrt(2))))) = 1.854074677                              (13) 

Once we have the value of this integral i.e. the value of the time for descending from A to B ( or going upward from B to A) and which is equal to = 1.854074677 seconds, we divide it by 10 to get one small time increment : Δt. Then using Maple command JacobiSN(Δt,`+`(`*`(`/`(1, 2), `*`(sqrt(2))))) we get sin (ΔΦ) corresponding to Δt.
From relation (7) we get θ which defines our pendulum position at Δt we finally get θ:
sin (ΔΦ) = x  = `/`(`*`(sin(`+`(`*`(`/`(1, 2), `*`(theta))))), `*`(sin(`+`(`*`(`/`(1, 2), `*`(alpha)))))) = `*`(sin(`+`(`*`(`/`(1, 2), `*`(theta)))), `*`(`/`(`+`(`*`(`/`(1, 2), `*`(sqrt(2)))))))
                                                         → sin(`+`(`*`(`/`(1, 2), `*`(theta)))) = `+`(`*`(`/`(1, 2), `*`(sqrt(2), `*`(sin, `*`(`ΔΦ`))))).                                      (14)
This is angle θ of the pendulum as counted from
point B since ΔΦ is in fact from 0 to Φ = ΔΦ corresponding to time from t = 0 to t = Δt.
We repeat the same process for 2Δt to get Δ'Φ from which we get θ counted from
point B to Δ'Φ corresponding to 2Δt and so on till we reach the last one i.e. 10Δt = T[alpha = `+`(`*`(`/`(1, 2), `*`(Pi)))]= 1.854074677 corresponding to θ = `+`(`*`(`/`(1, 2), `*`(Pi)))always counted from point B.
                                                                               
Note

For the geometrical relations between the angles α,θ & Φ see :
1- Figure-2 below,

2- Greenhill's
"Elliptic Functions and Their Applications". 

Image 

 

 

Maple Program for plotting the pendulum bob 

As stated above we use Maple command JacobiSN to get sin (ΔΦ) from which we get angle θ. 

The following short program will give us just angle θ at successive time increments:
Δt,2Δt.......10Δt = 1.854074677. Then we get all x&y coordinates of the bob at these values of θ as a list.
To get the left half of the swing we just make the sign of all x-coordinates negative and set them in a 2d list.
Plotting is done for both lists with Maple command display.
 

> restart:with(plots):with(plottools):
 

> LIST1:=[[0,0]];LIST2:=[];k:=sqrt(2)/2;dt:=0.1854074677;
 

> for i to 10 do
sn:=JacobiSN(dt*i,(1/2)*sqrt(2));
theta:=2*arcsin(sn*k);
X:=9.8*evalf(sin(theta));Y:=9.8-9.8*evalf(cos(theta));
LIST1:=[op(LIST1),[X,Y]];LIST2:=[op(LIST2),[-X,Y]];
od:LIST1;LIST2;A:=pointplot(LIST1):B:=pointplot(LIST2):display([A,B],scaling = constrained,color=blue);
 

 

 

 

 

 

 





Plot_2d
 

Maple Program for animating the pendulum 

Here we have to pay attention to the way we have to animate the pendulum. Care should be exercised when animating the pendulum over the upswings and downswing and whether it is on the right half or the left half.
We need to have the pendulum start from point A then down to B and again up into the left half. From there down to B and again up to A. The coordinates are arranged in the lists from B up to A so we need to reverse the order of plotting. Furthermore the pair [0,0] should be suppressed from the 2d list so that no hanging is observed when the bob is passing this point from either direction.
The following procedure will take care of all these critical points
 

 

> `:=`(pendulum, proc (t0) local cir, t, ray, X1, Y1; if `and`(`<`(11, t0), `<`(t0, 22)) then `:=`(t, `+`(t0, `-`(11))); `:=`(X1, LIST2[t][1]); `:=`(Y1, LIST2[t][2]) end if; if `<`(t0, 12) then `:=`(t, ...
`:=`(pendulum, proc (t0) local cir, t, ray, X1, Y1; if `and`(`<`(11, t0), `<`(t0, 22)) then `:=`(t, `+`(t0, `-`(11))); `:=`(X1, LIST2[t][1]); `:=`(Y1, LIST2[t][2]) end if; if `<`(t0, 12) then `:=`(t, ...
`:=`(pendulum, proc (t0) local cir, t, ray, X1, Y1; if `and`(`<`(11, t0), `<`(t0, 22)) then `:=`(t, `+`(t0, `-`(11))); `:=`(X1, LIST2[t][1]); `:=`(Y1, LIST2[t][2]) end if; if `<`(t0, 12) then `:=`(t, ...
`:=`(pendulum, proc (t0) local cir, t, ray, X1, Y1; if `and`(`<`(11, t0), `<`(t0, 22)) then `:=`(t, `+`(t0, `-`(11))); `:=`(X1, LIST2[t][1]); `:=`(Y1, LIST2[t][2]) end if; if `<`(t0, 12) then `:=`(t, ...
`:=`(pendulum, proc (t0) local cir, t, ray, X1, Y1; if `and`(`<`(11, t0), `<`(t0, 22)) then `:=`(t, `+`(t0, `-`(11))); `:=`(X1, LIST2[t][1]); `:=`(Y1, LIST2[t][2]) end if; if `<`(t0, 12) then `:=`(t, ...
`:=`(pendulum, proc (t0) local cir, t, ray, X1, Y1; if `and`(`<`(11, t0), `<`(t0, 22)) then `:=`(t, `+`(t0, `-`(11))); `:=`(X1, LIST2[t][1]); `:=`(Y1, LIST2[t][2]) end if; if `<`(t0, 12) then `:=`(t, ...
`:=`(pendulum, proc (t0) local cir, t, ray, X1, Y1; if `and`(`<`(11, t0), `<`(t0, 22)) then `:=`(t, `+`(t0, `-`(11))); `:=`(X1, LIST2[t][1]); `:=`(Y1, LIST2[t][2]) end if; if `<`(t0, 12) then `:=`(t, ...
`:=`(pendulum, proc (t0) local cir, t, ray, X1, Y1; if `and`(`<`(11, t0), `<`(t0, 22)) then `:=`(t, `+`(t0, `-`(11))); `:=`(X1, LIST2[t][1]); `:=`(Y1, LIST2[t][2]) end if; if `<`(t0, 12) then `:=`(t, ...
`:=`(pendulum, proc (t0) local cir, t, ray, X1, Y1; if `and`(`<`(11, t0), `<`(t0, 22)) then `:=`(t, `+`(t0, `-`(11))); `:=`(X1, LIST2[t][1]); `:=`(Y1, LIST2[t][2]) end if; if `<`(t0, 12) then `:=`(t, ...
`:=`(pendulum, proc (t0) local cir, t, ray, X1, Y1; if `and`(`<`(11, t0), `<`(t0, 22)) then `:=`(t, `+`(t0, `-`(11))); `:=`(X1, LIST2[t][1]); `:=`(Y1, LIST2[t][2]) end if; if `<`(t0, 12) then `:=`(t, ...
`:=`(pendulum, proc (t0) local cir, t, ray, X1, Y1; if `and`(`<`(11, t0), `<`(t0, 22)) then `:=`(t, `+`(t0, `-`(11))); `:=`(X1, LIST2[t][1]); `:=`(Y1, LIST2[t][2]) end if; if `<`(t0, 12) then `:=`(t, ...
`:=`(pendulum, proc (t0) local cir, t, ray, X1, Y1; if `and`(`<`(11, t0), `<`(t0, 22)) then `:=`(t, `+`(t0, `-`(11))); `:=`(X1, LIST2[t][1]); `:=`(Y1, LIST2[t][2]) end if; if `<`(t0, 12) then `:=`(t, ...
`:=`(pendulum, proc (t0) local cir, t, ray, X1, Y1; if `and`(`<`(11, t0), `<`(t0, 22)) then `:=`(t, `+`(t0, `-`(11))); `:=`(X1, LIST2[t][1]); `:=`(Y1, LIST2[t][2]) end if; if `<`(t0, 12) then `:=`(t, ...
`:=`(pendulum, proc (t0) local cir, t, ray, X1, Y1; if `and`(`<`(11, t0), `<`(t0, 22)) then `:=`(t, `+`(t0, `-`(11))); `:=`(X1, LIST2[t][1]); `:=`(Y1, LIST2[t][2]) end if; if `<`(t0, 12) then `:=`(t, ...
`:=`(pendulum, proc (t0) local cir, t, ray, X1, Y1; if `and`(`<`(11, t0), `<`(t0, 22)) then `:=`(t, `+`(t0, `-`(11))); `:=`(X1, LIST2[t][1]); `:=`(Y1, LIST2[t][2]) end if; if `<`(t0, 12) then `:=`(t, ...
`:=`(pendulum, proc (t0) local cir, t, ray, X1, Y1; if `and`(`<`(11, t0), `<`(t0, 22)) then `:=`(t, `+`(t0, `-`(11))); `:=`(X1, LIST2[t][1]); `:=`(Y1, LIST2[t][2]) end if; if `<`(t0, 12) then `:=`(t, ...
`:=`(pendulum, proc (t0) local cir, t, ray, X1, Y1; if `and`(`<`(11, t0), `<`(t0, 22)) then `:=`(t, `+`(t0, `-`(11))); `:=`(X1, LIST2[t][1]); `:=`(Y1, LIST2[t][2]) end if; if `<`(t0, 12) then `:=`(t, ...
`:=`(pendulum, proc (t0) local cir, t, ray, X1, Y1; if `and`(`<`(11, t0), `<`(t0, 22)) then `:=`(t, `+`(t0, `-`(11))); `:=`(X1, LIST2[t][1]); `:=`(Y1, LIST2[t][2]) end if; if `<`(t0, 12) then `:=`(t, ...
`:=`(pendulum, proc (t0) local cir, t, ray, X1, Y1; if `and`(`<`(11, t0), `<`(t0, 22)) then `:=`(t, `+`(t0, `-`(11))); `:=`(X1, LIST2[t][1]); `:=`(Y1, LIST2[t][2]) end if; if `<`(t0, 12) then `:=`(t, ...
`:=`(pendulum, proc (t0) local cir, t, ray, X1, Y1; if `and`(`<`(11, t0), `<`(t0, 22)) then `:=`(t, `+`(t0, `-`(11))); `:=`(X1, LIST2[t][1]); `:=`(Y1, LIST2[t][2]) end if; if `<`(t0, 12) then `:=`(t, ...
`:=`(pendulum, proc (t0) local cir, t, ray, X1, Y1; if `and`(`<`(11, t0), `<`(t0, 22)) then `:=`(t, `+`(t0, `-`(11))); `:=`(X1, LIST2[t][1]); `:=`(Y1, LIST2[t][2]) end if; if `<`(t0, 12) then `:=`(t, ...
`:=`(pendulum, proc (t0) local cir, t, ray, X1, Y1; if `and`(`<`(11, t0), `<`(t0, 22)) then `:=`(t, `+`(t0, `-`(11))); `:=`(X1, LIST2[t][1]); `:=`(Y1, LIST2[t][2]) end if; if `<`(t0, 12) then `:=`(t, ...
`:=`(pendulum, proc (t0) local cir, t, ray, X1, Y1; if `and`(`<`(11, t0), `<`(t0, 22)) then `:=`(t, `+`(t0, `-`(11))); `:=`(X1, LIST2[t][1]); `:=`(Y1, LIST2[t][2]) end if; if `<`(t0, 12) then `:=`(t, ...
 

> `:=`(k_vals, seq(k, k = 1 .. 41)); -1; `:=`(to_animate, [seq(pendulum(k), k = k_vals)]); -1; display(to_animate, insequence = true, scaling = constrained); 1
`:=`(k_vals, seq(k, k = 1 .. 41)); -1; `:=`(to_animate, [seq(pendulum(k), k = k_vals)]); -1; display(to_animate, insequence = true, scaling = constrained); 1
`:=`(k_vals, seq(k, k = 1 .. 41)); -1; `:=`(to_animate, [seq(pendulum(k), k = k_vals)]); -1; display(to_animate, insequence = true, scaling = constrained); 1
 

Plot_2d
 

 

Legal Notice: ? Maplesoft, a division of Waterloo Maple Inc. 2008. Maplesoft and Maple are trademarks of Waterloo Maple Inc. Neither Maplesoft nor the authors are responsible for any errors contained within and are not liable for any damages resulting from the use of this material.  This application is intended for non-commercial, non-profit use only. Contact the authors for permission if you wish to use this application in for-profit activities.