Solving second order ODE representing equation of motion of damped oscillating pendulum using MATLAB

Aim

To solve the differential equation representing the oscillatory motion of a damped pendulum

Governing equation

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

`theta`= angular displacement of the pendulum

b = damping coefficient

m = mass of pendulum ball (kg)

g = acceleration due to gravity (m/`s^2`)

l = length of string

Objective

The objective of this program is to solve the given second order differential equation and use the solution to simulate the motion of the pendulum governed by this differential equation. Apart from this, we would also plot the variation of displacement and velocity with respect to time.

In MATLAB, a `n^(th)` order differential equation is solved by breaking it into n linear differential equation. For example, a `2^(nd)` order ifferential equation is solved by forming 2 linear differential equations and solving them to get the final result.

Considering the given equation

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

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

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

`implies (d theta_2)/dt + b/m*theta_2 + g/l*sintheta_1 = 0`

`or, (d theta_2)/dt=-b/m*theta_2-g/l*sintheta_1`

Thus, above system of equations are used to solve the given differential equations using ode function of MATLAB. This is done with the help of a function defined to solve the linear differential equations formed followed by using ode function to get the results.

Program

%program to solve ODE
clear all
close all
clc

%inputs
b = 0.05; % damping coefficient
g = 9.81; %acceleration due to gravity
l = 1; %length of pendulum
m = 1; %mass

%initial conditions
theta0 = [0;3]; %displacement,velocity

% time span
t_span = linspace(0,20,500);

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

%plotting
figure(1)
plot(t,results(:,1)) %time vs displacement
hold on
plot(t,results(:,2)) %time vs veocity
xlabel('time')
ylabel('plot')
legend('displacement','velocity');

%pendulum
ct = 1
for i=1:length(results(:,1)) 
    THETA = results(i)
    figure(2)
    x0 = 0;
    y0 = 0;
    x1 = l*sin(THETA);
    y1 = -l*cos(THETA); %negative due to downward direction
 
    p = plot([x0,x1],[y0,y1],'linewidth',2,'color','b')
    axis([-1,1,-1.2,0.1]);
    start_point = [0,0]; %origin of pendulum
    origin_point = viscircles(start_point,0.02); %marking pendulum orgin
    ball = viscircles([x1,y1],0.05); %pendulum mass
    title('motion pf pendulum');
    grid on
    pause(0.001)
    M(ct) = getframe(gcf);
    ct = ct+1;
end

%Creating movie
movie(M)
V = VideoWriter('pendulum.avi','Uncompressed AVI')
open(V)
writeVideo(V,M)
close(V)

First, the required inputs are provided and the initial condition (displacement, velocity) is defined. The time interval (`t_(span)`) is provided as dicrete values in the desired interval.

ode 45 function

It is an inbuilt solver which solves ordinary differential equations. However, it can only solve `1^(st)` order differential equations due to which higher order equations are splitted into number of single order equations. Inside it, a function is defined which describes the parameters affecting the differential equations and a time interval is provided over which the system of equations are integrated. The initial conditions are also provided.

The variable t provides the time. Results is a 3D array containing the values of displacement and velocity at different point of time in column 1 and 2 respectively.

Given below is the function created to pass to ode45

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

For the plots, displacement and velocity obtained in the `1^(st) and 2^(nd)` column of results matrix are each plotted against time.

Output

Simulation of damped oscillating pendulum

Given below is the google drive link to the animation created showing the oscillations of a damped pendulum.

Conclusion

The amplitude of oscillation of a damped pendulum keeps on decreasing with time and finally the pendulum comes to rest. The time required to come to rest is less for higher damping coefficient, that is, higher degree of dampness.

The velocity of pendulum is maximum at the mean position and 0 at the amplitude.