Week 3 - Solving second order ODEs
Code Explanation:
First define the basic parameters of the pendulum-like the length, mass, damping coeff etc.
then we solve the governing equations using ode45 inbuilt matlab functions and the we plot the results obtained as time vs dispalcement and angular vel.
to create the animation we hinge the pendulum at x=0 and y=0 and the final position of the pendulum is calculated using the length of pendulum and the displacement of the bob at a particular instant of time.
Animation link:https://youtu.be/SzxGbZeA9UI
Errors while coding:
various errors were encountered while writing the for loop and in the calculation of x1,y1 like i had used column 2 also in the calculation of y1 which resulted in correct animation.
hold off command helped to eliminate the trace of path of the pendulum.
clear all
close all
clc
% define basic parameters
b=0.8; % Damping coeff
g=9.81; % gravity
l=1; % length of pendulum in 'm'
m=1; % mass of pendulum in kg
% initial cond
theta_0=[0;3]; % initial angle o deg and ang vel is 3 rad/s
% time points
t_span=linspace(0,20,200) % for 20 secs
%solve ode
[t,results]=ode45(@(t,theta) ode_func(t,theta,b,g,l,m), t_span,theta_0);
plot(t,results(:,1))% plot first column and all rows
hold on
plot(t,results(:,2))% plot second column and all rows
xlabel('Time(s)')
ylabel('plot')
x0=0;
y0=0;
ct=1;% counter to capture each frame for the animation
for i=1:length(t_span)
% define initial coordinates of the pendulum rod about which it will rotate
x0=0;
y0=0;
% define final coordinates of the bob based on the result array
x1=l*sin(results(i,1));
y1=-l*cos(results(i,1));
figure(2)
plot([x0 x1],[y0 y1],'linewidth',3,'color','b');
hold on
plot(x1,y1,'marker','o','MarkerSize',15,'markerfacecolor','r')
axis([-2 2 -2 2]);% define length of axis
hold off
grid on
legend('Pendulum rod',"bob")
M(ct)=getframe(gcf);
ct=ct+1;
end
% steps to generate an animation video file
movie(M)
videofile=VideoWriter('pendulum.avi');
open(videofile)
writeVideo(videofile,M)
close(videofile)
function[dtheta_dt]=ode_func(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