Solving ODE to simulate the motion of a Simple Pendulum

Objective: -

The objective is to write a program in Matlab to simulate the motion of a simple pendulum.

Solution: -

For this current challenge, we are going to solve the following ODE to simulate the motion of a pendulum.

`(d^2theta)/dt^2 + (b/m)((d theta)/dt) + (g/l)sintheta = 0`

Where, b = damping coefficient

m = mass of the pendulum

g = acceleration due to gravity

l = length of the pendulum

The given inputs for the simulation are:

l = 1 metre

m = 1 kg

b = 0.05

g = 9.81 m/s2

Angular Displacement and Angular Velocity are 0 and 3 rad/sec respectively at time t = 0. The given time span for the simulation is for about 20 seconds which with 500 intervals.

Solving the ODE:

Assume `theta = theta_1`,

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

`(d^2theta)/dt^2 = (d^2theta_1)/dt^2 = (d theta_2)/dt`

Comparing the above parameters with the given ODE, we have

`(d theta_1)/dt = theta_2` which describes the Angular Dispalcement

`(d theta_2)/dt = -(b/m)theta_2 - (g/l)sintheta_1` which described the Angular Velocity.

The above assumptions were made into a seperate function to calculate the ODE.

MATLAB Code:

%Input parameters

g=9.81; % gravity in m/s2
l=1; % length of the pendulum in m
m=1; % mass of the ball in kg
b=0.05; % damping coefficient

theta_0 = [0 3]; % Angular Dispacement, Angular Velocity

t_span = linspace(0,20,500);

%ODE using anonymous function
[t ,results] = ode45(@(t,theta) ode_func(l,b,g,m,t,theta),t_span,theta_0);

figure(1)                         
plot(t,results(:,1),'color','r')
hold on 
plot(t,results(:,2),'color','b')
grid on;
xlabel('Time(s)');
ylabel('Results');
legend({'Angular Displacement','Angular Velocity'});
axis([0 10 -5 5])

x=0;
y=0;

ct=1;
for i=1:length(results(:,1))
    
    a=1*sin(results(i,1));
    b=-l*cos(results(i,1));
    
    %Plotting the fixture of the pendulum
    figure(2)
    plot([-1 1],[0 0],'linewidth',4,'color','b');
    hold on
    
    %Plotting the string of the pendulum
    line([x a],[y b],'linewidth',2,'color','r');
    hold on;
    
    %Plotting the ball of the pendulum
    plot(a,b,'-o','markers',15,'markerfacecolor','g');
    axis([-1 1 -1.5 0.5]);
    grid on
    
    hold off
    M(ct)=getframe(gcf);
    ct=ct+1;
end

video = VideoWriter('Pendulum.avi', 'Uncompressed AVI')
open(video)
writeVideo(video,M)
close(video)

The function for Second order ODE:

function[dtheta_dt] = ode_func(l,b,g,m,t,theta);

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

Results: -

The above figure shows the plots for Angular Displacement and Velocity when time t = 0.

The above video shows the motion of the simple pendulum.