Week 3 - Solving second order ODEs

Matlab Code for solving second order Equation.

Initially, the function is Created, where the Second-order Equation is transformed into a simple Equation for solving in Matlab

Transformation of 2nd order ODE is clearly understood from the video session given by the instructor and please refer to the notepad attached.

function [dtheta_dt] = ode_func(t,theta,b,g,l,m)
theta1 = theta(1)
theta2 = theta(2)
dtheta1_dt = theta2  %  Euation 1 
dtheta2_dt = -(b/m)*theta2-(g/l)*sin(theta1); % Equation 2.
dtheta_dt = [dtheta1_dt;dtheta2_dt];
end

The Function is called in the main code to solve the Second-order using the inbuilt Matlab ODE45 Function.

clear all
close all
clc

%% Inputs 

b = 0.05; %  Damping Coefficient
g = 9.81; %  Gravitional Force in 'M/S2'
l = 1;    %  Length of the pendulum in 'M'
m = 1;    %  Mass or BOB at the end of Pendulum in 'KG' 

% initial condition

theta_0 = [0;3]; % angular velocity=3 rad/sec

t_span  = linspace(0,20,100); % Simulated the motion between 0-20 sec


%% solve ODE Syntax

% t,X :Solution Vectors ,Matlab returns X for each time t
% ode_func : funtion that returns dtheta2_dt = -(b/m)*theta2-(g/l)*sin(theta1)
% t_span : inital time and final time of the pendulum
% theta_0 : final value untill which the solution is resolved.

[t,X]= ode45(@(t,theta) ode_func(t,theta,b,g,l,m),t_span,theta_0);

T = 2*pi*sqrt(l/g); % time taken for one Oscillation. 

%% plotting the results 

figure(1)
plot(t,X(:,1),LineWidth=3)
hold on
plot(t,X(:,2),LineWidth=2)
xlabel('Time (t) in Seconds')
ylabel('angular velocity in rad/Sec')
title('Second Order ODE')
%% Animation
l1 = 1.5; % First Arm Link
l2 = 1.2; % Second Arm Link

theta1 =90;
theta2 =0;

x = 0.0;
y = 0.0;

ct = 1;

    x = 0; % Initial x coordinate
    y = 0; % Initial y coordinate

    x1 = cosd(theta1)*l1;  % x1 coordinate for first link
    y1 = sind(theta1)*l1;  % y1 coordinate for first link

    x2 = x1+cosd(theta2)*l2; % x2 coordinate for second link
    y2 = y1+sind(theta2)*l2; % y2 coordinate for second link

    % for loop used to oscillate the pendulum with time;t

  for i = 1:length(t)  

  x3 =l*sin(X(i,1))+x2  % x3 coordinate + forward motion of the pendulum in Rad
  y3 =-l*cos(X(i,1))+x2 % y3 Coordinate + reverse motion of the pendulum in Rad
           
  figure(2)
  plot([x x1],[y,y1],[x1 x2],[y1 y2],[x2,x3],[y2,y3],LineWidth=3) % PLotting with Arm link
  axis([0 2.5 0 2]);  % Axis set to have clear picture
  hold on             %To have multiple points in same plot
  plot(x1,y1,'O','markersize',5,MarkerFaceColor='yellow') % Joint 1
  plot(x2,y2,'O','markersize',5,MarkerFaceColor='red') % Joint 2
  plot(x3,y3,'O','markersize',20,MarkerFaceColor='black') % Joint 3 as shpere
  title('Damping of Simple Pendulum')
  legend('1st Arm Support','2nd Arm support','Length of Pendulum','joint 1','joint2','mass',Location='best')
  grid on  % turns on the grid
  hold off % Deletes mulitple point in the same figure
  M(ct)=getframe(gcf); % Collecting the frame in vector"M"
  ct= ct+1; % Increasing the counter to add to create Movie

  end

  %% creating Movie.
movie(M)
videofile = VideoWriter('simple pendulum-2','Uncompressed AVI')
open(videofile)
writeVideo(videofile,M)
close(videofile)