Week 3 - Solving second order ODEs

AIM: To write a program that solves the Ordinary Differential Equation (ODE).

OBJECTIVES: 1) To write a program to solve 2^nd order ODE which represents the equation of motion of a simple pendulum with damping.

2) To write a program that will simulate the pendulum motion between 0 to 20 sec to create an animated movie.

EQUATION

The primary equation use in this topic is as following;

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

Where,

g = gravity in m/

l = length of the pendulum in m,

m = mass of the ball in kg

b = damping coefficient

This is a 2^nd order ODE that is used in this topic to determine the pendulum motion.

THEORY :

A simple pendulum described ideally as a point mass suspended by a massless string from some point about which it is allowed to swing back and forth in a place. A simple pendulum can be approximated by a small metal sphere that has a small radius and a large mass when compared relatively to the length and mass of the light string from which it is suspended. If a pendulum is set in motion so that is swings back and forth, its motion will be periodic.

Fig: Simple Pendulum

Simple pendulum motion helps us in the following manner;

1) To study the simple harmonic motion

2)To learn the relationships between the period, frequency, amplitude, and length of a simple pendulum.

INPUT PARAMETERS :

Damping coefficient , b = 0.05

Acceleration of gravity, g = 9.81 m/

Length of the pendulum , l = 1 m

Mass attached , m = 1 kg

Velocity = 3 rad/s

PROGRAM CODE:

1) Main Program

Programming Language: MATLAB

clear all 
close all
clc

b = 0.05;
g = 9.81;
l = 1;
m = 1; 

%Initial condition
theta_0 = [0 3];

% time points
t_span = linspace(0,20,500);

 % solve ODE
 [t , results] = ode45(@(t,theta) ode_func(t,theta,b,g,l,m),t_span,theta_0);
 
 %plotting
 
 figure(1)
 
 plot(t,results(:,1))
 hold on
 plot(t,results(:,2))
 title('Simple Pendulum consisting Damping')
 legend('Angular displacement','Angular velocity in ras/s')
 ylabel('plot')
 xlabel('time')
 
 ct = 1;
 
 for i = 1:length(results(:,1))
     
     %Initial points
     x0 = 0;
     y0 = 0;
     
     
     %Final locations
     x1 = l*sin(results(i,1));
     y1 = -l*cos(results(i,1));
     
     %Plotting
     figure(2)
     
     plot([-1 1],[0 0],'linewidth',8,'color','b');  
     line([x0 x1],[y0 y1],'linewidth',5,'color','g')
     hold on 
     plot(x1,y1,'o','markersize',15,'markerfacecolor','r','color','r')
     axis([-2 2 -2 1]);
     title('Pendulum Movement')
     pause(0.03)
     hold off;
     
     M(ct) = getframe(gcf);
     ct = ct+1;
 end
 
 %Making Movie
 movie(M)
 videofile = VideoWriter('Pendulum.avi','Uncompressed avi')
 open(videofile)
 writeVideo(videofile,M)
 close(videofile)
     
     
     

2) Primary Function Program

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];

end

CODE EXPLANATION:

1. After the use of command ‘close all’, ‘clear all’, ‘clc’, all the input variables are defined: damping coefficient (b), gravity (g), length of the pendulum( l ), mass(m).

2. Here 2nd order ODE is converted into the first-order ODE, hence we provided two initial conditions, displacement, and velocity whose values are stored in an array called theta_0.

theta_0 = [0 3]

where,

displacement = 0 m

velocity = 3 rad/s

3. In order to solve ODE, we need time limit which is defined in the next step using linspace command whose values are stored in an array as follows

t_span = linspace (0,20,500)

here we are using 500 values in between 0 to 20 sec.

4. Before using the solver to solve ODE, a separate ode function called ode_func is created to accept the values of t, theta, b, g, l, m respectively. The file is saved with the same name ode_func to call it in the main program.

Also, this secondary program is used to split the higher-order into first order and created two variables namely dtheta1_dt & dtheta2_dt whose values are stored in single array dtheta_dt to create column matrix.

5.Finally to solve given 2^nd order ODE , we used default MATLAB solver ode45 in following format including all input parameters.

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

6.plot command is used to generate the graph for angular velocity and angular displacement. Also proper labels are given to x axis & y axis.

7. Now, ‘for loop’ is used to generate pendulum and to give the motion further to visualize the effect.

Here the condition is initialized by variable i & to start motion initial coordinates (x0, y0) as well as final coordinates (x1, y1) are provided to trace the swing area.

According to the variable i, pendulum takes various positions w.r.t to angle, velocity, time, etc.

8. Now figure(2) is used to generate a figure window.

9. To create a baseline, plot command is used & subsequently line command is used to create a pendulum arm. Then hold on command is used to keep the current plot while adding a new elements to the current plot.

10. To create a pendulum ball input ‘o’ is given for the marker shape of a circle. Limits are set for both the axes to give the required position of pendulum assembly. A pause of 0.03 is given.

11. To store individual frame ‘getframe’ function is used which will store in variable M. The counter is increased by 1 and end command is used to complete the for loop which terminates set of instructions.

12.movie(M) function is used to create a movie & a new uncompressed avi file is created to store it using VideoWriter command.

13. At last open, writeVideo &close functions are used to complete the video process for operating it.

OUTPUT:

1)Graph of effect of Damping

2)Pendulum motion animation:

Below is the link foe video;

https://drive.google.com/file/d/1K6eqkeIY9dzIDB_lla8bXcMTpD156dt3/view?usp=sharing

CONCLUSION:

The output for angular displacement & angular velocity is obtained for a simple pendulum. Also, animated simulation is created with the help of MATLAB program.