Week 3 - Solving second order ODEs

Second order ODE soving using matlab

In mathematics, an ODE is containing one or more functions of one independent variable and the derivative of thode functions. An ODE is 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 equati0ns are classified as Linear , Non-linear ,Homogeneous and Nonhomogeneous equations . Based on the order of the equation , they are also classified as First order ODE, second order ODE, or N’th order ODE. In engineering, ODE is uded to describe the transient behaviour of a system and one of its example is simple pendulum.

SIMPLE PENDULUM:-

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

The dynamic behaviour of the pendulum bob is important to understand the motion of simple pendulum. When the pendulum bob is given initial push while it is in rest ,it experiences a force side ways.It is displaced from its equilibrium position and start to oscillate about 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ɵ/dt)^2+(g/l)*sinɵ=0;

Where

ɵ=angular displacement of the pendulum bob

g=acceleration due to gravity in m/s^2

l=length of the pendulum string in m

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

(dɵ/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 ODE for example Pendulum. By using this equation we can simulate the transient behaviour of a simple pendulum using matlab.

Derivation for code simplification-

Let angular displacement ɵ= ɵ1 and angular velocity will be ,

(dɵ/dt)=(dɵ1/dt)=ɵ2;

(dɵ/dt)^2=(dɵ1/dt)^2=(d/dt)(dɵ/dt)=(dɵ2/dt);

Substituting these in the second order ODE for a simple pendulum, we get,

(dɵ/dt)^2+(b/m)(dɵ/dt)+(g/l)sinɵ=0;

(dɵ1/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 ,

(dɵ1/dt)=ɵ2;

(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 w.r.t time.

MATLAB Program:

Inputs-

The values of length of the pendulum string ‘l’is taken as 1m , mass of the pendulum bob ‘m’ is taken as 1kg 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 timet=0.

%second order ODE of a simple pendulum and its animation
clear  all
close all
clc

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

%initial condition
theta_0=[0;3];

%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,theta_0);
%plotting angular displacement and angular velocity w.r.t time
plot(t,results(:,1))
hold on
plot(t,results(:,2))
axis([0 20 -3.5 3.5]);
xlabel('time');
ylabel('angular dispalcement & angular velocity');
grid on
title('damping effect in simple pendulum');

figure(1);
ct=1;
for i=1:length(t)
    %initial position of the pendulum
    x0=0;
    y0=0;
    %final position of the pendulum
    x1=1*sin(results(i,1));
    y1=-1*cos(results(i,1));
    figure(2);
    %simulation of motion of pendulum
    %base line
    plot([-1 1],[0 0],'linewidth', 3 ,'color','g');
    %pendulum line
    plot([x0 x1],[y0 y1],[-0.5 0.5],[0 0],'linewidth',2,'color','b');
    hold on;
    %pendulum ball
    plot(x1,y1,'marker','o','markersize',25,'markerFacecolor','r');
    axis([-1.5 1.5 -1.5 0.5]);
    title('motion of a simple pendulum with damping');
    pause(0.1);
    hold off;
    M(ct)=getframe(gcf);
    ct=ct+1;
end
%animation of motion of simle pendulum
movie(M)
videofile=VideoWriter('simple_pendulum.avi','Uncompressed AVI');
open(videofile);
writeVideo(videofile,M);
close(videofile);

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

• ‘clear all’ command is used to clear the data stored in work space.
• ‘close all’ command is used to close all plots and figures.
• ‘clc’ command is used to clear the command window.
• 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’.
• The time span of the motion of simple pendulum ‘t_span’ is taken to be 20 seconds with 500 equal intervals.
• 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,m,l)
  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-is an inbuilt Matlab function for ODE solvers, it compares 4^th and 5^th order methods to estimate error and determine step size. For differential equations with smooth solutions, ode45 is often more accurate then ode23.
• [t,results] contains the results which we need to plot w.r.t time in which results(1) contains angular displacement values in radian and results(2) contains angular velocity values in radian/sec. By using this we can calculate different values of angular displacement and angular velocity w.r.t time. Thus 'results' is an array.
• Using these results , the graph of angular displacement and angular velocity w.r.t time is plotted.
• Then counter'ct' is initialized . Here 'ct' is a variable that simply increments by a specific step value in every iteration of loop.
• To provide continuous motion to the pendulum , a 'fpr' 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,xy0,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 the circle or bob of pendulum and the counter is incremented by 1.
• The 'pause' command is used to achieve better animation speed.
• Then the 'getframe' command is used to capture the motion of lines and stored in variable M.
• By using 'movie' command , Matlab is going to make a movie with all the frames stored in variable M.
• 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.
• Then the video file is opened after the plot , the movie is run and after completion the video stops and closes.

RESULTS:-

1. As seen in thr graph, the periodic curve of Angular displacement and Angular velocity reduces w.r.t time . This 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 animation of motion of the simple pendulum , it can be concluded that, as the damping coefficient 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.
3. The link for the animation of simple pendulum coded in the above matlab software is given below.
4. https://youtu.be/FRRmBJxmaYE