Week 4.1 - Solving second order ODEs

% Solving 2nd order ode for pendulum 
clear all
close all
clc

% Input data
L = 1    % meter   
m = 1   % kg
b = 0.05
g = 9.81 % m/s^2

% time span for pendulum
t_span = linspace(0,20,500); % total time taken is 20s where 500 inbetween values are considered

% Intial condition
theta_0 = [0;3]  % 1st coloumn represents displacement while 2nd represents velocity

% Solving ode by ode45
[t,results] = ode45(@(t,theta) ode_func(t,theta,b,m,g,L),t_span,theta_0); 

% Plotting for analysing pandulum transient behavior and pendulum animation

figure(1) % Creating figure window for plotting
plot(t,results(:,1)) % Plotting a graph for time vs results(:,1) which is displacement
hold on              % Hold the figure window 1 for furture plotting
plot(t,results(:,2)) % Plotting a graph for time vs results(:,2) which is velocity
xlabel('Time')                  % Named x-axis as Time
ylabel('Displacement/Velocity') % Named x-axis as Displacement/Velocity
grid on  % make grid visible
hold off % Hold off the figure window 2

ct=1; % Assign value 1 to counter
for i=1:length(results(:,1)) % Assign velue to i 
  
  % Origin or starting position of pendulum which is fixed
  x0=0;
  y0=0;
  
  % Varied value for end point of pendulum
  x1=L*sin(results(i,1));
  y1=-L*cos(results(i,1));
  
  figure(2) % Creating figure window for plotting an animation of pendulum 
  hold on   % hold the figure window 2
  
  plot([-1 1],[0 0],'linewidth',8);     % Plotting fixed roof for pendulum\'s fixed point
  plot([x0 x1],[y0 y1],'linewidth',3); % Plotting link or thred by which pendulum\'s ball is held
  plot(x1,y1,'marker','o','markersize',30,'markerfacecolor','b');     % Plotting the ball of pendulum
  plot(x0,y0,'marker','o','markersize',9,'markerfacecolor',[1 1 1]);  % Plotting the fixed point of pendulum
  line([x0 x0],[y0 -1.5],'linestyle','--','color','r','linewidth',3); % Plotting of centre line 
  
  axis('tight');  % To adjust axis as our requriments
  pause(0.00003); % update the figure at interval of 0.00003s
  
  grid on % Make the grid visible
  
  text = sprintf('test%d.png',ct); % Creat .png file for each image created by plotting
  saveas(gca,text);                % save all .png file 
  
  clf();   % Clearing the current graphic image 
  
  ct=ct+1; % Increment in figure counter
 

end % end for loop


OUTPUT:-


L =

     1


m =

     1


b =

    0.0500


g =

    9.8100


theta_0 =

     0
     3

>> 

function [dtheta_dt] = ode_func(t,theta,b,m,g,L) % creating function named ode_func
  theta1=theta(1); % Assign value theta(1) to theta1
  theta2=theta(2); % Assign value theta(2) to theta2
  
  % find derivatio of theta1 and theta2
  dtheta1_dt = theta(2); 
  dtheta2_dt = -(b/m)*theta2-(g/L)*sin(theta1);
  
  % creating array of dtheta_dt for function which is coloumn matrix
  [dtheta_dt] = [dtheta1_dt; dtheta2_dt];
  
end % end function