Week 3 - Solving second order ODEs

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:-

An ordinary differential equation (ODE) is an equation with ordinary derivatives (and NOT the partial derivatives). A differential equation is an equation having variables and a derivative of the dependent variable with reference to the independent variable. A differential equation contains at least one derivative of an unknown function, either an ordinary derivative or a partial derivative. In particular, an ordinary differential equation has ordinary derivatives. Here the ordinary differential equations would be commonly referred to as only differential equations.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 a mechanical arrangement that demonstrates periodic motion. The simple pendulum comprises a small bob of mass ‘m’ suspended by a thin string secured to a platform at its upper end of length L.

The simple pendulum is a mechanical system that sways or moves in an oscillatory motion. This motion occurs in a vertical plane and is mainly driven by the gravitational force. Interestingly, the bob that is suspended at the end of a thread very light somewhat we can say it is even massless. The period of a simple pendulum can be made extended by increasing the length string while taking the measurements from the point of suspension to the middle of the bob. However, it should be noted that if the mass of the bob is changed it will the period remains unchanged. Period is influenced mainly by the position of the pendulum in relation to Earth as the strength of gravitational field is not uniform everywhere.

Meanwhile, pendulums are a common system whose usage is seen in various instances. Some are used in clocks to keep track of the time while some are just used for fun in case of a child’s swing. In some cases, it is used in an unconventional manner such as a sinker on a fishing line. In any case, we will explore and learn more about the simple pendulum on this page. We will discover the conditions under which it performs simple harmonic motion as well as derive an interesting expression for its period.

In general the motion of a simple pendulum is considered to be simple harmonic motion with constant amplitude and is given by the equation,

`frac{d^2theta}{dt^2}+frac{g}{l}sintheta=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,

`frac{d^2theta}{dt^2}+frac{b}{m}(frac{d theta}{dt})+frac{g}{l}sintheta=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=theta1` and angular velocity will be,

`frac{d theta}{dt}=frac{d theta1}{dt}=theta2`

`frac{d^2theta}{dt^2}=frac{d^2theta1}{dt^2}=frac{d}{dt}(frac{d theta1}{dt})=frac{d theta2}{dt}`

Substituting these in the Second Order ODE for a simple pendulum, we get,

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

`frac{d^2theta}{dt^2}+frac{b}{m}(frac{d theta1}{dt})+frac{g}{l}sintheta1=0`

`frac{d theta2}{dt}+frac{b}{m}theta2+frac{g}{l}sintheta1=0`

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

`:.frac{d theta1}{dt}=theta2`

`:.frac{d theta2}{dt}=-frac{b}{m}theta2-frac{g}{l}sintheta1`

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

`frac{d}{dt}[(theta1),(theta2)]=[(theta2),(-frac{b}{m}theta2-frac{g}{l}sintheta1)]`

where,

`theta1`=Angular Displacement of the pendulum bob;

`theta2`=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.

MAIN CODE:-

clear all
close all
clc

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

% initial conditon
theta_0 = [0;3];

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

%solving ODE 
[t, results] = ode45(@(t,theta) Testfunction(t,theta,b,g,l,m),t_span, theta_0);
 
%Plotting Angular displacement and Angular velocity w.r.t Time
figure(1)
plot(t,results(:,1))
grid on
hold on
plot(t,results(:,2))
title('Damping Effect in Simple Pendulum')
xlabel('Time in sec')
ylabel('Angular displacement and Angular velocity')
legend('Angular displacement','Angular velocity')

% create animation
ct = 1;
theta_1 = results(:,1);

for i = 1:length(results(:,1))
    %initial coordinates

   x_start = 0;
   y_start = 0;
   theta_2 = theta_1(i);
   % final cordinates
   x_end = 1*sin(theta_2);
   y_end = -1*cos(theta_2);

   % animation plotiing
   figure(2)
   % pendulum base
   plot([-1 1],[0,0],'LineWidth',4,'Color','r')
   hold on

   % string of pendulum
   line([x_start x_end],[y_start y_end],'linewidth',6,'color','b')
   hold on

   % bob of pendulum
   plot(x_end, y_end,'Marker','o','MarkerSize',28,'MarkerFaceColor','b')
   axis([-2 2 -2 2])
   title('Motion of a Simple Pendulum with Damping')
   pause(0.03)
   hold off
   M(ct)=getframe(gcf);
   ct = ct+1;
end
% Animation of motion of Simple Pendulum
movie(M)
videofile = VideoWriter('simple_pendulum.avi','Uncompressed AVI');
open(videofile);
writeVideo(videofile,M);
close(videofile);

EXPLANATION OF CODE:-

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 `theta1`= 0, angular velocity `theta2`= 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 0 to 20 seconds with 450 equal intervals.
4. To demonstrate the relation between all the input values in second order ODE, a different function code is written which is as follows,

FUNCTION CODE:-

function [dtheta_dt] = Testfunction(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 solver in MATLAB used to integrate differential equations. To solve the differential equation we use an anonymous function. Anonymous functions are those functions which are not stored as a function file, but it is a data type which uses the 'function handle' data type.For differential equations with smooth solutions, ode45 is often more accurate than ode23.

5. [t,results] 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. 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.
9. The ‘marker’ command is used to form a circle or the bob of the pendulum and the counter is incremented by 1.
10. The ‘pause’ command is used to achieve better animation speed.
11. Then ‘getframe’ command is used to capture the motion of the lines and stored in a variable M.
12. By using ‘movie’ command, Matlab is going to make a movie with all the frames stored in variable M
13. 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.
14. 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.

ANIMATION:-

https://youtu.be/TD0OrUbzqDY

Project Link:-

https://skill-lync.com/student-projects/solving-second-order-ode-using-matlab-14