Week 3 - Solving second order ODEs

Simulating the simple pendulum motion using 2nd order ODE

Introduction to Simple pendulum motion: It consists of mass 'm' hanging by string of length 'l' from fixed support or fixed point. The mass 'm' called as pendulum bob ,when the bob is displaced to the inital angle and released, the pendulum bob moves to and fro in the form of periodic motion .It is closely related to the simple harmonic motion.

.

Where the motion of the pendulum depends on the Newton's second law for rotational systems and by following the law ,we obtain a second order Ordinary Differential Equation that describes the momentary position of the pendullum bob or (system) w.r.t time.

Aim: Write a program to simulate the pendulum motion using the given data and ODE equation.

Objective: To write a program that solves the below ODE. It represens the equation of motion of a simple pendulum with damping and to simulate the motion of the pendulum bob ,animate the motion of the pendulum by utilizing the given data below.

In the above equation,

where,

g = gravity in m/s2,

L = length of the pendulum in m,

m = mass of the bob in kg,

b=damping coefficient.

In-order to solve the second order Ordinary Differential Equation ,we have to convert the equation into simpler form of First order ODE by th following process,

Lets consider,

1. `theta=theta1` then
2. `d theta/dt =d (theta1)/dt`,
3. say `d (theta1)/dt=theta2` consider this as eq-1
4. Now ,let's say that `d^2theta/dt^2=d^2(theta1)/dt^2` ,as mentioned from the above point.1

5. so to simplify , `d^2(theta1)/dt^2=d/dt(d (theta1)/dt)=d/dt(theta2)` ,as mentioned from point .3 consider this as eq-2.

So ,Now coming to the second order ODE equation,

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

Now substituting the values of Eq-1 and point.4 in the ODE equation , we obtain the below equation

• `d^2(theta1)/dt^2+b/md(theta1)/dt+g/lsintheta1`Now when we furthur simplify this ODE equation into first order ODE ,by substituting the values of eq-1 and eq-2 ,

• `d/dttheta2+b/mtheta2+g/lsintheta1` , Hence a First order ODE equation had been aquired.

From the above second order ODE equation , we have obtained series of two first order ODE which is

`d(theta1)/dt=theta2` and

`d(theta 2)/dt=-b/mtheta2-g/lsintheta1`

Now represent the series of two first order ODE in the form of matrice ,so that to run in matlab and to get the output results

`d/dt[[theta1],[theta2]]=[[theta2],[-b/mtheta2-g/lsintheta1]]`

where the left hand side represents the `theta`(i.e `theta1,theta2`)and the right hand side of the matrice represents the 2 coloumns of solution vectors.

now we will solve the simulation of the pendulum motion with the help of above simplified series of two first order ODE and the input data values

L=1 metre,

m=1 kg,

b=0.05.

g=9.81 m/s2.

Simulate the motion between 0-20 sec, for angular displacement=0,angular velocity=3 rad/sec at time t=0.

Code of the project:

clear all
close all
clc

% inptuts for ode calculation and plotting
b=0.05;  %damping co-efficient
g=9.81; %acceleration due to gravity in 'm/s^2'
l=1;  %length of the pendulum string im 'm'
m=1;   %mass of the bob in 'kg'

%inital  boundary conditions of theta, where angular displacement=0,angular velocity=3 rad/sec at time t=0.
theta_0=[0,3];

%time points
t_span=linspace(0,30.02,200);

%solving ode equation
[t,results]=ode45(@(t,theta) ode_function(t,theta,b,g,m,l),t_span,theta_0);

%PLotting the graph for Angular displacement and Angular velocity w.r.t Time
figure(1);
plot(t,results(:,1));
hold on
plot(t,results(:,2));
ylabel('plot');
xlabel('time');
title('Angular displacement and Angular velocity w.r.t Time');
legend('Displacement Vs Time','Velocity Vs Time');
grid on

%inputs for simulating the pendulum
displacement=results(:,1);

ct=1;
for i=1:length(displacement)
    Disp=displacement(i);
    
   %co-ordinates 
    x0=0;
    y0=0;
    
    x1=l*sin(Disp);
    y1=-l*cos(Disp);
    
    %plotting
    figure(2);
    plot([-1 1],[0 0],'linewidth',3,'color','b')
    hold on
    plot([0 0],[0 -1],'linewidth',2,'color','g')
    hold on
    %plotting of the pin
    plot([x0], [y0],'.','markersize',30,'markerfacecolor','y')
    plot([x0 x1],[y0 y1] ,'linewidth',2,'color','m')
    axis([-1.5 1.5 -1.5 1])
    plot([x1],[y1],'o','markersize',20,'markerfacecolor','r')
    title('Simulation of Simple Pendulum')
    pause(0.00001)
    hold off
   M(ct)=getframe(gcf);
     ct=ct+1;
end
 movie(M);
 videofile=VideoWriter('second_order_ode.avi','uncompressed avi');
 open(videofile);
 writeVideo(videofile,M);
 close(videofile);





% function for ode45 formula
function [dtheta_dt]=ode_function(t,theta,b,g,m,l)
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​

Procedure:

1. Initially we enter the given input data which is required to calculate the ode equation

   % inptuts for ode calculation and plotting
   b=0.05;  %damping co-efficient
   g=9.81; %acceleration due to gravity in 'm/s^2'
   l=1;  %length of the pendulum string im 'm'
   m=1;   %mass of the bob in 'kg'​

   then , provide the initial boundary conditions and time points

   %initial  boundary conditions of theta, where angular displacement=0,angular velocity=3 rad/sec at time t=0.
   theta_0=[0,3];
   
   %time points
   t_span=linspace(0,20,200);

2. Now to solve the ode equation , we use the command called ode45 which is a pre-defined function in Matlab which integrates the initial boundary conditions w.r.t to time points by calling a function,where we have to define the function with required input data to solve the ode formula ,shown hereby

   % function for ode45 formula
   function [dtheta_dt]=ode_function(t,theta,b,g,m,l)
   theta1=theta(1);       
   theta2=theta(2);
   dtheta1_dt=theta2;   %the formula is breakdown into series of simple first order ODE 
   dtheta2_dt=-(b/m)*theta2-(g/l)*sin(theta1);
   dtheta_dt=[dtheta1_dt; dtheta2_dt];     %represented in matrice form
   end
   ​

   As the explanation behind the formula is expressed bove in terms of equation before entering the Code of the project

this Function is called to ode 45 in terms of the below format

%solving ode equation
[t,results]=ode45(@(t,theta) ode_function(t,theta,b,g,m,l),t_span,theta_0);

where we obtain results w.r.t time represented in the form [t,results] and then try to execute the programm where we obtain to 2coloumns of results where the 1st coloumn represents the Displacement and the 2d coloumn represents the velocity at span of 200 values using the linspace command.

3.Now to plot the graph between time and plots ,

PLotting the graph for Angular displacement and Angular velocity w.r.t Time
figure(1);
plot(t,results(:,1));   %using element wise operator to extract all the valuees of 1st coloumn
hold on
plot(t,results(:,2));     %using element wise operator to extract all the valuees of 2nd coloumn
ylabel('plot');
xlabel('time');
title('Angular displacement and Angular velocity w.r.t Time');
legend('Displacement Vs Time','Velocity Vs Time');
grid on

we obtain a plotted graph hereby,

In the plotted graph , the displacement point has not reached to the resting point at the inital position between the given timespan 0-20sec ,say t_span=linspace(0,20,200); so its not satisfying the periodic motion of a simple pendulum ,which can be shown below

AS you can observe that the bob has not reached to its resting positon at the inital point.

4.This error can be explained by taking two cases as examples and explain the factors affectiing,

Case 1 , where the initial boundary conditions of theta, where angular displacement=0,angular velocity=3 rad/sec at time t=0.

and ,increase the time span from 0-20sec to 0-30.02 sec where the displacement curve rests at initial point and meets at zero satisfying the condition of simple pendulum and we obtain the results as below:

Case2, let's increase the initial boundary conditions from the existing values , where angular displacement=0,angular velocity=5 rad/sec

say, theta_0=[0,5];

t_span=linspace(0,20,200); where time is between 0-20sec

then we obtain the following plot,

As you can observe that the displacement curve has rested at th point -1 on the axis, thus not fullfilling the condition of simple pendulum motion.

If we increase the timespan from 0-20 to 0-30sec , the result is as follows,

At this condition even it failed to meet at the resting point,as it require more timespan to reach the initial condition.

5.From both the examples , we came to know that increase in the boundary conditions results in increase in the time span and area of dispalcement also increases, and changes in the mass of the bob also reflects in the results of timespan .

6. In order to achieve the challenge , i had changed my time plots as follows

%time points
t_span=linspace(0,30.02,200);

And execute the program , we got the following plot,

7. Inorder to simulate the motion of a simple pendulum we introduce for loop,

the values of angular displacement are stored in a array called results(:,1)

Now assign these array of values to the new arrayfo called 'Displacement'

 %inputs for simulating the pendulum
displacement=results(:,1);

8.Input the for loop command ,where

for i=1:length(displacement) (where 'i' is called loop counter variable which can store single value at a time and '1:length(displacement)' this term is reffered as to run the loop from 0 to 200 values as theta array consists of 200values.Each and every time when the loop counter changes , it executes the commands inside the loop .)
Disp=displacement(i); (where 'Disp' variable stores single value for every loop counter variable changes , it assigns the value to the variable called "Disp")

Assign the coordinates to locate the position of the pendulum,as follows

ct=1;
for i=1:length(displacement)
    Disp=displacement(i);
    
   %co-ordinates 
    x0=0;
    y0=0;
    
    x1=l*sin(Disp);
    y1=-l*cos(Disp);
    

9.Now start plotting according to the coordinates by using plot command,

where 'hold on' command is used to hold all the figures after every counter loop run

and the plotting coordinates are provided for the support position and fixed pin ,length of the string and the bob.

the bob is plotted using command called marker and color of the bob is given as markerfacecolor, when run the command ,i came with an error

this error can be overcome by introducing a symbol representing the marker command as follows below

  %plotting
    figure(2);
    plot([-1 1],[0 0],'linewidth',3,'color','b')
    hold on
    plot([0 0],[0 -1],'linewidth',2,'color','g')
    hold on
    %plotting of the pin
    plot([x0], [y0],'.','markersize',30,'markerfacecolor','y')  %position of the fixed pin
    plot([x0 x1],[y0 y1] ,'linewidth',2,'color','m')             %position of the pendulum
    axis([-1.5 1.5 -1.5 1])
    plot([x1],[y1],'o','markersize',20,'markerfacecolor','r')      %position of the bob
    title('Simulation of Simple Pendulum')
    pause(0.00001)
    hold off
    M(ct)=getframe(gcf);
     ct=ct+1;

10. Now stich all the plots using the following commands,to animate the simulation of the simple pendulum

movie(M);
 videofile=VideoWriter('second_order_ode.avi','uncompressed avi');
 open(videofile);
 writeVideo(videofile,M);
 close(videofile);
 

Conclusion: Hence ,hereby to animate the motion of the pendulum using ODE function is achieved