Week 4.1 - Solving second order ODEs
AIM:
To write a program in Matlab\Octave that will simulate the pendulum motion.
GIVEN DATA:
In the above equation,
g = gravity in m/s2,
L = length of the pendulum in m,
m = mass of the ball in kg,
b=damping coefficient.
use,
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.
1)CODE:
clear all
close all
clc
% Variable defenition
b=0.05; % damping coefficient
g=9.81; % gravity in m/s2
l=1; % length of the pendulum in m
m=1; % mass of the ball in kg
% initial condition
theta_0=[0;3];
% Time points
t_span=linspace(0,20,500);
% Solving ODE
[t,result]=ode45(@(t,theta) ode_sk(t,theta,b,g,l,m),t_span,theta_0);
ct=1;
theta=result(:,1);
for i = 1:length(t_span)
Theta1 =theta(i);
x0 = 0;
y0 = 0;
x1 = sin(Theta1);
y1 = -cos(Theta1);
% plotting
plot([x0 x1],[y0 y1]);
hold on
axis([ -2 2 -2 2]);
plot(x1,y1,Marker='o',MarkerFaceColor='r',MarkerSize=15);
hold off
xlabel('Time')
ylabel('Displacement & Velocity')
title('Pendulum')
grid on
M(ct) = getframe(gcf);
ct = ct+1;
end
movie(M);
videofile = videowriter('pendulum.avi','Uncompressed AVI');
open(videofile);
writeVideo(videofile,M);
close(videofile);
Function file:
function [dtheta_dt] = ode_sk(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
2) OUTPUT:
theta1 =
0
theta2 =
3
theta1 =
4.0190e-05
theta2 =
3.0000