Week 3 - Solving second order ODEs using Matlab

Aim- To write a program that solves the following ODE. This ODE represents the equation of motion of a simple pendulum with damping. And also write a program to simulate the transient behaviour of a simple pendulum and to create an animation ot its motion.

Theoy- A simple pendulum consists of a point mass m hanging on a string of length L from one support. One pulls the mass back to some starting angle θ0θ0 and releases it than the goal is to find the angular position as a function of time. The way the pendulum moves depends on the Newton's Second Law of Motion. When this law is written down, we get a second order differential eqauation that describes the position of the ball w.r.t time. The simple pendulum is described by the below mentioned ODE-



Given that
g = 9.81 (m/s^2)
L = 1 meter
m = 1 kg
b = 0.05
Simulate the motion between 0-20 sec.
For angular Displacement = 0
Angular velocity = 0.3 rad/sec
at time t=0.

Matlab Code
1) Main code

%Matlab Code of Solving second order ODE's 
clear all
close all 
clc
%Inputs
b=0.05;
g=9.81;
l=1;
m=1;
% Initial conditions
theta_0 = [0;3];
%Time points
t_span = linspace(0,20,400);
%solve ode
[t,results] = ode45(@(t,theta) ode_func(t, theta, b,g,l,m),t_span, theta_0);
plot(t,results (:,1))
hold on
plot(t,results(:,2))
ylabel('plot')
xlabel('time')
fig=1;
for i = 1:length(t)
    
%initial position of the pendulum
    
    x0 = 0;
    y0 = 0;
    
    %final position
    
    x1 = l*sin(results(i,1));
    y1 = -1*cos(results(i,1));
    
%plotting  
    figure(2)
    %base line
    plot([-1,1],[0,0],'linewidth',6,'color','k');
    % pendulum line
    line([x0 x1],[y0 y1],'linewidth',4,'color','b');
    hold on;
    %pendulum ball
    plot(x1,y1,'o','markers',25,'markerfacecolor','g');
    grid on
    axis([-2 2 -2 2]);
    pause(0.01);
    hold off;
    M(fig) = getframe(gcf);
    fig= fig+1;
    end
% for genarating and storing the video file in avi format
  movie(M)
  videofile = VideoWriter('second_order_ode_simple_Pendulum.avi','Uncompressed AVI');
  open(videofile)
  writeVideo(videofile,M)
  close(videofile)

2)Function Code

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

Explanation of code -
1) To solve the ODE an ode function is required which accepts the values of t, theta,b,g,m,t_span,theta_0, so a fuction called ode_func is declared and the syntax foe the function is written in the part 1 of coding. The function consists of function declaration an ode_func is the function statement or name and it accepts the values of t,theta,b,g,L,m. In the next two lines the second order ode is split in to two parts and in the end both the terms dtheta1_dt,dtheta2_dt are made in to an array of single variablee dtheta_dt which is a column matrix.
2) Now going to the main code in the starting the input variables like b,g,L,m are declared.
3) As we are solving higher order odes they are split in to first order ode velocity,displacement, in the below line displacement is 0 and angular velocity is 3 and the are stored in a array theta_0.
4) To obtain solution we need to intrigate the ode which requires values for time therefore the time points are generated using linspace command which generates 400 values between 0 to 20.
5) [t,results] = ode45(@(t,theta) ode_func(t, theta, b,g,l,m),t_span, theta_0);
In the above line ode45 is a default solver in matlab and it accepts the values ode t,theta as in inputs ode_func is our function which we have written earlier. t_span is the time span and theta_0 is the initial condition consisting velocity & displacement values as an array.
6)
plot(t,results (:,1))
hold on
plot(t,results(:,2))
ylabel('plot')
xlabel('time')
The above lines are used to plot the result values in a graph and hold on command is used to generate a second plot in the same figure window, the generated plot shows the displacement & time.
7)
fig=1; % It generates the figure window that generates the plot
for i = 1:length(t) % Initiation of for loop with index i

%initial position of the pendulum

x0 = 0;
y0 = 0;

%final position

x1 = l*sin(results(i,1));
y1 = -1*cos(results(i,1));
In the above lines the condition for initial position of the pendulum & the various positions the pendulum takes place w.r.t to time & angle is given

%plotting
figure(2) % It generates a figure window 2

%base line
plot([-1,1],[0,0],'linewidth',6,'color','k');
It explains the base line for the pendulum which acts as a hinge point with linw width 6 and color black

% pendulum line
line([x0 x1],[y0 y1],'linewidth',4,'color','b');
It generates the pendulum arm with blue color in the same figure window

hold on; % It keeps the current plot while adding new element to the plot

%pendulum ball
plot(x1,y1,'o','markers',25,'markerfacecolor','g'); % It generates the pendulum ball in green color of 'o' shape
grid on % It switches on grid

axis([-2 2 -2 2]); % Axis ares for thr plot
pause(0.01); %Pause value for the animation
hold off;
M(fig) = getframe(gcf); % captures the frame generated into animation
fig= fig+1; % loop encounter increment used for multiple generation of plot
end % end of for loop

% for genarating and storing the video file in avi format
movie(M)
videofile = VideoWriter('second_order_ode_simple_Pendulum.avi','Uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)
The above lines are used to generate the animation & to store it avi format.

SCREENSHOTS OF THE RESULTS

RESULTS - Values for displacement in 1st and velocity in 2nd column.

CONCUSION - The governing equation for simple pendulum with damper is solved & the result showing position & angular velocity is plotted. Animation is made using Matlab coding to show the motion of simple pendulum governed by a second order ode.

NOTE - THE GOOGLE DRIVE LINK FOR THE ANIMATION IS GIVEN BELOW -

https://drive.google.com/file/d/1-7Xt_fE6KnCtu-cVxojSTMP18VW6gsET/view?usp=sharing