Solving Second Order ODE using MATLAB

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

OBJECTIVE :

• To write a program to simulate the transient behaviour of a simple pendulum
• To create an animation of the motion of simple pendulum using MATLAB.

THEORY :

Ordinary Differential Equation (ODE) –

In mathematics, an ordinary differential equation (ODE) is a differential equations containing one or more functions of one independent variable and the derivatives of those functions. An ordinary differential equation (ODE) is an equation that involves some ordinary derivatives (as opposed to partial derivatives) of a function. The term ordinary is used in contrast with the term partial differential equation which may be with respect to more than one independent variable. 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 Ordinary Differential Equation, Second Order Ordinary Differential Equation or ‘nth’ Order Equation. In engineering, ODE is used to describe the transient behaviour of a system and one of its example is a simple pendulum.

Simple Pendulum –

A simple pendulum is defined to have an object that has a small mass, also known as the pendulum bob, which is suspended from a light wire or string, such as shown in the figure below.

The dynamic behaviour of pendulum bob is important to understand the motion of simple pendulum. When the pendulum bob is given an initial push while it is at rest, it experiences a force sideways. It is displaced from its equilibrium position and starts to oscillate about the pivot point swinging back and forth. 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/l sin⁡θ=0`

Where,

`theta` = angular displacement of the pendulum bob;

g = acceleration due to gravity in m/s²;

`l` = length of the pendulum in m.

In practice, the motion of a simple pendulum is subjected to a restoring force along with friction and air drag which reduces their 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 continuously. 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 bob in kg.

This differential equation is of order ‘2’ and thus it is called as Second Order Ordinary Differential Equation (ODE) for a simple pendulum. By using this equation we can simulate the transient behavior of a simple pendulum using Matlab.

Derivation for code simplification –

Let angular displacement `theta` = `theta_1` and 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`

Substituting these 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 θ)/(dt^2 )+b/m ((dθ_1)/dt)+g/l sin⁡θ_1=0`

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

This is a first order differential equation, hence we can write,

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

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

Writing these equations in the form of a matrix, we have,

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

where,

`θ_1`= Angular Displacement of the pendulum bob;

`θ_2` = Angular Velocity of the pendulum bob.

By solving the above equation in Matlab, we can get the values of angular displacement and angular velocity with respect to time.

MATLAB Program :

Inputs –

The values of length of the pendulum '`l`' is taken as 1 m, mass of the pendulum bob ‘m’ is taken as 1 kg and damping coefficient ‘b’ is taken as 0.05. The motion of simple pendulum is simulated between 0-20 sec, for angular displacement '`θ_1`' = 0, angular velocity '`θ_2`' = 3 rad/sec at time t=0.

%Second Order ODE of a Simple Pendulum and its Animation
clear all
close all
clc

%Inputs
b = 0.05;     %damping coefficient
g = 9.81;     %acceleration due to gravity in m/s2
l = 1;        %length of the pendulum in m
m = 1;        %mass of the pendulum bob in kg

%Initial condition
theta_0 = [0;3];

%Time span
t_span = linspace(0,20,600);

%Solving ODE
[t,y] = 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
plot(t,y(:,1),'linewidth',2);
hold on;
plot(t,y(:,2),'linewidth',2);
axis([0 20 -4 4]);
xlabel('Time (sec)');
ylabel('Angular displacement and Angular velocity');
grid on
title('Damping Effect in Simple Pendulum');
legend('Angular displacement','Angular velocity');
pause(10);

ct = 1;
for i = 1:length(t_span);    
    x0 = 0;
    y0 = 0;
      
    x1 = l*sin(y(i,1));
    y1 = -l*cos(y(i,1));
    
    %Simulation of motion of Pendulum
    plot([x0 x1],[y0 y1],[-0.5 0.5],[0 0],'linewidth',3)
    hold on 
    plot(x1,y1,'marker','o','markerSize',25,'markerFacecolor','r')
    axis([-1.5 1.5 -1.5 0.5])
    grid on
    title('Motion of a Simple Pendulum with Damping')
    pause(0.01)
    M(ct) = getframe(gcf);
    ct=ct+1;
    hold off
end

%Animation of motion of Simple Pendulum
movie(M)
videofile = VideoWriter('Simple_Pendulum.avi','Uncompressed AVI')
open(videofile)
writeVideo(videofile,M)
close(videofile)

1. The code starts with clear all, close all and clc where,

• ‘clear all’ command is used to clear the data stored in the workspace.
• ‘close all’ command is used to close all plots and figures.
• ‘clc’ command is used to clear the command window.

2. The given input values of b, g, `l` and m are entered and the initial condition for angular displacement '`θ_1`' = 0, angular velocity '`θ_2`' = 3 rad/sec are entered to the variable ‘theta_0’.
3. The time span of the motion of simple pendulum ‘t_span’ is taken to be 20 seconds with 600 equal intervals.
4. To establish the relation between all the input values in second order ODE, a different function code is written which is as follows,

   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

• This function code is based on the matrix which we derived. Here the first value of LHS row matrix is defined as ‘theta1’ and second value is defined as ‘theta2’.
• Then all the corresponding equations and relations are entered. This function which is named as ‘ode_func’ is then called in the main code with a special syntax ‘ode45’ .
• ode45 – ode45 is an inbuilt Matlab function of choice among the ODE solvers. It compares 4th and 5th order methods to estimate error and determine step size. For differential equations with smooth solutions, ode45 is often more accurate than ode23.

5. [t,y] contains the results which we need to plot with respect to time in which y(1) contains angular displacement values in radian and y(2) contains angular velocity values in radian. By using this we can calculate different values of angular displacement and angular velocity with respect to time. Thus ‘y’ is an array.
6. Using these results, the graph of angular displacement and angular velocity with respect to time is plotted.
7. Then counter ‘ct’ is initialized. Here ‘ct’ is a variable that simply increments by a specific step value in every iteration of a loop.
8. To provide continuous motion to the pendulum, a ‘for’ loop is used where all the values of ‘t_span’ 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 the link or string of the pendulum. This line and its motion is thus plotted using suitable axis. The ‘marker’ command is used to form a circle or the bob of the pendulum and the counter is incremented by 1.
9. The ‘pause’ command is used to achieve better animation speed.
10. Then ‘getframe’ command is used to capture the motion of the lines and stored in a variable M.
11. By using ‘movie’ command, Matlab is going to make a movie with all the frames stored in variable M
12. To save this movie, the ‘videofile’ command is used which has the option to store the video in an Uncompressed ‘avi’ format along with the file name.
13. Then the video file is opened after the plot, the movie is run and after completion the video stops or closes.

Program Output :

RESULTS :

1. As seen in the graph, the periodic curve of Angular Displacement and Angular Velocity reduces with respect to time. This is mainly due to the damping factor. It can also be stated as the simple pendulum is subjected to a restoring force due to damping factor that will accelerate it back towards the equilibrium position. This is why the motion of simple pendulum is called as damped simple harmonic motion.
2. From the animation of motion of the simple pendulum, it can be concluded that, as the damping coefficient is increased, the motion of the pendulum reduces correspondingly and approaches the equilibrium position. Also increase in damping coefficient reduces the time taken by the pendulum to reach its equilibrium position and vice-versa.

ERRORS :

1. Invalid function name “sin.(y(i,1))” –Solution – The Dot ‘.’ Operator is incorrect as there is only one value of ‘y’ and not an array of values. Removing the dot operator gives the correct solution.

2. Error in movie M(ct) command –Solution – The counter ‘ct’ must be incremented after the M(ct) = getframe(gcf) command. If the counter is incremented before, then the movie will contain uninitialized frames.

ANIMATION :