Week 3 - Solving second order ODEs

Aim : To stimulate the motion of a simple pendulum with damping

Principle : The following ODE represent the eaquation of a motion of a simple pendulum with damping

where

b = damping coeffiecient

L = length of the pendulum

m = mass of the pendulum

g = 9.81m/s^2

since the above equation is in second order we have convert it into first order

assume

thetha = thetha1

dthetha1/dt = thetha2

hence the above equation becomes

dthetha2/dt + b/m thetha2 + g/l sin(thetha1)

the above equation in array form

d/dt[thetha1,thetha2] = [thetha2, -b/m thetha2 - g/lsin(thetah1)

asumming the time is 0 to 20

Angular displacement = 0

Angular velocity = 3 at t=0

Program 1

clear all
 close all
 clc

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

 theta_0 = [0;3];

 t_span = linspace(0,20,500);
 [t, results] = ode45(@(t,thetha) odefunct(t,thetha,b,g,l,m),t_span,theta_0)

 plot(t,results(:,1));
 hold on
 plot(t,results(:,2));
 xlabel('time');
 ylabel('plot'); 
hold off 

%for animation of the pendulum
ct =1 
 for i=1:length(results(:,1))
     x0=0
     y0=0
     x1=l*sin(results(i,1))
     y1=-l*cos(results(i,1))

     figure(2)
     plot([-0.7,0.7],[0,0],'linewidth',3,'color','b')
     hold on
     plot([x0 x1],[y0,y1],'linewidth',3,'color','r')
     hold on
     plot(x1,y1,'o','markers',18,'markerfacecolor','r')
     axis([-1 1 -2 0.5])
     pause(0.03)
     hold off
      M(ct)=getframe(gcf)
ct =ct+1;
 end

% for creating the movie
 movie(M)
 videofile = VideoWriter('pendulum_animation.avi','Uncompressed AVI')
     open(videofile)
     writeVideo(videofile,M)
     close(videofile)

Program 2

function [dthetha_dt] = odefunct(t, thetha, b,g,l,m)
thetha1 = thetha(1);
thetha2 = thetha(2);
dthetha1_dt = thetha2;
dthetha2_dt = [-(b/m)*thetha2-(g/l)*sin(thetha1)]
dthetha_dt = [dthetha1_dt;dthetha2_dt];
end

Explanation of the code

clear all remove the stored variables

close all remove the plots ,graph

clc clear all the values in command window

linspace command is since the pendulum moves through time 0 to 20sec

It generates 500 values btw 0 and 20

Inorder to solve the differential equation there is an inbulit solver called ode45. It is an standard solver for ordinary differential equations of fifth-order (ode45). The ode45 function applies Runge-Kutta formulae with the time step variable for easy computation.

Call function 'odefunct' is used to calculate [ t, results]

Plot command is used to plot the graph

For loop is used since time has an array of 500 values

since we have to fix one end of pendulum to a base

plot([-0.7,0.7],[0,0],'linewidth',3,'color','b')

plot([x0 x1],[y0,y1],'linewidth',3,'color','r')

first line creates fixed line from where which pendulum is fixed

second line creates the line of a pendulum

since pendulum has a mass at the end

plot(x1,y1,'o','markers',18,'markerfacecolor','r') is used

where 'o' is the markers

18 is the size

'r' is the color

To create a video file of the movie command movie(M) is used as shown in the program above

Error Faced :

All the frames are shown in the single graph. The graph looks cluttered. To rectify this Hold off command is used

OUTPUT

Graph

Output Screen

The link for the animation of the pendulum :https://youtu.be/xXGYWXXWiZ0