Week 3 - Solving second order ODEs
clear all
close all
clc
% Inputs
b = 0.05; % damping coefficient
g = 9.81; % gravitatinal acceleration in meters per sec^2
l = 1; % Length of the pendulum in meters
m = 1; % Mass of the bob in kg
% initial conditions
theta_0 = [0;3];
% time span
t_span = linspace(0,20,300);
% ODE
[t,results] = ode45(@(t,theta)ode_fun(t,theta,b,g,l,m),t_span,theta_0);
plot(t,results(:,1))
hold on
plot(t,results(:,2))
ylabel('Angular Velocity (rad/sec)')
xlabel('time in sec')
legend('Displacement','Angular Velocity')
title('Motion of Pendulum')
theta = results(:,1);
c = 1;
for i = 1:length(theta)
THETA = theta(i);
%intial points
x0 = 0;
y0 =0;
% Position of pendulum
x1 = sin(THETA);
y1 = -1*cos(THETA);
figure(2);
% string Plot
plot([x0 x1],[y0 y1],'color','b');
hold on
% pendulum bob plot
plot(x1,y1,'o','markers',15,'MarkerFacecolor','r')
hold off
axis([-3 3 -2 2]);
pause(0.03)
M(c) = getframe(gcf);
c = c+1;
end
% movie
movie(M);
videofile = VideoWriter('pendulum.avi','Uncompressed avi');
open(videofile)
writeVideo(videofile,M)
close(videofile)
ode function code
function [dtheta_dt] = ode_fun(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
