Week 3 - Solving second order ODEs

%program for simulating the pendulum motion


clear all
close all
clc

%inputs

b = 0.05;   %damping co-efficient
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,200);

%Solve ODE
[t,results] = ode45(@(t,theta)ode_func(t,theta,b,g,l,m),t_span,theta_0);
hold on
plot(t,results(:,1),'linewidth',3,'color','b')
hold on
plot(t,results(:,2),'linewidth',3,'color','r')
xlabel('time')
ylabel('amplitude')
grid on

legend(['Displacement(rad)','Velocity(rad/sec)'])
ct=1;

for i=1:length(t_span)
    x=l*sin(results(i,1));
    y=-l*cos(results(i,1));
    
    %plotting
    
   figure(2)
   
   
   %hold on
   
   plot(0,0,'o','Markersize',10,'color','b','Markerfacecolor','b') %allocting co-ordinates for fixed support
   plot([0 x],[0 y],'linewidth',2,'color','b') %allocting co-ordinates for pendulum string
   hold on
   plot(x,y,'-o','markers',10,'markerfacecolor','r') %allocting co-ordinates for BOB         
   grid on                
 
  %limits of axis
  axis([-2 2 -2 2]);
  
  %collecting frames of figures to create animation
  M(ct)=getframe(gcf);
  hold off
  ct=ct+1;
end

%for creating an animation of pendulum

movie(M)
Videofile = VideoWriter('pendulum_animation.avi','Uncompressed AVI')
open(Videofile)
writeVideo(Videofile,M)
close(Videofile)





%code for making function to solve ode

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