Solving Second Order Differential Equation in MATLAB for Simulating the motion of the pendulum

Objective: Write program that simulates the motion of the pendulum by solving the equation given below

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.

Input parameters :

L=1 metre,

m=1 kg,

b=0.05.

g=9.81 m/s2.

We are using ODE 45 solver which is used to solve the differential equation of order one. So, reducing the order of the given equation.

`theta = theta_1`

so, `(d theta)/ dt = (d theta_1)/dt = theta_2`

`(d theta_2)/dt + (b/m)*theta_2 +g/l sin theta_1=0`

In this way we have reduced the order of ou problem now we can solve this with ODE45 solver.

Ode45 works Runga Kutta (R-K) 4th order numerical method.

[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 ycorresponds to a value returned in a column vector t.

ode45 accepts that function in which there are only dependent and independent variables as input arguments but in our case their more input arguments that need to given to the solver. So, for that, we have one function name 'ode_func' and values giving thee input arguments in this function.

and while assigning the input argument to the ODE45 solver, we are using anonymous function [@(t, theta)]and feed all the input parameters. By looking to an anonymous function solver will accept the argument and after entering there is no way back then solver to solve.

Main Code

The code below will plot the graph showing relationship between angular velocity and angular displacement with time

clear all
close all 
clc
%Input of the program 
b=0.05;                                     %damping_coefficient
g=9.81;                                     %Accegralation due to gravity
l=1;                                        %Length of the string
m=0.1;                                      %Mass of the ball attached

%initial condition
theta_omega = [0;3]                             %Angular velocity

%time points
t_span = linspace(0,20,500)                 %Time_span_divided into 500 intervals

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

%plotting
figure(1)
plot(t,results(:,1))                        %Plot of angular displacement vs time
hold on
plot(t,results(:,2))                        %Plot of angular velocity vs time
legend('Angular Displacement in radians','Angular velocity in rad/sec');  %legends
ylabel('Plot')
xlabel('time')

Code of ode_Func:

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 of the above code :

Observation: This plot is clearly showing the response as damped vibration system.

Code for animating the motion of the pendulum:

%Marking the poitions of the bob in 2D space
ct = 1;                                     %Counter
for i = 1:length(results(:,1))
x0 = 0;                                     %Fixed poitions(End attached to the support)
y0 = 0;
x1 = l*sin(results(i,1));                   %Point traced by the bob as it moves in 2D Space
y1 = -l*cos(results(i,1));
figure(2)
%Plot joining the points traced by the bob acting as link cooneccting bob and fixed support
plot([x0 x1],[y0 y1],'-o','MarkerSize',10,'MarkerFaceColor','r','MarkerEdgeColor','r','linewidth',2)
axis([-1.3 1.3 -1.3 0.5])
hold on
plot([-1 1],[0 0],'linewidth',4,'color','g')        %Horizontal line to show support
text(-0.5, 0.4, 'Simulation of Pendulum','color','r')
text(0.15, 0.1, 'Fixed Support')
text(x1+0.1,y1+0.1,'Ball')
R(ct) = getframe(gcf);
hold off
ct = ct + 1;
end

%Visualizing the motion of pendulum
movie(R)
RR_Video = VideoWriter('Pendulum_Challenge.avi','Uncompressed AVI')
open(RR_Video)
writeVideo(RR_Video,R)
close(RR_Video)

%Completed

About Code

To calculate the position of the ball at the very next instant. we have used for loop taking values of theta at each instant and simple geometry is required then to plot the next point. In that loop x1, y1 changes their values according to the values of theta.

Output :

Thank You.