Solving Second Order ODEs and Simulating the Pendulum Motion using MATLAB

SOLVING SECOND ORDER ODEs

AIM

To solve a second order ordinary differential equation and simulate the pendulum motion using MATLAB.

GOVERNING EQUATIONS

Variables to be considered are damping coefficient (or b), angle made by the rod and vertical(or θ),mass of the rod (or m) and length of the string(or l).

`(d^2theta)/dt^2 + b/m(d(theta))/dt - g/Lsintheta = 0`-----(1)

Now, let `theta = theta_1`

So, we can say that `d(theta)/dt = d(theta_1)/dt = theta_2 (say)`----(2)

Therefore, `(d^2theta)/dt^2 = (d^2theta_1)/dt^2 = d(theta_2)/dt`----(3)

Substituting (2) and(3) in (1), we obtain:

`(d^2theta_1)/dt^2 +b/md(theta_1)/dt + g/Lsintheta_1 = 0`

or, `(d(theta_2))/dt + b/mtheta_2 + g/Lsintheta_1=0`

Since the equation is reduced to a lower order, we can also say that,

`(d(theta_2))/dt = -b/mtheta_2-g/Lsintheta_1`

OBJECTIVES OF THE PROJECT

1. To solve the second order ODE with graphs for velocity and displacement using MATLAB.
2. To simulate the motion of the simple pendulum using MATLAB.

FLOW OF THE PROJECT

• Lower the order of the ODE using substitution method. Prepare a function for the same.
• Set the geometrical and physical parameters of the pendulum like length, mass, damping coefficient as well as gravity.
• Set the initial conditions and timespan for which you want to run the simulation.
• Use the inbuilt ODE solver named ODE45 to solve the ODE which is reduced by calling our function that is made in step 1 and obtain values of displacement and velocity and store it in an array.
• Plot the graph of displacement v/s time and velocity v/s time.
• Create simulation of the positions of the pendulum bob and arm using movie command.

MATLAB CODE

Code for the function program

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:

• The steps mentioned in the GOVERNING EQUATIONS section have been written into code, the only addition being that the final derivative of theta1 and theta2 is stored in an array because of the presence of multiple values.

Code for the main program (using the function program)

clear all
close all
clc

%Input parameters damping coefficient, gravitational force, length and mass
%of pendulum
b = 0.05; g = 9.81; l = 1; m = 1;

% Initial condition
theta_0 = [0;3];  %[angular dispalcement;velocity]

% time points
t_span = linspace(0,20,500);

% solve ODE
[t, results] = ode45(@(t,theta) pendulum_damped_func(t,theta,b,g,l,m),t_span, theta_0);

figure(1)
plot(t,results(:,1)) %plot all rows from first column
hold on  %retain all data from first plot
grid on
plot(t,results(:,2)) %plot all rows from second column
axis([0 20 -3 3]) %resize graph to show X and Y axis from/to x,x(0,20) and y,y(-3,3)
ylabel('Plot')
xlabel('time')

x0=1; y0=1; ct=1;  %setting a hinge point
coordinate=results(:,1);  %all rows from first column which are angle values 
for i=1:length(coordinate)
figure(2) %open a new window,new figure
    x1=1-l*sin(coordinate(i)); %coordinates from trigonometry
    y1=1-l*cos(coordinate(i)); %coordinates from trigonometry
    plot([0 2],[1 1],'linewidth',4,'color','k'); hold on
    plot([1,x1],[1,y1],[x1,x1],[y1,y1]  ,'o','markers',20,'markerfacecolor','r');%plotting the bob using marker and facecolor
    hold off
    axis([0 2 -0.5 1.5])
    grid on
    pause(0.02)
    
    M(ct)=getframe(gcf)
    ct = ct+1;
end

movie(M)  %starting to create a video file
videofile = VideoWriter('Animation_of_Pendulum.avi', 'Uncompressed AVI'); %creating videofile with specific name and avi extension
open(videofile)  %open the video file
writeVideo(videofile,M) %writing plot images to video file
close(videofile) %close the video file

Explanation:

• Physical parameters of the pendulum like damping coefficient, gravitational force, mass of rod and length of string are input.
• Initial condition for the angle between the rod and vertical (as a function of angular displacement and velocity) and the time span are set.
• ode solver is is used by applying the ode45 function with main inputs as t and theta.
• pendulum_damped_func is called with inputs as (t,theta,b,g,l,m).
• the results are stored in an array of [t,results].
• Figure (1) is set to plot all the rows from the first column and the hold and grid functions are turned on tor retain all data from the first plot. Also, rows from the second column are also plotted while resizing the axis.
• y-axis is labelled as ‘plot’ and x-axis is labelled as ‘time’.
• To create animation, a for loop is used to obtain x and y-coordinates using trigonometry. x0 and y0 are hinge points since they are at (1,1).
• pause function is used to give a pause for every iteration of looping of values.
• Finally, movie command is use to create a video file and then it is opened, written and closed to obtain the required animation video.

ERRORS FACED:

1. Used a semi-colon after the M command in the for loop and it showed an error saying

Dot indexing is not available for variables of this type.

Corrected it by removing the semi-colon.

2. Typo in line 32 where I forgot to type the coordinate for y1 correctly, due to which the bob was not visible. Corrected that by adding the left out term.

OBSERVATIONS:

We have obtained a plot of Displacement (blue color) and Velocity (red color) against time and the reduction of velocity and angle displacement can also be observed

CONCLUSION:

• The reduction in velocity and Angle displacement is due to damping coefficient introduced in the Second order ODE of a simple pendulum.
• The second order ODE is successfully solved.
• Youtube link for the animation - https://youtu.be/uwTP2jN7rcg

REFERENCES:

- Figure 1 - Khan Academy