Solving second order ODE using MATLAB

AIM : To write a program in MATLAB in order to solve the second order ODE of a motion of a simple pendulum with damping.

OBJECTIVE:

1. To write a program in MATLAB that simulates the transient behaviour of the simple pendulum
2. To create an animation of the motion of a simple pendulum using MATLAB.

THEORY:

1. In mathematics, an ordinary differential equation (ODE) is a differential equation containing one or more functions of one independent variable and the derivatives of those functions. The term ordinary is used in contrast with the term partial differential equation which may be with respect to more than one independent variable.
2. The differential equations are classified as linear, Non-linear, Homogenous and Non-Homogenous equations. Based on the order of equation, they are also classified as first order differential equation, second order differential equation or n^th order differential equation.
3. In engineering, ODE is used to describe the transient behaviour of a system and one of it's example is simple pendulum.

Simple pendulum –

A simple pendulum consists of a mass m hanging from a string of length L and fixed at a pivot point P. When displaced to an initial angle and released, the pendulum will swing back and forth with periodic motion, as shown in the below figure

The dynamic behaviour of pendulum bob is important to understand the motion of a simple pendulum. When the bob is given an initial push, when it is at rest, it experiences a force sideways. It is displaced from it's equilibrium position and starts to oscillate about the frictionless pivot point. In general, the motion of a simple pendulum is considered to be simple harmonic motion with constant amplitude and is given by the equation

`(d^2θ)/(dt^2)`+ `g/Lsinθ` `= 0`

where, `θ =`angular displacement of the pendulum bob,

`g =`acceleration due to gravity

`L =`length of the pendulum

In practice, the motion of a simple pendulum is subjected to a restoring force along with friction and air drag which reduces the amplitude of swing with time and bring it back towards the equilibrium position. When the motion of an oscillator reduces due to an external force, the oscillator and its motion are said to be damped. These periodic motions of gradually decreasing amplitude are damped simple harmonic motion. In the damped simple harmonic motion, the energy of the oscillator dissipates contionously. But for a small damping, the oscillations remain approximately periodic. The differential equation for the motion of a simple pendulum with damping is as follows:

`(d^2θ)/(dt^2)+b/m (dθ)/(dt)+g/L sinθ=0`

where,

b - damping coefficient

m - mass of the pendulum

The above differential equation is of order 2 and thus called as an second order ordinary differential equation. By using the above equation, we can simulate the transient behaviour of a simple pendulum using MATLAB.

To solve the 2nd order ODE, we need to create series of first order ODE's.

CALCULATION:

let angular displacement `θ=θ_1`

and the angular velocity will be `(dθ)/(dt)=(dθ_1)/(dt)=θ_2`

`(d^2θ)/(dt^2)=(d^2θ_1)/(dt^2)=(d)/(dt)((dθ_1)/(dt))=(dθ_2)/(dt)`

By substituting the above values in the second order ODE for a simple pendulum, we get

`(d^2θ)/(dt^2)+b/m (dθ)/(dt)+g/L sinθ=0`

`(d^2θ_1)/(dt^2)+b/m (dθ_1)/(dt)+g/L sinθ_1=0`

`(dθ_2)/(dt)+b/m θ_2+g/L sinθ_1=0`

The above equation is a first order differenetial equation, hence we can write as

``therefore```(dθ_1)/(dt)=θ_2`

`therefore``(dθ_2)/(dt)=-b/m θ_2-g/L sinθ_1`

By writing the above equations in the form of a matrix, we get

`(d)/(dt)[[θ_1],[θ_2]]=[[θ_2],[-b/mθ_2-g/Lsinθ_1]]`

where,

`θ_1=`Angular displacement of the pendulum bob

`θ_2=`Angular velocity of the pendulum bob

By solving the above equations in MATLAB, we can get the values of angular displacement and angular velocity with repect to time

MATLAB CODE:

Given that,

Length of the simple pendulum (L) = 1 metre,

Mass of the pendulum (m) = 1 kg,

Damping coefficient (b) = 0.05

Acceleration due to gravity (g) =9 .81 m/s2.

We need to Simulate the motion between 0-20 sec, for angular displacement `(θ_1)`= 0, angular velocity`(θ_2)`= 3 rad/sec at time t=0.

% CODE 1: Function to compute 2nd 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

% CODE 2: Second order ODE of a simple pendulum
clear all
close all
clc

% inputs
b = 0.05; % damping coefficient
g = 9.81; % acceleration due to gravity in m/s^2
l = 1;    % length of the pendulum in 'm'
m = 1;  % mass of the pendulum bob in 'kg'

% initial condition
theta_0 = [0;3]; % angular displacement is '0' and angular velocity is '3'

% time points
t_span = linspace(0,20,600); % pendulum oscillates between 0-20 seconds

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

% plotting angular displacement and angular velocity w.r.t time
subplot(3,1,1)
plot(t,results(:,1),'linewidth',2,'color','b')
axis([0 20 -4 4]);
xlabel('Time (sec)');
ylabel('Angular displacement');
legend('Angular displacement');
grid on

subplot(3,1,2)
plot(t,results(:,2),'linewidth',2,'color','r')
axis([0 20 -4 4]);
xlabel('Time (sec)');
ylabel('Angular velocity');
legend('Angular Velocity');
grid on

subplot(3,1,3)
plot(t,results(:,1),'linewidth',2,'color','b')
axis([0 20 -4 4]);
hold on
plot(t,results(:,2),'linewidth',2,'color','r')
xlabel('Time (sec)');
ylabel('Angular displacement and Angular velocity');
grid on
title('damping effect in simple pendulum');

% CODE 3: Simulation of the motion of pendulum
value = results(:,1);  % The value is the results obtained from the first column

ct = 1;
for i =1:length(t_span);
    VALUE = results(i,1);
    
    % initial conditions
    x0 = 0;
    y0 = 0;
    
    % final conditions
    x1 = l*sin(VALUE);
    y1 = -l*cos(VALUE);
    
    plot([x0 x1],[y0 y1],[-0.5 0.5],[0 0],'linewidth',3)
    hold on
    plot(x1,y1,'marker','o','markersize',30,'markerfacecolor','g')
    axis([-1.5 1.5 -1.5 0.5])
    grid on
    title('Motion of a simple pendulum with damping')
    pause(0.02)
    M(ct) = getframe(gcf)
    ct =ct+1;
    hold off
end
% Animation of the pendulum
movie(M)
videofile = VideoWriter('Simple_Pendulum','uncompressed AVI')
open(videofile)
writeVideo(videofile,M)
close(videofile)

EXPLANATION:

1. Inorder to establish a relation between all the given input values in the second order ODE, a function code is written which is shown in the above code as code 1
2. This function code is based on the matrix which we have derived, Here the first value of LHS row matrix is defined as 'theta1' and second value is defined as 'theta2'.
3. The function which is named as 'ode_func' is then called in the main code with the syntax 'ode45', which is an inbuilt matlab function of the ODE solvers.
4. The code 2 starts with clear all, close all and clc where, 'clear all' command is used to clear the data stored in workspace, 'close all' command is used to close all the previou plots and figures, 'clc' command is used to clear the command window.
5. The given input values of b, g, l and m are entered and the initial condition for angular displacement `θ_1=0`, and angular velocity,`θ_2=3`rad/sec are entered to the variable 'theta_0'.
6. The time span of the motion of simple pendulum 't_span' is taken to be 20 seconds with 600 values using the linspace command, which creates a linearly spaced vector.
7. [t,results] contains the results which we need to plot with respect to time in which results(1) contains angular displacement values in radians and results(2) contains angular velocity values in rad/sec. By using this, we can calculate different values of angular displacement and angular velocity with respect to time. Thus, 'results' is an array.
8. Using the above results, a graph has been plotted between angular displacement vs time, angular velocity vs time, and subplot command has been used inorder to compare all the three graphs at same place.
9. In code 3, counter 'ct' is intialized. Here, 'ct' is a variable that simply increments by a specific step value in every iteration of arovide loop.
10. To provide continous motion to the pendulum, 'for' loop has is used where all the values of t_spn are assigned to a variable 'i' one at a time. The coordinates of the points of line (x0,y0,x1,y1) and their corresponding relations are entered to form a link or the string of a pendulum. This line and it's motion is thus plotted using suitable axis.
11. 'Marker' command is used to form a circle or the bob of a pendulum and the counter is incremented to 1.
12. 'pause' command has been used to achieve the animation speed so that we can observe better.
13. 'getframe' command is used to capture the motion of lines and stored in a variable M.
14. 'Movie' command is used to make a movie with all the frames which are stored in variable M.
15. 'VideoFile' command is used to save the movie and also has an option to store in an uncompreeses 'avi' format along with the file name.
16. Then the video file is opened after the plot, the movie is run and after completion, the video closes and been saved.

OUTPUT:

Animation:

https://youtu.be/sK6wJFEe3uQ

ERRORS:

The errors which i have faced while writing the code for second order ODE of a simple pendulum.

The above error is due to not using the hold off command, due to which all the motions which are captured are overlapped.

RESULTS:

1. From the graph, we can conclude that the periodic curve of angualr displacement and angular velocity reduces with respect to the time, which is mainly due to the damping factor, that is why the motion of a simple pendulum is called as damped simple harmonic motion.
2. From the animation of the simple pendulum, it can be concluded that, as the damping coefficient increases, motion of the pendulum reduces correspondingly and reaches the equilibrium position.