Week 3 - Solving second order ODEs

Aim: To simulate the Motion of a Pendulum using Matlab by solving the second-order Differential Equation.

Theory:

• The motion of a Pendulum under gravity in the atmosphere is defined by using is defined using a Second-order Differential Equation shown below:

`(d^2theta)/dt^2 + b/m * (d theta)/dt +g/L sin(theta)=0`

Where:

b= damping coefficient

m= mass of the pendulum

g= Acceleration due to gravity

L= length of Pendulum

It is further simplified into set of 2-1st order differential equation in matrix form:

• In order to solve a 2^nd order differential equation, the equation needs to be broken down into a set of 1^st order differential equation as the solvers in Matlab requires 1^st order differential equations.
• Out of the different ODE solvers available in Matlab, ode45 is the most common option being suitable for almost every type of problem.
• Syntax of ode45:
  • [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. Each row in the solution array y corresponds to a value returned in column vector t.
    • Reference to: https://www.mathworks.com/help/matlab/ref/ode45.html
  • odefun in the syntax is where the function and its variables are defined using a particular syntax shown in the code.

Main Code:

%---------------Section1-----------------
clear all
close all
clc
 
%--------------section2------------------
%Declaring constants
b=0.02;
g=9.81;
l=1;
m=0.1;
 
%--------------Section3----------------------
theta_0= [0;3];                %initial conditions
 
t_span=linspace(0,20,500);     %limits of integration
 
[t,result]=ode45(@(t,theta) ode_func(t,theta,b,g,l,m), t_span, theta_0);
 
%-------------Section4-----------------------------
figure (1)
plot(t,result(:,1));         %plotting displacement vs time
hold on;
plot(t,result(:,2));         %plotting velocity vs time
 
%------------Section5--------------------------------
[a]=result(:,1);              %displacement magnitude along the arc. [angle]=a./l;                 %angle of the pendulum string wrt vertical.
 
figure (2)
for i=1:length(a)
   
    figure (2)
    %plotting the fixed point and the ball end of the pendulum
    plot([0,-l.*sin(angle(i))],[0,-l.*cos(angle(i))],'-*r');
    hold on
    
    plot([-2,2],[0,0],'-b');               %plotting fixed wall
    axis([-2 2 -2 2]);
   
    hold off
    
    M(i)= getframe;          %capturing plot frames from every iteration
end
 
%------------Section6------------------------------------------
movie(M);
file= VideoWriter('Pendulum.avi','uncompressed AVI');
open(file);
writeVideo(file,M);
close(file)

Defined Function:

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

Main Code Explanation:

Section1:

• The commands clear all the variables from the workspace, text from the command window and also all the figures generated from the previous set of codes.

Section 2:

• The constants required are initiated.

Section3:

• The initial conditions i.e. displacement and Velocity magnitude are initiated in a column matrix theta_0.
• The time span (t_span) of 20 seconds with 500 points in between is initiated which would serve as the limits for the integration of the differential equations.
• Ode45 is used in the above-mentioned Syntax.
  • @(t, theta) represents the declaration of the independent and dependent variable in the differential equation.
  • ode_func is a function defining the terms in the answer matrix shown in equation 1. The function is explained later.
  • The Initial condition column vector (theta_0) and the time span (t_span) row vector is given as an input.
  • The result of the ode function [t, result] gives the Timespan in column matrix format (t), and ‘result’ contains a 2 column matrix with the first column containing the values of displacement and second column containing values of velocity respective to time (t).

Section 4:

• The values of Displacement and Velocity from the result matrix are plotted against the time value from the t matrix.

Section 5:

• The value of theta is calculated using the relation of the length of the arc which is basically the displacement value from the result matrix.
• The For loop has a counter corresponding the no. of time steps. The origin is taken as the pivot of the pendulum and the position of the ball is plotted using the relation shown above and a negative sign is used since the ball lies in the Negative Y quadrant.
• The fixed support is plotted as a line and also the axes for the plot are specified.
• Every plot in a single pass of the For loop is captured and stored in the matrix M using a ‘getframe’ command.

Section 6:

• The ‘Movie’ function plays the movie defined by the input array M which constitutes frames from ‘getframe’ command.
• The next set of code is used to create a video file named ‘Pendulum’ in an uncompressed AVI format, open it, Write the video into the file and close it.

Defined function explanation:

• The function ode_func is defined using the variables dependent and independent and the constants required.
• Inside the function, the relations based on equation 1 are stated and result is a 2 column matrix including displacement and Velocity magnitudes.

Output Graph :

Result Video :

The animation of the pendulum is uploaded on Youtube :

https://youtu.be/6tOMVmDieR8