Computer lab assignment 2022

Assessment No: 3
Assessment Type: Individual Computer Lab Assignment
Due Date: Friday, 11th October 2019
Student ID:
Name:

Introduction, Background and Project Overview
A dynamic problem can be split in two steps: obtaining the equation of motion and solving it. To obtain
the equation of motion we use common dynamics theory, such as Newton’s second law. The equation of
motion can then be solved through analytical or numerical techniques. So far, you have mostly looked at
analytical solutions to problems, however, even for simple dynamics problems these can rapidly become
complex and unwieldy. In actuality, the presence of a true analytical solution is rare in many real-world
problems and rely on simplifications and assumptions (e.g. neglecting drag). Numerical solutions on the
other hand can be readily used to solve simple and complex dynamics problems. For numerical solutions
the main challenges become computational resources and the elegance of the model.
Elastic pendulum (also called spring pendulum or swinging spring) is a physical system where a piece of
mass is connected by a spring. Compared with the ideal pendulum, not only angle, but also the length of
string (spring) is changing in the process, which makes the problem nonlinear, and extremely difficult to
solve analytically. An example in reality for elastic pendulum is the bungee jump. In this assignment, we
will first try to get the analytical solution for elastic pendulum with multiple approximations. Then finite
difference method will be adopted to get the numerical solution. In the end, we will try to build up a
simple numerical model for the whole bungee jump process.
Figure 1. The schematic graph for an elastic
pendulum. The grey curve is the trajectory of the
mass.
Table 1. Physical quantities in the system.
Property Symbol Quantity Unit
Item
Total mass m 0.1 kg
Initial angle θ0 0.05 rad
Initial angular
velocity ˙θ0 0.01 rad/s
Initial radius
velocity ˙r0 0.1 m/s
Spring
Free length l0 20 m
Initial elongation dr0 2 m
Elastic
coefficient k 20 N/m
For the following questions:
 This elastic pendulum is swinging in one plane, so two dimensional coordinate system is
adequate to describe the system.
 The equation of motion for a harmonic oscillator is in the form of :
¨θ ( t )+ a ⋅(θ ( t )−C)=0( a>0)
a and C are constants, which can vary in different equation of motion, and need to be found by
rearranging. Such form of equation has the analytical solution like:
θ ( t )= A cos ( ωn ⋅t−ϕ ) +C
Where
Angular frequency ωn= √ a
Amplitude A= √ (θ0 −C)2+¿ ¿

Phase ϕ =atan ¿
 Small angle theorem: when θ<5o , sinθ≈ θ, cosθ ≈ 1.

Equation of Motion
1. Draw the free body diagram for the elastic pendulum model in the Fig 1. Write down the
equations of motion in polar coordinate systems, and Cartesian coordinate system.
Answer within this box. Change the size of the box if needed.

Figure 2: The free body diagram
Polar coordinates ⃗
r =r ^et
v= dr
dt
¿ ˙r ^er +r ˙θ ^eθ
¿ vr ^er + vθ ^eθ ⃗ a= dv
dt
¿ ( ¨r −r ˙θ2 ) ^er + ( r ¨θ+2 ˙r ˙θ ) ^eθ +ar ^er + aθ ^eθ
vr { ( magnitude change ) ¨r
direction change ˙r ˙θ
vθ {magnitude change r ¨θ+ ˙r ˙θ
direction change r ˙θ2
(¨r −r ˙θ2) ^er (t) + a. (θ ( t )−C) =0
In Cartesian coordinates, the equation becomes
Tan θ = y
x
¨θ = 2 xy
( x2+ y2 ) 2
Therefore in Cartesian coordinates the equation becomes
2 xy
( x2+ y2 )2 + a.(atan y
x –c) = 0

Analytical Solution
2. To get the analytical solution, for a tough spring and a small mass, we can use the following
assumption:
 The elastic pendulum does harmonic oscillation in θ and r direction.
 In r direction, the oscillation of the spring is the dominate motion, and terms relevant with θ
can be ignored.
 In θ direction, the change of the spring length can be ignored, terms with r can be ignored.
Arrange the equation of motion in polar coordinate system to the harmonic oscillation form: (1)
use small angle theorem first, (2) then delete extra non-linear terms for approximation. Get the
analytical solution of elastic pendulum with the parameters in the table 1.
Answer within this box. Change the size of the box if needed.

⃗r =r ^et ,v= dr
dt
¿ ˙r ^er +r ˙θ ^eθ
¿ vr ^er + vθ ^eθ ⃗ a= dv
dt
¿ ( ¨r −r ˙θ2 ) ^er + ( r ¨θ+2 ˙r ˙θ ) ^eθ +ar ^er + aθ ^eθ
vr { ( magnitude change ) ¨r
direction change ˙r ˙θ
vθ {magnitude change r ¨θ+ ˙r ˙θ
direction change r ˙θ2
( ¨r −r ˙θ2) ^er (t) + a. (θ ( t ) −C) =0
θ ( t ) = A cos ( ωn ⋅t−ϕ ) …
Since it is known that A= √ (θ0 −C)2+ ¿ ¿, ωn= √ a and ϕ =atan ¿, the exact solution now becomes;
θ ( t )= √(θ0−C)2 +¿ ¿
replacing C with zero and Given that θ0 =0.005 and ˙θ0 = 0.1 rad/s
θ ( t )=¿
θ ( t )= √(0.05)2 +¿ ¿
∴
= ( √ 0.0025+ 0.01
a ) cos ( √ a⋅ t−atan ( 0.1
0.05∗√ a ))

Finite Differences
3. Write down the equation of motion in Cartesian coordinate system from Question 1 in the finite
difference form. State the initial conditions for finite difference method. Show all working.
Answer within this box. Change the size of the box if needed.

Expressing the equation in terms of x and y
xi+1 = x(ti+ Δt ¿
yi+1 = y( ti+ Δt ¿
The equation of motion can be rewritten as follows;
2 x (ti+ Δt ) y (ti+ Δt)
¿ ¿ ¿ + a.( tan−1 y (ti+ Δ t)
x(ti + Δt ) –c) = 0

MATLAB Code
4. Show the MATLAB code for solving the finite difference equations obtained in Question 3.
Answer within this box. Change the size of the box if needed.
(Include initialisation based on Question 3. Do not include the code part responsible for the plotting.)

val=2; %
delta_t=val/5; % the change in time
% the equation is as follows;
........a(arctan((y(tn+dt)/x(tn+dt))-c)+(2*x*(tn+dt)*y*(tn+dt))*...
...((x*(tn*+dt))^2+(y*(tn*+dt))^2)^2)^-2 +
...=0
c=val/2-1; % the excess term
a=val/2; % the initial value of constant
length=val; % the free length of the spring
delta_r=val/10; % the initial increment of the radius
for angle=0:36:180
y_axis=length*cosd(angle);
x_axis=length*sind(angle);
for t_i=0:1:5
eq_mot= a*(atan((y_axis*(t_i+delta_t))/(x_axis*(t_i+delta_t))-c))+...
((2*x_axis*(t_i+delta_t)*y_axis*(t_i+delta_t))*...
((x_axis*(t_i*+delta_t))^2+(y_axis*(t_i*+delta_t))^2).^-2) ; % the e
quation of motion
end
end

Discretization and Convergence:
5. The accuracy of numerical solution is completely dependent on the timestep dt. The accuracy and
the computational cost both increase indefinitely as dt → 0. Because of this, in practical
application of computational resources, we are only interested in dt which will provide results
that do not change/benefit from any further increase in resolution.
1) Estimate the dt value according to the frequency in Question 2. There should be at least 40
points in one period in each direction.
2) Plot the trajectory of the elastic pendulum from the numerical solution (x-y plot), with at least
2 periods in each direction.
Answer within this box. Change the size of the box if needed.

FDM
z=1;
% THE TERMS FROM THE TABLE AS FOLOWS
Length=2; % Spring length
mass=0.4; % bob mass
delta_t = 0.001; % small change in time
g=9.981; % force of gravity
% empty vectors
zeta= []; % the angle made when swinging
speed = []; % the velocity
period= []; % the time vector
% Initializing Vectors
speed(z)=5;
zeta(z)=0;
period(z)=0;
s(z)=0;
u=period(1);
while cos(zeta(z))-(speed(z))^2/(Length*g)>eps
z= z+1;
period(z)= u+(z-1)*delta_t;
speed(z) = speed(z-1)+delta_t*(g*(sin(zeta(z-1))-...
mass*cos(zeta(z-1)))+mass*(speed(z-1))^2/Length);
s(z)=s(z-1)+ delta_t*(speed(z)+speed(z-1))/2;
zeta(z)= zeta(z-1)+delta_t/...
(2*Length)*(speed(z)+speed(z-1));
end
s(z) = (s(z)+s(z-1))/2;
speed(z) = (speed(z)+speed(z-1))/2;
zeta(z)=((zeta(z)+zeta(z-1))/2);
s2 = Length*zeta(z);
zeta(z) = zeta(z)*180/pi;
disp(['Period = %5.4f sec\n',num2str((period(z)))])
disp(['angle = %5.4f deg\n',num2str((zeta(z)))])
%
plot(period,speed)

Numerical solution VS Analytical solution:
6. 1) Plot the x history (x-t plot) from numerical solution and analytical solution against each other
on the same figure. Are they similar with each other? How about y history (y-t plot)?
2) Increase the mass and decrease the elastic coefficient in Table 1 as you like. Plot the trajectory
from the numerical and analytical solution (x-y plot) with at least two periods in each direction.
Which solution is more authentic and close to reality? Why?
Answer within this box. Change the size of the box if needed.

val=2; %
delta_t=val/5; % the change in time
% the equation is as follows;
........a(arctan((y(tn+dt)/x(tn+dt))-c)+(2*x*(tn+dt)*y*(tn+dt))*...
...((x*(tn*+dt))^2+(y*(tn*+dt))^2)^2)^-2 +
...=0
c=val/2-1; % the excess term
a=val/2; % the initial value of constant
length=val; % the free length of the spring
delta_r=val/10; % the initial increment of the radius
z=1;
% THE TERMS FROM THE TABLE AS FOLOWS
Length=2; % Spring length
mass=0.4; % bob mass
delta_t = 0.001; % small change in time
g=9.981; % force of gravity
% empty vectors
zeta= []; % the angle made when swinging
speed = []; % the velocity
period= []; % the time vector
% Initializing Vectors
speed(z)=5;
zeta(z)=0;
period(z)=0;
s(z)=0;
u=period(1);
for angle=0:36:180
y_axis=length*cosd(angle);
x_axis=length*sind(angle);
for t_i=0:1:5
eq_mot= a*(atan((y_axis*(t_i+delta_t))/(x_axis*(t_i+delta_t))-c))+...
((2*x_axis*(t_i+delta_t)*y_axis*(t_i+delta_t))*...
((x_axis*(t_i*+delta_t))^2+(y_axis*(t_i*+delta_t))^2).^-2) ; % the e
quation of motion
end
end

while cos(zeta(z))-(speed(z))^2/(Length*g)>eps
z= z+1;
period(z)= u+(z-1)*delta_t;
speed(z) = speed(z-1)+delta_t*(g*(sin(zeta(z-1))-...
mass*cos(zeta(z-1)))+mass*(speed(z-1))^2/Length);
s(z)=s(z-1)+ delta_t*(speed(z)+speed(z-1))/2;
zeta(z)= zeta(z-1)+delta_t/...
(2*Length)*(speed(z)+speed(z-1));
end
s(z) = (s(z)+s(z-1))/2;
speed(z) = (speed(z)+speed(z-1))/2;
zeta(z)=((zeta(z)+zeta(z-1))/2);
s2 = Length*zeta(z);
zeta(z) = zeta(z)*180/pi;
disp(['Period = %5.4f sec\n',num2str((period(z)))])
disp(['angle = %5.4f deg\n',num2str((zeta(z)))])
%
plot(period,speed)

Numerical Solution Application:
7. Based on the finite difference code in Question 4, develop a new code for bungee jump
simulation. Use the following set ups to adapt the code into the bungee
system:
 When the distance between the person and the fixed end is less
than the free length of the cord, the person does free fall/ projectile
motion.
 When the distance between the person and the fixed end is more
than the free length of the cord, the bungee cord applies the tension
force and a damping force: FB=T + FD=−k∗( L(t )−L0 )−5∗v y
 To avoid entanglement of the cord, the distance between the two
ends of the cord is x 0=2 m at the beginning.
 The mass of the cord can be ignored. The mass of the person M =
60 kg.
The elastic coefficient k = 20 N/m. The free length of the cord is L0 = 20 m.
The initial velocity of the person is 0.
1) Plot the trajectory of the person in the Cartesian coordinate in the first 60 seconds.
2) Find the minimum distance to the ground and the cliff for the person to bungee jump safely.
3) Show the velocity magnitude and acceleration magnitude history of the person in the first 60
seconds.
Hint: use if command in Matlab to switch between two equations of motions.
Answer within this box. Change the size of the box if needed.

%Goal: Model the Spring Pendulum in 3-D with RK4
format 'long'
mass = 0.1; %Mass in kg
gravity = 9.81; %Gravity in m/s^2
k = 20; %Spring constant in N/m
Length = 2; %Unstretched length of spring in m
w = sqrt(k/mass); %Angular frequency of ossilations in rad/s
x_0 = 3; %Initial x-position in m
y_0 = 1; %Initial y-position in m
z_0 = -2.2; %Initial z-position in m
vx_0 = .1; %Initial x-velocity in m/s
vy_0 = .2; %Initial x-velocity in m/s
vz_0 = .3; %Initial x-velocity in m/s
delta_t = .001; %Time step in s
time = 20; %Total time to run simulation in s
%Vectors store positions and velocities
x_p = []; %Stores x-position
y_p = []; %Stores y-position
z_p = []; %Stores z-position
x_velocity = []; %Stores x-velocity

y_velocity = []; %Stores y-velocity
z_veocity = []; %Stores z-velocity
r_position = []; %Stores radial position
velocity = []; %Stores total velocity
for i = 1:time/delta_t
r = sqrt(x_0^2 + y_0^2 + z_0^2); %Radial position
%velocity in x direction
dvx1 = delta_t*(-w^2*(1-Length/r)*x_0);
dv_x2 = dvx1;
dv_x3 = dv_x2;
dv_x4 = dv_x3;
v_xdirection = vx_0 + (dvx1 + 2*dv_x2 + 2*dv_x3 + dv_x4)/6;
%Solves the y-velocity
dvy1 = delta_t*(-w^2*(1-Length/r)*y_0);
dvy2 = dvy1;
dv_y3 = dvy1;
dv_y4 = dvy1;
vy = vy_0 + (dvy1 + 2*dvy2 + 2*dv_y3 + dv_y4)/6;
%Solves the z-velocity

dvz1 = delta_t*(-w^2*(1-Length/r)*z_0 - gravity);
dvz2 = dvz1;
dvz3 = dvz1;
dvz4 = dvz1;
v_zdirection = vz_0 + (dvz1 + 2*dvz2 + 2*dvz3 + dvz4)/6;
%Solve the x-position
d_x1 = delta_t*vx_0;
dx2 = d_x1; dx3 = d_x1; dx4 = d_x1;
x_position = x_0 + (d_x1 + 2*dx2 + 2*dx3 + dx4)/6;
%Solve the y-position
d_y1 = delta_t*vy_0;
d_y2 = d_y1;
d_y3 = d_y1;
d_y4 = d_y1;
y_position = y_0 + (d_y1 + 2*d_y2 + 2*d_y3 + d_y4)/6;
%Solve the z-position
d_z1 = delta_t*vz_0;
dz2 = d_z1;dz3 = d_z1;dz4 = d_z1;
z_position = z_0 + (d_z1 + 2*dz2 + 2*dz3 + dz4)/6;

%Appends vectors to save all positions and velocities
x_p = [x_p,x_0]; y_p = [y_p,y_0];
z_p = [z_p,z_0];x_velocity = [x_velocity,vx_0];
y_velocity = [y_velocity,vy_0];z_veocity = [z_veocity,vz_0];
r_position = [r_position,r];velocity = [velocity,sqrt(vx_0^2 + vy_0^2 + vz_0^2)]
;
%Set the current positions and velocities as the original to iterate.
x_0 = x_position;y_0 = y_position;
z_0 = z_position;vx_0 = v_xdirection;
vy_0 = vy;vz_0 = v_zdirection;
end
%period vs position
plot(1:time/delta_t,r_position)
title('radial position against time')

Simple harmonic motion from a mass and spring system and that from a pendulum helps demonstrate
about natural frequency. Vibration is an important part of understanding the nature of materials and their
frequency.
The following equipment will contribute to the success of the lab exercise
i. Data studio and a computer
ii. Table clamp
iii. Pendulum bob
iv. Mass hunger
v. Stop switch
vi. Right angled rod clamp
vii. Spring
viii. meter stick
ix. balance scale
x. support rod
xi. short rode
xii. Set of mass
xiii. String
xiv. Stop watch
xv. Motion sensor
xvi. Force sensor

Procedure
Procedure
Simple pendulum
Mass of the bob was weighed and recorded. The bob was suspended from the bench by 1m string
Time t taken for bob to make a number of complete oscillations
Steps 1 and 2 were repeated using 0.8m , 0.6m, 0.4m and 0.2m.
Spring and mass
A small mass was hang from the spring. The masses were added progressively. In a table the masses were
recorded with the corresponding increase in height.A new mass was added and the spring allowed to settle. The
new height is recorded along with the weight.
The nest stage, data was collected using computer. Sensors were used to determine the height of the spring.
Complex harmonic motion is studied which is characterized by
a. The system obeys hook`s law. This is to mean that displacement of the mass is directly proportional to
displacement
Mv1=(m1+m2)v2
( m+ M ) g h2 +½ ¿
b. The position of the mass when oscillating is sinusoidal in nature. That is to mean that if the graph of time
against position is plotted, it will be a sine curve.
c. The last trait of the system is that energy is conserved. The energy is not wasted through friction or air
resistance.
There are some systems that do not meet one or more of these requirements and therefore are not classified as
simple harmonic motion.
Natural frequency is oscillating system is given by
ω= √ k
m
Period of the system is given by
T = 2 π √ k
m
Pendulum on the other hand has a restoring force that do not obey hook`s law, the gravity
F =-mg sin ϴ
When dealing with a small angle
F =-mg ϴ = -mg(x/l)
Mg/l is the effective k
T = 2 π √ l
g
Background
The
Procedure
Elastic pend
ulum
23
EGB211 Computer lab assignment

Mass of the bob was weighed and recorded. The mass of the steel ball was measured and recorded. The ball was
placed on the pendulum catcher and the distance between pivot and the balance pivot of the pendulum and the
ball.The pendulum was reconnected
The ball was shot such that it was trapped by the pendulum and the maximum angle ϴ3 recorded.
The set up is as shown in the figure below;
Figure 3: The elastic pendulum set up
The table below indicates the values of each of the components of the set up.
Property Symbol Quantity Unit
Item
Total mass m 0.1 kg
Initial angle θ0 0.05 Rad
Initial angular velocity ˙θ0 0.01 rad/s
Initial radius velocity ˙r0 0.1 m/s
Spring
Free length l0 20 m
Initial elongation dr0 2 m
Elastic coefficient k 20 N/m
The oscillators equation of motion is as follows
¨θ( t) + a ⋅ (θ (t) – C ) = 0 ( a > 0) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Equation 1
8. where a and C are constants.
9. the constants can be found by rearranging the equation. Such form of equation has the analytical solution like:
θ(t)=A cos (ωn⋅t – ϕ ) +C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .Equation 2
10. Where
11. Angular frequency wn is
ωn= √ a , . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .Equation 3
12. Amplitude A is given by
A= √(θ0 −C)2+ ¿ ¿ . . . . .. .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .Equation 4
13. And Phase angle is given by
ϕ =atan ¿ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Equation 5
14. when θ<5o , sinθ≈ θ, cosθ ≈ 1. This is known as the small angle theorem.
discussion
Oscillator has a periodic back and forth motion. Many systems in life can be approximated using oscillator model.
Many systems in nature display harmonic oscillation with natural or set frequency which depends on many factors in
the system.
24
EGB211 Computer lab assignment