Week 3 - Solving second order ODEs

clear all 
close all
clc

b = 0.05; %DAMPING COEFFICIENT
g = 9.81; %gravitational force
l =1;   %length of pendulum
m = 1;  %mass 

%initial conditions as theta 0 and 5

theta_0 =[0;3];

%time points here 100 values are taken betweem 0 to 20
t_span = linspace(0,20,500);

%solve ode and ode function 
%ode45 id inbuilt solver working on Runge-Kutta method, here t is indepent
%where as theta is dependent varaible whos values are to be found.
%[ t , y ] = ode45( odefun , tspan , y0 ) , where tspan = [t0 tf] , integrates the system of differential 
%equations y ' = f ( t , y ) from t0 to tf with initial conditions y0 
%in this case odefun is for (t,theta,b,g,l,m) and t_span with its itnial
%value 0 and final value 5 
%@ is used to define that my ode is fuction t and theta
[t,results] = ode45(@(t,theta) ode_func(t,theta,b,g,l,m),t_span,theta_0); %calling the the ode_fuction 
figure(1)
%plotting the results
subplot(2,1,1)
plot(t,results(:,1))
xlabel('time')
ylabel('angular displacement')
hold on
subplot(2,1,2)
plot (t,results(:,2))

ylabel('angular velocity') %labeling the axes
xlabel('time')
grid on 

 title('VELOCITY & DISPLACEMENT CURVE') %providing title to plot
 

%simulation of pendullum
figure(2)
%value of i is set from 1 : length of results so that loop is executed that
%many times
ct = 1;
for i = 1:length(results(:,1))
    x0 = 0; 
    y0 = 0;
    x1 = 1*sin(results(i,1));
    y1 = -1*cos(results(i,1));
    
     %providing support line for pendulum 
    plot([-1,1],[0,0],'linewidth',3,'color','r');
    
 title('SIMULATION OF PENDULUM') %providing title to plot
  %providing rod of  pendulum 
 line([x0 x1],[y0 y1],'linewidth',2,'color','b'); 
   
   hold on
   %plot for  creating pendulum bob  
   plot(x1,y1,'o','markers',20,'markerfacecolor','g') ;
%contraints of axes
   axis ([-1.5 1.5 -1.5 1.5 ]); 
   pause(0.01); %pause is provide to smoothen the animation
   hold off
   M(ct) = getframe(gcf); %used to capture the frames at different angles
   text  = sprintf("pendulum%d.png",ct);
   saveas(1,text); %saving
   ct = ct+1;
  

end

 %command to produce animation
  movie(M)
 %saving the animation 
  videofile = VideoWriter('pendulum_animation.avi','uncompressed AVI');
  open(videofile);
 %playing the saved animation
  writeVideo(videofile,M);
  close(videofile)
    

......................................................................................



% coversion of second order ode into 2 odes
% theta = o ..symbolic representation
% (d^2*o/dt^2) + b/m(do/dt) + (g/l)sino = 0 .....eq1
% 
% let o = o1
% 
% do/dt = d01/dt = o2
% 
% hence, do^2/dt^2 =  do1^2/dt^2 = d/dt( do1/dt) = d/dt(o2)
% eq1 becomes
% 
% do2/dt + (b/m)o2 +(g/l)sino1 = 0
% 
% i.e do1/dt = o2 
% 
% do2/dt = -(b/m)o2 - (g/l)sino1





function [dtheta_dt] = ode_func(t,theta,b,g,l,m)

%representing the numerics in matlab format
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

https://youtu.be/amRics4KR3o