Week 4.1 - Solving second order ODEs

Solving Second order ODE’s-Simple Harmonic Motion of Pendulum

Aim:

To solve an ordinary Differential Equation which represents the equation of motion of simple pendulum with damping and to simulate the motion of the pendulum using MATLAB.

Objective:

To write a program that solves the following ODE. This ODE represents the equation of motion of a simple pendulum with damping

`(d^2theta)/(dt^2)+b/m(d theta)/dt+g/lsin theta =0`

Introduction:

In Engineering, ODE is used to describe the transient behavior of a system. A simple example is a pendulum

The way the pendulum moves depends on the Newtons second law. When this law is written down, we get a second order Ordinary Differential Equation that describes the position of the "ball" w.r.t time.

Differential Equation:

A Differential equation is the one involved with one or more derivatives including the study of changes.

A derivative or slope is the ratio of change in ‘Y’ to change in ‘X’ respectively.

Differentiating the function y=f(x), leads to the ‘First order derivatives.

dy/dx

Types of Differential Equations:

1. Ordinary Differential Equation (ODE)
2. Partial Differential Equation (PDE)

Application of Differential Equation:

• Used in the study of changes where it helps in analyse the dynamic problems involving flow studies and behavioural studies
• Helps in computational simulations of real-world problems
• Highly applied in fluid dynamics, structural analysis & thermal studies.

Solving method of Differential Equations:

1. Taylor Series:

The Taylor series of a function is an infinite sum of terms that are expressed in terms of the function’s derivatives at single point.

It will be used in numerical approximations of any function.

2. Runge Kutta method (Rk):

RK method is one of the most familiar of approximation which is an extended Euler method derived from the Taylor series

It overcomes the demerits of excluding the calculations of higher order derivatives in turns. It can be used for any functions having high accuracies.

• 1^st order RK / Euler Method

Euler method is nothing but the first order RK method involving only the first order derivatives

It is formed from the Taylor series as follows:

y_Euler = y_k + ∆x f(x_k , y_k )

• 4^th order RK Method:

ODE45 uses an RK 4^th order method to perform numerical integration at each point to find the solution

`dy/dx = f(x,y)`

With initial condition (x0, y(x0))

4^th order RK Method

`y(x0+h) = y(x0)+1/theta(F1+2F2+2F3+F4)`

Where,

`F1=hf(x0,y0)`

`F2=hf(x0+h/2,y0+(F1)/2)`

`F3=hf(x0+h/2,y0+(F2)/2)`

`F4=hf(x0+h,y0+F3)`

Ode45 Syntax:

[t, y] = ode45(odefun, t_span, y0)

Approach Flow:

1. Discretisation of the ODE:

`(d^2theta)/(dt^2)+b/m(d theta)/dt+g/lsin theta =0`

Displacement,

`theta=theta1`

Velocity, `(d theta)/dt=theta2` ------ 1

`(d^2theta)/dt^2=-b/(ml)theta2-g/lsintheta1=0`------ 2

2. Initialising the Variables:

Damping coefficient (b) = 0.05

Acceleration due to gravity (g) =9.81m/s²

Length of string (l) = 1m

Mass of Bob (m) = 1kg

Initial Displacement = 0rad

Initial Velocity = 3rad/s

Time span = 0 to 20 sec.

Program Code explanation:

1. Input of different parameters (damping Coefficient, mass, Acceleration due to gravity and length of the string of pendulum) are taken.

% Inputs

% Damping coefficient(b)

b = 0.05;

% gravity in m/s^2 (g)

g = 9.81;

% Length of pendulum in meters(l)

l = 1;

% Mass of Bob in kg (m)

m = 1;

2. Initial conditions are defined. At time t=0, angular displacement = 0, angular velocity = 3 rad/sec.

% Initial condition

theta_0 = [0;3];

3. Input for time span or period over which pendulum oscillates is taken.

% Time points

t_span = linspace(0,20,500);

4. Separate program is written by using ‘function’ in MATLAB to calculate dtheta/dt, where is column vector given by theta = [theta1; theta2]. Here theta1 represents angular displacement and theta2 represents angular velocity. After doing so we have values of dtheta1/dtand dtheta2/dt.

% Function defining the ODE

function [dtheta_dt] = ODE_function(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

5. Now to calculate theta1 and theta2, ‘ode45’ command is used in MATLAB.

By using ode45 command, dtheta1/dt` and dtheta2/dt is solved to get angular displacement and angular velocity. Here 'results' is an array with 2 columns and rows equal to number of time values (here 500).

% Solve ODE

[t, results] = ode45(@(t,theta)ODE_function(t,theta,b,g,l,m),t_span,theta_0);

6. Pivot point of pendulum is defined (x0, y0) = (0, 0) and plotted on x-y axis

To extract the first column of a result and assigned it to 'Displacement'

% Plot the Results

plot(t, results(:,1))

hold on

plot(t, results(:,2))

ylabel('PLOT')

xlabel('TIME')

legend('Displacement','Angular velocity')

grid on

hold off

displacement = results(:,1);

velocity=results(:,2);

%fixed coordinates

x0=1;

y0=1;

7. To create images for each time value and stored it in array ‘M’.

For loop is used from i=1: length(t_span) where t_span is no of times the angular displacement and angular velocity is to be calculated from time t = 0 to 20.

x and y coordinate of mass of a pendulum is calculated for each time as x1 and y1

Line is created from (x0,y0) to (x1,y1) which represents string of a pendulum.

Bob or mass of a pendulum is created by using ‘marker’ command with suitable size of marker and colour.

Tittle is given to image with 'title' command.

'hold off' command is used to erase image from previous time step.

Images are captured for each time step and stored in array ‘M’.

ct=1;

for i = 1:length(displacement)

%coordinates

x1=1-l*sin(displacement(i));

y1=1-l*cos(displacement(i));

figure(2)

%plot fixed link

plot([0 2],[1 1],'linewidth',5,'color','b')

hold on

plot([1 x1],[1 y1],[x1 x1],[y1 y1], 'O','markersize',10,'markerfacecolor','r','linewidth',4)

plot(x0,y0,'-O','markersize',3,'markerfacecolor','b')

title({'Simple harmonic motion'; 'of a pendulum'}, 'color','r')

%axis setting

axis([0 2 -0.5 2])

pause(0.03)

% to change the frame of the movie

M(ct)= getframe(gcf);

ct=ct+1;

hold off

end

8. Animation is created using 'movie(M)' command:

movie(M)

videofile=VideoWriter('simple pendulum.avi','uncompressed AVI');

open(videofile)

writeVideo(videofile,M)

close(videofile)

Main Program Code:

Part 1 : Defining function:

% Function defining the ODE
function [dtheta_dt] = ODE_function(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

Part 2 : Main code:

% Solving the second order ODE & Show the Animation
clear all
close all
clc
% Inputs
% Damping coefficient(b)
b = 0.05;
% gravity in m/s^2 (g)
g = 9.81;
% Length of pendulum in meters(l)
l = 1;
% Mass of Bob in kg (m)
m = 1;

% Initial condition
theta_0 = [0;3];

% Time points
t_span = linspace(0,20,500);

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

% Plot the Results
plot(t, results(:,1))
hold on
plot(t, results(:,2))
ylabel('PLOT')
xlabel('TIME')
legend('Displacement','Angular velocity')
grid on
hold off

displacement = results(:,1);
velocity=results(:,2);

%fixed coordinates
x0=1;
y0=1;
ct=1;

for i = 1:length(displacement)
    
    %coordinates
    x1=1-l*sin(displacement(i));
    y1=1-l*cos(displacement(i));
    
    figure(2)
     %plot fixed link
     plot([0 2],[1 1],'linewidth',5,'color','b')
     hold on                              
     plot([1 x1],[1 y1],[x1 x1],[y1 y1],'-O','markersize',10,'markerfacecolor','r','linewidth',4)
     plot(x0,y0,'-O','markersize',3,'markerfacecolor','b')
     title({'Simple Harmonic Motion'; 'of a Pendulum'}, 'color','r')
     %axis setting
     axis([0 2 -0.5 2])
     pause(0.03)
     % to change the frame of the movie
     M(ct)= getframe(gcf);                 
     ct=ct+1;
     hold off
end

movie(M)                                                                
videofile=VideoWriter('simple pendulum.avi','uncompressed AVI');
open(videofile)                                                         
writeVideo(videofile,M)                                                 
close(videofile)

Results:

1. Plot for Displacement & Angular Velocity:

2. Animation video link uploaded

https://drive.google.com/file/d/1BX5k3xBV8R9PUjQBVlhG9jn5BWK8Wcaj/view?usp=sharing

Conclusion:

The reduction in velocity and Angle displacement is due to damping coefficient introduced in the Second order ODE of a simple pendulum. Governing differential equation of simple pendulum is solved using MATLAB program. Animation is working correctly. Hence it can be concluded that solution obtained by program is correct.