Week 3 - Solving second order ODEs
% program to solve 2nd order ode for the pendulum
close all
clear all
clc
% inputs from the challenge
L = 1; % length of pendulum
m =1; % mass of ball
b =0.05; % damping coefficient
g = 9.81; % gravitational acceleration
% initial condition
theta_0 = [0;3];
% time points
t_span = linspace(0,20,500);
% ode solution
[t results] = ode45(@(t,theta)ode_func(t;theta,b,g,l,m),t_span,theta_0);
% plot results
figure(1)
plot(t,results(:,1))
hold on
plot(t,results(:,2))
xlabel('time')
ylabel('plot')
% creating animation
ct = 1; % counter varaiable to keep track of frame number
for i=1; length(results(:,1))
x0 = 0;
y0 = 0;
x1 = 1*sin(results(i,1));
y1 =-1*cos(results(i,1));
figure(2)
plot([-1 1],[0 0],'linewidth',10,'color','g') % plot fixed support
axis([-2 2 -2 2])
hold on
line([x0 x1],[y0 y1],'linewidth',5,'color','b') % plot string position
plot(x1,y1,'0','markersize',25,'markerfacecolor','r') % plot ball position
hold off
M(ct) = getframe(gcf); % stores current frame
ct = ct+1;
end
% create movie
movie(m)
videofile =VideoWritier('pendulum_simulation.avi','uncompressed AVI'); % create videofile
open(videofile); % open videofile
WriteVideo(videofile,M) % write frames to videofile
close(videofile) % close videofile
animation link :
ode_func :
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
note : ode_func converts second order DE into two single order DE