Solving second order ODEs using MATLAB

AIM : Write a program in Matlab that will simulate the pendulum motion that solves the given ODE which represents the equation of motion of a simple pendulum with damping.

INTRODUCTION : an ordinary differential equation (ODE) is a differential equation containing one or more functions of one independent variable and the derivative of those functions.

In engineering ODE is used to describe the transient behavior of a system. A simple example is a pendulum. The way the pendulum moves depends on the Newtons second law. When this law is written down, we get a second order Ordinary Differential Equation that describes the position of the \"ball\" w.r.t time.

Similarly, there are hundreds of ODEs that describe key phenomena and if you have the ability to solve these equations then you will be able to simulate any of these systems.

The main objective here is to write a program that solves the following ODE. This ODE represents the equation of motion of a simple pendulum with damping.

In the above equation,

g = gravity in m/s2,

L = length of the pendulum in m,

m = mass of the ball in kg,

b=damping coefficient.

Θ = displacement

The idea behind solving ODE both analytically and conceptually is when given an higher order ODE it can be broken down into simple set of first order ODE so that the problem becomes simpler to solve and execute.

Upon solving the given ODE, we covert it into two first order ode which we obtain as shown in the below fig

Now with the given data,

i.e L=1 metre,

m=1 kg,

b=0.05.

g=9.81 m/s2.

Simulate the motion between 0-20 sec, for angular displacement=0,angular velocity=3 rad/sec at time t=0.

When we do this we will get the position of the pendulum as a function of time. Using this information we need to animate the motion of the pendulum.

PROGRAM

%program to simulate a simple pendulum by solving second order ODE
clear all
close all
clc

% input parameters
l=1;
m=1;
b=0.05;
g=9.81;

% initial condition for a given angular displacement and angular velocity
theta_0 = [0;3];

% time points for 100 values in between
t_span = linspace(0,20,100);

% performing ODE45 command to solve second order ODE
[t, results] = ode45(@(t,theta) ode_func(t, theta, b,g,l,m),t_span, theta_0);

% plots of time vs angular displacement & time vs angular velocity
figure(1)
plot(t,results(:,1))
hold on
plot(t,results(:,2))
ylabel(\'Plot\')
xlabel(\'time\')


% using for loop command for different values of time of the simple pendulum
figure (2)
ct=1;
for i=1:length(results(:,1))
   THETA=results(i,1);
   
%Coordinates
x0=0;
y0=0;

x1=1*sin(THETA);
y1=-l*cos(THETA);

plot([-1.5 1.5],[0 0])
hold on
plot([x0 x1],[y0 y1],\'color\',\'b\')

% plotting of mass by creating a sphere 
plot(x1,y1,\'-o\',\'markersize\',30,\'color\',\'r\')
hold off
axis([-1.5 1.5 -1.5 5.5])
pause(0.0001)

M(ct) = getframe(gcf);
ct = ct+1
grid on
end 

% movie making
movie(M)
videofile=VideoWriter(\'second_ode_pendulum.avi\',\'uncompressed AVI\');
open(videofile)
writeVideo(videofile,M)
close(videofile)

calling function command :

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

Code explanation :

Step 1: Whenever we start a program it is always good to start with the command stating the title of the program which is to be executed.

we then proceed with clear all, close all and clc are entered, which helps in clearing all the computed data in the editor window, command window and workspace so that all the previous data would be cleared to run the new program.

Step 2 : Now the input parameters of the given simple pendulum are entered i.e,

L=1 metre,

m=1 kg,

b=0.05.

g=9.81 m/s^2.

Step 3 : Now to solve the two first order ode , we need to specify the inintial conditions i.e the angular displacements and angular velocity as shown in the program.

Step 4 : With the time limit of 0 to 20 sec for the motion of simple pendulum , we can assume 100 intervals in between.

By using linspace command we can get for each value in time of 0 to 20 sec with 100 intervals which creates an image that represents the position of the pendulum.

Step 5 : Now the ODE 45 command is used to solve the second order ODE, for this to happen we need to call for the function, which in the program is called by \"ode_func\" which is shown in 2nd program, this will help us to simply the program and ODE is solved in here & when it is called it to the main function assigning its input partameters, it provides us the required output.

In the function code, the two first order ode solutions are written and theta1 and theta2 values are assigned to theta(1) & theta(2) respectively to call them in the main function using by the word \"theta\".

Step 6 : Next the plot command is used in the main program to obtain a plot of time vs angular displacement & time vs angular velocity as shown in the program.

Step 7 : Now using the for loop command for obtaining various angles of displacements of the simple pendulum for different values of time in the ODE is created by giving the initial condition x0 and y0, and then x1 y1 conditions are specified for the pendulum.

Step 8 : Next with the help of marker command we create a sphere, by selecting a random color code as shown in the program.

Step 9 : Now the axis is specified and with the help of pause command we can pause the movement.

Step 10 : Then to create a movie and save our plot so that the animation would be like a continuous form, we have a command called ct and let us set the counter to be 1 to save movie, and the frames of various positions of the pendulum will be then stored to ct = ct+1; which will be then called in by getframe command to call all the plotted frames which can be displayed with the given values of pause making it like a movie.

Step 11 : Creating movie – To create a movie, we use a movie command, for all the values obtained for the M variables.

In order to this to happen we need a place/file to store our video, and to do this we use videowriter command to create a new file for storing our movie animations.

Now we open the videos for writing, using open command and then the values obtained from variable M will be written by the command ‘writevideo’

At the end, after the file is written we need to close the file by using the close command.

OUTPUT :

YOUTUBE LINK :

link : https://youtu.be/pTv5M_UOQA8

Error faced :

Error faced during the program was not considering each images getting during the movement of pendulum with different time intervals and this has been notified and corrected in the program and executed in the above program.

Result : Towards the end we can know how to write a program in Matlab that will simulate the pendulum motion that solves the given ODE which represents the equation of motion of a simple pendulum with damping. Similarly, there are hundreds of ODEs that describe key phenomena and if we have the ability to solve these equations then we will be able to simulate any of these systems and solve the real time problems by applying this method.