Challenge

Week 3 - Solving second order ODEs

VVishnu Visakhan16 Nov 2021

% Program to solve 2nd order ODE using "ode45"
% Program to plot Angular velocity vs Time
% Program to plot Angular displacement vs Time
% Program to Animate the Simple Pendulum
clear
close all
clc

% Inputs
  b=0.05;          % Damping coefficient (constant)
  g=9.81;          % Gravity in m/m^2
  L=1;             % Length of the pendulum in m
  m=1;             % Mass of the ball in kg

% Initial condition
  theta_0=[0;3];   % Initial conditions

% Time points
  t_span=linspace(0,20,500);  % Creating a time interval of 20 seconds

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

% Plotting for tracing of angular velocity and angular displacement
  subplot(2,1,1)                  % Sub polt
  plot(t,results(:,1),'LineWidth',1,'Color','r');
  xlabel('Time in sec');          % Labeling x-axis
  ylabel('Angular displacement'); % Labeling x-axis
  grid on;                        % Enabling the grid on
  subplot(2,1,2)                  % Sub polt
  plot(t,results(:,2),'LineWidth',1,'Color','b');
  xlabel('Time in sec');          % Labeling x-axis
  ylabel('Angular velocity');     % Labeling x-axis
  axis([0 20 -3 3])
  grid on;

% Animation of Simple Pendulum
% Inputs assumed
  X0=0;
  Y0=0;
  ct=1;

for i=1:length(t_span)
    X1=L*sin(results(i,1));
    Y1=-L*cos(results(i,1));
    figure(2)
    plot([-1 1],[0 0],'LineWidth',6,'Color','black');    % Support bar
    line([X0 X1],[Y0 Y1],'linewidth',1.5,'color','b');   % Link of pendulum
    hold on;
    plot(X0,Y0,'O','MarkerSize',5,'MarkerFaceColor','y');  % Fixed point
    plot(X1,Y1,'O','MarkerSize',15,'MarkerFaceColor','r'); % Ball of Simple pendulum

    title('Damping of simple pendulum'); % Assigning title to the Plot
    % Description for the plot
    legend('Support rod','Length of string','Fixed point','Bob of pendulum');
    axis([-3 3 -3 3]);
    
    hold off;
    M(ct)=getframe(gcf);    % Array for storing individual frame
   
    ct=ct+1;
    
end
  movie(M);   % Reading the data stored in Array "M"
% Passing the frames to the videofile 
    videofile= VideoWriter('simple_pendulum'); 
% Opening the video file
    open(videofile);
% Writing the frames to the videofile named "simple_pendulum"
    writeVideo(videofile,M);
% Closing the video file
    close(videofile);

% Function Representing the solution of 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