Solving second order ODEs using Matlab - Simple Pendulum with Damping

The problem is to solve the given ODE which represents the 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.

Given Values:

L = 1 metre,

m = 1 kg,

b = 0.05,

g = 9.81 m/s^2.

Time interval = 0-20 sec

Initial Conditions:

angular displacement, `theta` = 0,

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

Solution:

i

As shown in the figure, the position of the bob (x,y) can be calculated with the formula given below:

x = l * sin`theta`

y = l * cos`theta`

where l is the length of the pendulum and `theta` is the angular displacement

Matlab Program:

clear
close all
clc

% Simple Pendulum Simulation
% By Arjun Damodar

% Given
b = 0.05;                    %damping coefficient   
g = 9.81;                    %gravity(m/s^2)
l = 1;                       %length of the pendulum (m)
m = 1;                       %mass of the bob (Kg)

%Initial Condition
theta_0 = [0,3];  %angular dispacement = 0, angular velocity = 3

%Time Interval
t_span = linspace(0,20,100); %t=20sec

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

% Displacement and Velocity Vector with respect to time plotting
plot(t,results(:,1))
hold on
plot(t,results(:,2))
xlabel('Time')
ylabel('Solution Vector')

%Visualization of Simple Pendulum

ct = 1;

for i = 1:20
    
    %Origin
    x0 = 0;
    y0 = 0;
    
	x = l*sin(results(i,1));
	y = -l*cos(results(i,1));
    
    figure(2)
    
	plot([-2,2],[0,0],'g','linewidth',10)
    hold on
    axis([-3, 3 -1.2,0.2])
   	plot([x0,x],[y0,y],'b','linewidth',5)
	plot(x,y,'o','markers',35, 'markerfacecolor','r')
	pause(0.2)
    M(ct) = getframe(gcf);
    hold off
    ct = ct+1;
end

% Animation
movie(M)
videofile = VideoWriter('simple pendulum matlab.avi','Uncompressed AVI');
videofile.FrameRate = 4
open(videofile)
writeVideo(videofile,M)
close(videofile)

    
% Function defining the Solution vectors (Breaking 2nd order ODE to 1st order ODE)

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

1. Engine parameters and initial conditions were assigned.
2. Time step between 0 to 20 sec is defined using 'linspace' command.
3. 'ode45' command is used to differentiate the function which is defined.
4. [t,results] represents the array where the solution of the ODE solver is stored where t denotes the time which means the differentiation of the function is happening with respect to time.
5. ode_func is the function defining the solution vectors where the 2nd order ODE is simplified to a 1st order ODE.
6. The ode solver is assigned with initial conditions, time limit and arguments along with the function which is needed to be differentiated as shown in Step 21.
7. Then the displacement and velocity vector is plotted with respect to time. 'results(:,1)' represents the displacement value of the column vector and likewise 'results(:,2)' represents the velocity values for the corresponding time.
8. Visualisation of the simple pendulum is done by creating the figure needed using the plot and marker command. marker command is used to draw the ball.
9. An array M(ct) is created to store the images or frames representing pendulum position which changes with respect to time and theta.
10. Animation of those stored images in the array is made using the 'movie' command.
11. Frame rate is set to 4fps.

Results:

Screenshot of the output:

Errors Faced:

1. Syntax Error: Syntax for marker command provided wrong by giving the value of marker size inside single quotation.

2. Logical Error: Used 'hold on' before the first plot making which ends up in getting a result as shown below: