Week 3 - Solving second order ODEs

AIM :

To write a matlab program to solve ordinary differential equation of a simple pendulum and to simulate the motion of pendulum.

GOVERNING EQUATION :

Ordinary differential equation of a simple pendulum is as follows

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

To split this second order ODE into two first order ODE,

`theta_1=theta`

so the second order ODE is splitted into two first order ODEs as

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

`(d theta_2)/dt+b/mtheta_2+g/L sintheta_1=0`

so `(d theta_1)/dt=theta_2` and `(d theta_2)/dt=-b/mtheta_2-g/Lsintheta_1`

`d/(dt)[[theta_1],[theta_2]]=[[theta_2],[-b/mtheta_2-g/Lsintheta_1]]`

OBJECTIVES :

• To form a differential equation that represents angular position of bob of pendulum on each moment.
• To split second order differential equation into first order differential equation using a suitable function. This function is to be defined.
• To solve the ordinary differential equation of simple pendulum by using matlab for angular position and angular velocity with given boundary conditions.
• To plot the position of bob on each moment using angular positions obtained on solving ODE.

CONTENT :

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 with respect to time. The Ordinary differential equation is already developed which is to be solved by using Matlab. To make the problem easier we are splitting the second order differential equation into two first order differential equations using a suitable function. This function is to be defined as follows,

function [dtheta_dt]=ode_func(t,theta,b,g,l,m)
  %DEfining function with name ode_func to split second order ODE into two
  %first order ODEs
  %theta is column matrix with first element theta1 and next as theta2
  theta1=theta(1,1);
  theta2=theta(2,1);
  %theta2=dtheta1/dt
  dtheta1_dt=theta2;
  dtheta2_dt=-(b/m)*theta2-(g/l)*sin(theta1);
  %function return a column matrix of elemnts dtheta1/dt and dtheta2/dt
  dtheta_dt=[dtheta1_dt;dtheta2_dt];
end

Now to solve differential equation, we need to use the function ode45 in matlab. In ode45 we need to specify the function to solve (here it is ode_func), interval of integration vector (here time span of 20 seconds from 0 second) and the initial conditions. Here the results obtained is a column matrix theta which consist of Angular position and Angular velocity. So we need to specify the initial angular position and velocity.

Here the damping coefficient taken is 0.05 with which it will be difficult to identify the damping of pendulum. Thus the damping coefficient taken as 0.4 in the program.The main program will be as follows,

clear all
close all
clc
%defining variable related to pendulum
%define damping coefficient
b=0.4;
%define acceleration due to gravity
g=9.81;
%Define length of pendulum and mass of bob
l=1;
m=1;
%Fixing Initial condition, displacement and velocity at time=0
theta_0=[0;3];
%Fixing the interval for which the analysis to be done
t_span=linspace(0,20,200);
%Calling ode45 function to solve ode of Pendulum
[t,results]=ode45(@(t,theta) ode_func(t,theta,b,g,l,m),t_span,theta_0);
figure(1)
%plotting angular displacement against time
plot(t,results(:,1),'color','r')
hold on
%plotting angular velocity against time
plot(t,results(:,2),'color','g')
grid on
%labelling axes of plot
xlabel('Time (sec)')
ylabel('Ang Disp & Ang vel')
legend('Angular displacement','Angular velocity')
%array initiator to make array to store plots
ct=1;
for i=1:length(results(:,1))
    %taking angular displacement of pendulum
    the=results(i,1);
    %calculating x and y coordinated of bob during oscillation
    x=l*sin(the);
    y=-l*cos(the);
    %Plotting positions of bob
    figure(2)
    %ploting a line to look like a support, where pendulum is fixed
    plot([-0.6,0.6],[0,0],'linewidth',5,'color','b')
    hold on
    %Plotting bob by using O marker tools
    plot(x,y,'marker','o','markersize',20,'markerfacecolor','g')
    hold on
    %Plotting thread of pendulum by connecting position of bob and
    %supported poin
    plot([0,x],[0,y],'color','r','linewidth',3)
    hold off
    %Fixing Axis limits
    axis([-1.2 1.2 -1.5 0.1])
    %Fixing time interval between adjescent plots
    pause(0.03)
    %Storing each plots into an array
    M(ct)=getframe(gcf);
    ct=ct+1;
end
%creating a movie with array of plots
movie(M)
%creating video file
videofile=VideoWriter('Pendulum.avi','Uncompressed AVI')
%opening videofile
open(videofile)
%Assign matrix of plots to this videofile
writeVideo(videofile,M)
close(videofile)

RESULTS :

• The second order ODE of a pendulum is solved to obtain angular position and velocity.
• The variations of Angular position and Velocity are plotted against time.

• With the angular positions obtained the animation of movement of pendulum is created and it is uploaded in YouTube.

https://youtu.be/8sPRRrO6rn4

CONCLUSIONS :

Inorder to identify damping of pendulum it is better to increase the time span or increase the damping coefficient. To make video file of animation, movie2avi command was used in earlier versions of Matlab, it is to be noted for the old and low spec computers which use these earlier versions.