Solving Ordinary Differential Equation that represents equation of motion of a simple pendulum with damping.

AIM:

To write a Matlab program to solve an Ordinary Differential Equation (ODE) that represents equation of motion of a simple pendulum with damping.

SIMPLE PENDULUM:

The simple pendulum is an idealized mathematical model of a pendulum. This is a weight (or bob) on the end of a massless cord suspended from a pivot, without friction. When given an initial push, it will swing back and forth at a constant amplitude. Real pendulums are subject to friction and air drag, so the amplitude of their swings declines. A pendulum is a weight suspended from a pivot so that it can swing freely. When a pendulum is displaced sideways from its resting, equilibrium position, it is subject to a restoring force due to gravity that will accelerate it back toward the equilibrium position. When released, the restoring force acting on the pendulum's mass causes it to oscillate about the equilibrium position, swinging back and forth. The time for one complete cycle, a left swing and a right swing, is called the period. The period depends on the length of the pendulum and also to a slight degree on the amplitude, the width of the pendulum's swing.

From the first scientific investigations of the pendulum around 1602 by Galileo Galilei, the regular motion of pendulums was used for timekeeping, and was the world's most accurate timekeeping technology until the 1930s. The pendulum clock invented by Christian Huygens in 1658 became the world's standard timekeeper, used in homes and offices for 270 years, and achieved accuracy of about one second per year before it was superseded as a time standard by the quartz clock in the 1930s. Pendulums are also used in scientific instruments such as accelerometers and seismometers. Historically they were used as gravimeters to measure the acceleration of gravity in geophysical surveys, and even as a standard of length. The word "pendulum" is new Latin, from the Latin pendulus, meaning 'hanging'.

GOVERNING EQUATION:

Figure.1

1. The way the pendulum moves depends on the Newtons second law. When this law is written down, we get a second order Ordinary Differential Equation that describes the position of the "bob" w.r.t time.
2. The ODE that define the equation of motion of the pendulum is,

`(d^2theta)/dt^2 +(b/m)*(d theta)/dt+(g/l)*sintheta = 0` ................. (i)

Where,

`theta`= angular displacment`

g = gravity in m/`s^2`,

L = length of the pendulum in m,

m = mass of the bob in kg,

b=damping coefficient.

3. Now in order to solve the above 2^nd order ODE , we can break it into two 1^st order ODE so that it will easy for us to write the code in Matlab.

So let,

`theta = theta_1` .................(ii)

Therefore,

`(d theta)/dt=(d theta_1)/dt`

Now let,

`(d theta_1)/dt=theta_2` ................. (iii)

Therefore,

`(d^2 theta)/dt^2=(d^2 theta_1)/dt^2=d/dt((d theta_1)/dt)=(d theta_2)/dt` .............(iv)

Now by substituting equation (ii), (iii), (iv) in equation (i), we will get,

`(d^2theta_1)/dt^2 +(b/m)*(d theta_1)/dt+(g/l)*sintheta_1 =0`

`(d theta_2)/dt+(b/m)*theta_2 +(g/l)*sin theta_1 =0`

`(d theta_2)/dt=-(b/m)*theta_2 -(g/l)*sin theta_1` .................... (v)

4. Now equation(iii) and equaton(v) are two first order equation and this is how we can break 2^nd order ODE into series of 1^st order ODE.
5. Now this two ODE’s are represented as matrix form so that it can be entered in Matlab easily.

`d/dt [(theta_1), (theta_2)] = [(theta_2), ((-b/m)*theta2-(g/l)*sintheta_1)]`

6. Now to find the co-ordinates of the two ends of pendulum we have used trigonometric equations.
7. Since one end of the pendulum is fixed to a fixed surface, its co-ordinate values are also fixed and they are considered as ‘(`X_0,Y_0`)’ which is equal to (0,0).
8. The co-ordinates of the other end of the pendulum which is position of the bob are considered as (`X_1,Y_1`) and they are calculated using following equations,

`X_1 =l*sin(theta)` …………… (vi)

`Y_1=-l*cos(theta)` ……………. (vii)

[ here ‘ – ‘ is used since the pendulum is below the x-axis according the assumptions]

Where,

`l` = length of string/rod

`theta` = angular displacement

9. Using above equations and matrix we are going to solve the ODE in Matlab.

OBJECTIVE:

1. To solve the ODE and find the values of angular displacement (i.e. `theta`) and angular velocity at each and every time interval‘t’ of the pendulum.
2. To create an animation video representing the motion of the pendulum.

PROGRAM:

CODE FOR SOLVING ODE:

close all
clear all
clc


%Inputs
b = 0.05;
g = 9.81;%(in m/s^2)
l = 1;%( in metres)
m = 1;%( in kilograms)

%Initial conditions
theta_0=[0;3]; %( in radians and rad/sec)

%Time points
t_span = linspace(0,20,500);%( In Seconds)

%Solve ODE
[t , results]= ode45(@(t,theta) ode_func(t, theta, b,g,m,l) ,t_span, theta_0);
 
%Plotting
plot(t,results(:,1)) %(GRAPH FOR DISPLACEMENT CURVE W.R.T TIME)
hold on
plot(t,results(:,2)) %(GRAPH FOR VELOCITY CURVE W.R.T TIME)
xlabel('Plot')
ylabel('Time')
legend('Dispalcement','Velocity')

EXPLANATION:

1. First the code is started with ‘close all’, ‘clear all’, ‘clc’ to close any previous plots, to clear the workspace and command window respectively.
2. Then the inputs are entered i.e. damping coefficient (b), acceleration due to gravity (g), length of the pendulum (l), mass of the bob (m).
3. Now here initial conditions i.e. initial angular displacement and initial angular velocity of the pendulum are given as inputs and the values are assigned to a variable ‘theta_0’ and values are given as a column vector.
4. The reason for entering the initial values in a column vector is that , as defined above 2^nd order ODE is broke down into two 1^st order ODE’s and the first ODE(i.e. equation (iii) ) is the angular displacement equation and the second ODE (i.e. equation (v) ) is the angular velocity equation.
5. In Matlab since these two equations are entered as column matrix so the initial values are also entered in a column matrix to equate the initial values of the angular displacement and angular velocity of the pendulum.
6. Now the time interval upto which the motion of the pendulum has to be calculated is defined by using ‘linspace()’ command. Here the time interval is taken from 0-20 seconds with 500 values between them. This means that we are going to get 500 values of angular displacements and angular velocities in the result for 500 values of time intervals.
7. Now to solve the ODE mentioned above, we have used a command called ‘ode45()’ which is in-built Matlab solver. Syntax for this command is

[tSol , ySol] = ode45(@(t,y) function_name ,[t0, t_end], y0)

Here,

tSol , ySol = Solution vectors, Matlab returns ySol for each time tsol.

function_name = Function that returns dy= F(t,y).

(t0,y0) = initial conditions representing.

t_end = final value until which solution is desired.

8. Here in the above syntax ,if an equation is dy/dt=F(t,y) , ‘t’ , ‘y’ are the dependent variable which return the value of dy/dt, so in ‘ode45’ the dependent variables are first enter as’, ‘@(t,y)’, here ‘@’ defines that we don’t know the value of ‘t’ , ‘y’ and that ‘t’, ‘y’ are input to the anonymous function F(t,y).
9. Then this anonymous F(t,y) is defined according to the problem. This is equivalent to defining a function that takes the input of ‘t’ , ‘y’ and gives the output dy/dt.
10. Now the function here can be entered directly or can be defined separately and then can be used in ‘ode45’ command.
11. Next the initial and final value of ‘t’ are to be entered for which you desire solution. So the input for ‘t’ is taken between ‘t0’ and ‘t_end’.
12. Initial value of ‘y’ is entered (i.e. y0) from where the solution will start and rest solution is calculated thereafter.
13. In the above code solution vector is [t , results] and since the ODE defined in our problem is dependent on ‘t’ and ‘θ’ we have enter as ‘@(t,theta)' .
14. And then a function is defined separately with the inputs and then it is used in ‘ode45’ command.

FUNCTION CODE:

function [dtheta_dt] = ode_func(t,theta,b,g,l,m)
 theta1 =theta(1);
 theta2 = theta(2);
 dtheta1_dt=theta2;
 dtheta2_dt= -(b/m)*theta2 - (g/l)*sin(theta1);
 dtheta_dt = [dtheta1_dt;dtheta2_dt];

15. In the above function entered we have used the equations (ii), (iii), (iv), (v) and the matrix which was derived above.
16. Here the fuction name is given as ‘ode_func’ and the inputs for this function are given as time (t), angular displacement (`theta`), damping coefficient (b), acceleration due to gravity (g), length of the pendulum (l), mass of the bob (m).
17. The left-hand side of the matrix equation,represents `theta` which is broken into `theta_1` and `theta_2` so here 'theta(1)' means that in the `theta`matrix first value is `theta_1`which is defined as 'theta1'.
18. Similarly for 'theta(2)' it is defined as 'theta2'.
19. Now this function is saved in the same folder where the main code is saved and then this function is called in the 'ode45' command.
20. Since initial and final time is defined as ‘time_span’ so it is entered directly in the ‘ode45’ command, and similarly initial conditions are defined in the variable ‘theta_0’ it is also entered directly.
21. Now when we execute the code we will get 500 values for ‘t’ and a (500x2) matrix for ‘results’ since ‘t_span’ is a row vector and ‘theta_0’ is a column vector as solution we have got a (500x2) matrix.

Figure.2

22. We can observe in Figure.2 , that in the Workspace window for ‘results’ we have (500x2) size matrix, and for ‘t’ we have (500x1) size matrix.
23. In the results the first column represents the values of angular displacements and second column represents values of angular velocities’
24. And to obtain a graph we have used ‘plot()’ command. First ‘plot()’ is used for generating a curve that represents how angular displacement changes with time and second ‘plot()’ command is used to generate a curve that represents how angular velocity changes with time.
25. A ‘hold on’ command is used to avoid the overwriting of the ‘plot()’ command.
26. Then 'legend()' command is used to show which curve represents angular displacement and angular velocity in the graph.

Figure.3

26. From above figure we can observe that the angular displacement and angular velocity are decreasing as the time increases that means the pendulum is coming back to rest.

CODE FOR ANIMATION:

close all
clear all
clc


%Inputs
b = 0.05;
g = 9.81; %(in m/s^2)
l = 1; %( in metres)
m = 1; %( in kilograms)

%Initial Conditions
theta_0=[0;3];

%Time Points( Initial and final time)
t_span = linspace(0,20,500);%( In Seconds)

%Solving ODE
[t , results]= ode45(@(t,theta) ode_func(t, theta, b,g,m,l) ,t_span, theta_0);

%Plotting
figure(1)
plot(t,results(:,1));
hold on
plot(t,results(:,2));
xlabel('Plot');
ylabel('Time');
legend('Dispalcement','Velocity');

%Animation
displacement= results(:,1);
ct=1;
for i=1:length(displacement)
    x0=0;
    y0=0;
    x1=l*sin(displacement(i));  %(X co-ordinate of the mass bob)
    y1=-l*cos(displacement(i)); %(Y co-ordinate of the mass bob)
    figure(2)
    plot([-1 1],[0 0] ,'linewidth',5)  %(Co-ordinates of the fixed support)
    hold on
    plot([x0 x1],[y0 y1],'linewidth',3); %(Joining fixed end and the mass bob)
    hold on
    plot(x1,y1,'o','MarkerSize',10,'MarkerFaceColor','r') %(For representing the bob)
    hold off
    axis([-1.5 1.5 -1.5 1.5])
    pause(0.001)
    M(ct)=getframe(gcf)
    ct=ct+1;
end


 movie(M)                      %(For creating animation video)
 videofile= VideoWriter('Simple_pendulum.avi','Uncompressed AVI');
 open(videofile)              %(Opening the vidoefile while running the code)
 writeVideo(videofile,M)      %(For writing the videofile)
 close(videofile)             %(Closing the videofile after running the code)

EXPLANATION:

1. To create an animation video we have used different positions of the pendulum at each time interval.
2. Position of the pendulum is determined by the values of angular displecment.
3. Since first column of the ‘results’ matrix are the values of displacement ,they are assigned to a variable ‘displacement’ .
4. Now we have used ‘FOR LOOP’ here to determine the different positions of the pendulum.
5. In ‘FOR LOOP’ first the it is defined that ‘i=1:length(displacement)’, which means that ‘i’ will vary from first value of angular displacement to last value of angular displacement. Here ‘i’ is a loop counter variable which takes a new value every time it is repeated in FOR LOOP.
6. Next fixed end point of the pendulum is defined i.e. ‘x0’, ‘y0’, since they are fixed ,values are also fixed.
7. Then the other end points of the pendulum are defined using the equation(vi) and equation(vii) as mentioned above.
8. Figure() command is used so that the animation is in different window and displacement and velocity graph is in different window.

ERROR OCCURED: Here first time I have used ‘figure(2)’ command outside the FOR LOOP, which lead to overwriting of figure(2) on figure(1), below picture represents the error. This is because if figure(2) is written outside FOR LOOP then for the first iteration it will consider it and then from the second iteration Matlab will consider only ‘figure(1)’ and it will overwrite it. So if you write it inside the FOR LOOP then Matlab will consider figure(2) only for every iteration.

Figure.4

You can observe that both the figure show same in both the figure windows.

9. Now to create a fixed support of the pendulum we have used ‘plot()’ command and defined points to get a horizontal fixed support.
10. Then ‘hold on’ command is used to avoid overwriting ‘plot()’ command.
11. Next ‘plot()’ command is used to join the fixed end and the position of the bob and then again ‘hold on’ command is used.
12. Next ‘plot()’ command is used to represent the bob like a circular mark in the figure , ‘MarkerSize’ means the size of the bob to be represented and ‘MarkerFaceColor’ means the color of bob.
13. Then ‘hold off’ command is used as for the next iteration, the first plot command has to executed again.
14. Then ‘axis()’ command is used to define the axis of the graph and ‘pause()’ command is used to pause the each position of the pendulum for specified time to get an animation.
15. Counter variable (ct) is defined so that iteration is goes on in FOR LOOP.
16. Each position is captured as frames using ‘getframe(gcf)’ command and it is assigned to a variable ‘M(ct)’ which creates a array of frames. Here (ct) is number of frames.
17. Every time when loop is executed current figure is saved in the array ‘M’. Then the FOR LOOP is ended .
18. After the end of the FOR LOOP we have entered command ‘movie(M)’. This command instructs the MATLAB to generate a movie using all the frames in ‘M’.
19. Now to create a AVI file, where MATLAB will save this video, we have used command ‘VideoWriter()’ where we have to enter file name and profile, in this command. In the above code it is specified accordingly. This command is assigned to a variable ‘videofile’.
20. Now to open this video while code is running,’ open(videofile)’ command is used ,since ‘videofile’ is the variable used so it is used inside the brackets.
21. Then, we are going to write in the ‘videofile’ , the frames that we captured using ‘writeVideo(videofile,M)’ command, where the ‘videofile’ and ‘M’ are assigned variables.
22. Now to close this file after the code has finished running, we have used ‘close(videofile)’ command.
23. Now when we run the code Matlab will create a video and it will save in the same folder where the code is saved.
24. Generated video link is shared for reference.

https://youtu.be/Yb4lFF96q2E

CONCLUSIONS:

1. Thus we have solved the ODE that represents equation of motion of a simple pendulum and created a animated video of the pendulum successfully.
2. We can solve any kind of ODE and ODE related problems using this approach, for example we can solve a third order ODE by breaking it into three first order ODEs and writing similar Matlab code to solve it.