Solving 2nd order ODEs using MATLAB

Aim:

To solve the 2nd order ODE using MATLAB. To represent the motion of the pendulum by solving its ODE of damping motion.

Ordinary Differential Equation:

The ordinar differential equations had a predominant role in the engineering world. The ODEs are used in various mathematical concepts including geometry, analytical mechanics and so on. The damping motion of the pendulum is represented by the ODE. The damping of simple pendulum is an example for the harmonic motion. It velocity and displacement plot of the pendulum differs sinusoidally. The damping motion of the pendulum has an enormous amount of applications. The clock, car shock absorber etc., were the application of simple pendulum.

The damping motion of the pendulum is represented by the following ODE.

where,

g = gravity in m/s2,

L = length of the pendulum in m,

m = mass of the ball in kg,

b=damping coefficient.

Input parameters:

L=1 metre,

m=1 kg,

b=0.05,

g=9.81 m/s2,

time span=0-20s,

angular displacement=0,

angular velocity=3 rad/sec at time t=0.

Procedure:

• Initially, the input parameters are assigned in the program
• The intial position from which the displacement and the velocity is assigned.
• The time period of the oscillation is assigned.
• The ODE can be solved in MATLAB using the variable 'ode45'.
• Once the 'ode45' variable is assigned, the two variables on the ODE is to be given. It is followed by calling the function, and then the time span and the initial position are assigned.
• Once the function is called, the solving procedure of ODE is done in the calling funtion program itself.
• Then, the velocity and displacement motion of the pendulum is plotted.

Animation of motion of the pendulum-

• Initiate a counter variable 'c' to collect all the frames for making the video
• Assigned the theta_0 or the initial position be equal to the result values of displacement.
• The for loop is operated from 1 to the length of the result values of displacement.
• The intial position of the pendulum is assigned as 0.
• The end position of the pendulum can be obtained from solving the sin and cos values of theta_1.
• The theta_1 values are changed corresponding to the changes in the for loop.
• The cos theta_1 value is assigned with a negative sign because the motion is from 180 to 360 (towards gravity).
• The baseline is plotted and then the string is plotted.
• The mass suspended on the pendulum is plotted using 'marker' command.
• The all frames are collected in the 'M' variable.
• The for loop is closed.
• The video file name is assigned and the video file is written on correct format.
• The video file was opened, wrote and then closed to get the motion video of the pendulum

clear all
close all
clc

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

%Initial values
theta_0=[0;3];

%Time period
t_span=linspace(0,20,500);

%Solving second order ODE
[t,result]= ode45(@(t,theta)ode_func(t,theta,b,g,l,m),t_span,theta_0);

%Result plot
figure(1)
plot(t,result(:,1));
hold on
plot(t,result(:,2));

%To make animation of the motion of the pendulum
c=1;
theta_0=result(:,1);
for n=1:length(result(:,1))
    
    %Intial position
    x0=0;
    y0=0;
    
    %End position
    theta_1=theta_0(n);
    x1=sin(theta_1);
    y1=-cos(theta_1);
    
    %Ploting the pendulum
    figure(2)

    %base line
    plot([-2 2],[0,0],'linewidth',2,'color','k')
    hold on

    %string of the pendulum
    plot([x0 x1],[y0 y1],'linewidth',3,'color','k')
    hold on

    %mass suspended on the pendulum
    line(x1,y1,'marker','o','markersize',15,'markerfacecolor','r')
    axis([-3 3 -3 3])
    hold off
    pause(0.03)
    M(c)=getframe(gcf)
    c=c+1;
end

%Animation creating process

movie(M)
videofile=VideoWriter('Pendulum.avi','Uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)

Calling function:

• The theta_1 and theta_2 values are assigned as theta(1) and theta(2).
• dtheta1_dt is assigned to be equal to theta_2.
• dtheta2_dt is equal to -(b/m)*theta2-(g/l)*sin(theta1).
• The dtheta1_dt and dtheta2_dt is assigned in a column matrix.
• Thus, the ODE is solved and it was called on the main 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

Output:

The displacement and velocity plot of the pendulum is representated below,

The image that represent the motion of the pendulum is shown below,

The motion of the pendulum is represented as video and its link is given below,

https://youtu.be/TJYBqt-tHWg

Errors:

1)

The " ' " sign is missed in the "marker" variable and later it was corrected.

2)

The "hold off" command is missed and later it was solved.

Conclusion:

The ODE for simple pendulum is solved and its motion is demonstrated using MATLAB.