Week 3 - Solving second order ODEs

Title:- Program in Matlab that will simulate the pendulum motion, use ode45 to solve the 2nd order differential equation for pendulum.

Aim:- Simulate the pendulum motion.

ODE45 :- Ode 45 is a inbuilt function present in matlab which generaly used to solve 1 order equation.

General syntax for calling ode45 is.

[t , theta] = ode45(@(t, theta) function(variables1, variables2), intigration limit, initial condition)

Objective:- write a program that solves the pendulum ode

Code:-

 %% program that solves the following ODE
 
 % given ode d^2(theta)/dt^2 +b/m*(dtheta/dt)+g/l*sind(theta) = 0
 
 clear all 
 close all
 clc
 
 %% given variables
 
l =1; % metre,
m=1; % kg,
b=0.05;
g=9.81; % m/s2

%% given initial condition

t_span =linspace(0,20,200);
theta_ini = [0;3];

%% solving ode

[t ,theta] = ode45(@(t, theta) ode_pendulum(t,theta,b,g,l,m), t_span,theta_ini);

angular_displacement = theta(:,1);
angular_velocity = theta(:,2);

%% plotting time  vs angular velocity &  time vs angular displacement
figure(1)
plot(t,angular_displacement);
hold on
plot(t,angular_velocity);
hold off
axis([0 10 -4 4])
legend('Angular Displacement','Angular Velocity')

%% plotting the pendulum on graph
p = 1;
for i = 1:length(angular_displacement)
    theta_displacement = angular_displacement(i);
    
    x0 = 0;
    y0 = 0; 
    
    x1 = l*sin(theta_displacement);
    y1 = -l*cos(theta_displacement);
    
    figure(2)
    plot([-1 1],[ 0 0],'linewidth',7)
    hold on
    plot([x0 x1],[y0 y1], 'linewidth' , 3)
    plot(x0,y0,'O','markersize',3,'markerfacecolor','y')
    plot(x1,y1,'O','markersize',20,'markerfacecolor','r')
    xlim([-2 2]);
    ylim([-1 1.5]);
    hold off
    pause(.001)
    M(p) = getframe(gcf);
    p=p+1;
   
end

%% Generating movie

movie(M)
videofile = VideoWriter('pendulum.avi;','Uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)

Function Code:-

function dtheta_dt = ode_pendulum(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 

Result:-

As result we are geeting time vs angular displacement & time vs Angular velocity graph.

We are getting an animation of pendulum.

Youtube Link :- https://youtu.be/1BAFakY-1MA

Explanation:- given Variables.

l =1; % metre,
m=1; % kg,
b=0.05;
g=9.81; % m/s2

given pendulum equation is 2nd oder equation in order to solve 2nd order ode we need to convert 2nd order ode into 1st order ode.

given equation

d^2(theta)/dt^2 +b/m*(dtheta/dt)+g/l*sin(theta) = 0 --------1st equation

let, theta = theta1

dtheta/dt = dtheta1/dt=theta2 --------- 2nd equation

d^2theta/dt^2 = d^2theta1/dt^2 = dtheta2/dt ---------3rd equation

put euqation 2 and 3 in equation 1 we get,

dtheta2/dt + b/m*theta2 + g/l*sin(theta1) = 0

dtheta2/dt = -b/m*theta2 -g/l*sind(theta1)

for inserting the intergration limit between 0 to 20 s i used linspace(0,20,150);, i use 150 values in between

for giving initial consition of theta, i use theta_ini =[0;3]

1st coloum of theta_ini represent angular displacement

2nd coloum of theta_ini represent angular velocity

i created a function which solves the ode

function dtheta_dt = ode_pendulum (t, theta, b,g,l,m), everything written is brakets are input variables.

then i define the the funtion as

theta1 = theta(1)

theta2 = theta(2)

dtheta1_dt = theta2;

dtheta2_dt = -(b/m)*theta2 -(g/l)*sin(theta1)

dtheta_dt = [dtheta1_dt;dtheta2_dt] is the output which store values

now by using end command i end the function and saved this file as 'ode_pendulum.m'

its mandatory that dunction file name should mataches with function name

in order to use ode45 which is inbuilt function.

1st we need to decide the output then we have to call the generated funtion inside the body of ode45

to do this,

[t, theta] = ode45(@(t, theta) ode_pendulum(t, theta,b,g,l,m), t_span, thata_ini)

by usind above equation we got output theta is a array which contain values of angular displacement as well as angular velocity

theta(:,1) represent angular displacement.

theta(:,2) represent angular velocity

in ordr to plot them with time i use

figure(1)
plot(t,angular_displacement);
hold on
plot(t,angular_velocity);
hold off
axis([0 10 -4 4])
legend('Angular Displacement','Angular Velocity')

for ploting the beam on which pendulum rest

plot([-1 1],[ 0 0],'linewidth',7)

for ploting the pendulum itself we need two point initial and final points

initial points are x0 = 0, y0= o;

final point are x1= l*sin(theta_diaplacement), y1 = l*cos(theta_dispalcement)

to draw hinges i use marker command.

plot(x0,y0,'O','markersize',3,'markerfacecolor','y')
plot(x1,y1,'O','markersize',20,'markerfacecolor','r')

after that i set the lim of x and y

after completing this we need to capture the frames for that i introduce a variables p= 1

M(p) =getframe(gcf)

and increment of frame is complusory to stiches these images in one video we need to give increment to the frame p=p+1

to create movie i use

movie(M)
videofile = VideoWriter('pendulum.avi;','Uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)

Error:- the only error that m getting is less oscillation of pendulum bc i use sind and cosd intead of sin and cos