Week 3 - Solving second order ODEs
clear all
close all
clc
% 2nd Order Differential (ODE) Order differential eqn
% Independent variable are as followings.
%g = gravity in m/s2,
%L = length of the pendulum in m,
%m = mass of the ball in kg,
%b=damping coefficient.
l=1 %metre
m=1 %kg
b=0.05
g=9.81 %m/s
function [dtheta_dt] = ODE_Function_of_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
% To Solve:Simulate the motion between 0-20 sec, for angular displacement=0,angular velocity=3 rad/sec at time t=0.
% Initial condition
angular_disp_vel = [0;3];
% Time points
t_span = linspace(0,20,500)
% Solve ODE
[t,results] = ode45(@(t,theta) ODE_Function_of_Pendulum(t,theta,b,g,l,m),t_span,angular_disp_vel)
% Plotting
figure(1)
plot(t,results(:,1),"Color",'b') % angular displacement
grid on
hold on
plot(t,results(:,2)) % angular velocity
xlabel('time')
ylabel('plot')
title('Simple pendulum')
figure(2)
ct=1;
for i = 1:length(t_span)
theta_a = results(i,1)
x1 = l*sin(theta_a);
y1 = -l*cos(theta_a);
plot([-0.5,0.5],[0,0],'LineWidth',3,'Color','y')
line([0 x1] ,[0 y1],'LineWidth',0.003,'color','b')
hold on
plot(x1,y1,'marker','o','MarkerSize',10,'MarkerFaceColor','g')
axis([-2 2 -2 1.5])
grid on
pause(0.00005)
hold off
M(ct)=getframe(gcf)
ct = ct +1
end
movie(M)
videofile = VideoWriter('ODE_equation_pendulum.m','Uncompressed AVI')
open(videofile)
writeVideo(videofile,M)
close(videofile)
Explanation:
Given data :
L=1 metre,
m=1 kg,
b=0.05.
g=9.81 m/s2
% Initial condition
angular_disp_vel = [0;3];
It mean angular displ and velocity starts from 0 and 3
% Time points
t_span = linspace(0,20,500)
It means integration of time from 0 to 20 sec
% Solve ODE
[t,results] = ode45(@(t,theta) ODE_Function_of_Pendulum(t,theta,b,g,l,m),t_span,angular_disp_vel)
It means ode45 is a inbuilt function in matlab in which we are using ODE_Function_of_Pendulum as a sub program for soliving 2 nd order diffrential eqn
figure(2)
ct=1;
for i = 1:length(t_span)
theta_a = results(i,1)
It means for the length of t_span = 500
x1 = l*sin(theta_a);
y1 = -l*cos(theta_a);
Coordinates for pendulum at different angles is shown above
plot([-0.5,0.5],[0,0],'LineWidth',3,'Color','y')
line([0 x1] ,[0 y1],'LineWidth',0.003,'color','b')
hold on
plot(x1,y1,'marker','o','MarkerSize',10,'MarkerFaceColor','g')
axis([-2 2 -2 1.5])
grid on
pause(0.00005)
hold off
M(ct)=getframe(gcf)
ct = ct +1
end
movie(M)
videofile = VideoWriter('ODE_equation_pendulum.m','Uncompressed AVI')
open(videofile)
writeVideo(videofile,M)
close(videofile)
Above steps are same for every animation videos .
Errors:
Typing errors and i forget to put hold on in plotting
Conclusion:
We can conclude that as we increase damping coefficient other dependent variable will decrease (i.e angular velocity and angular displacement)