Week 4.1 - Solving second order ODEs

IM - To solve second order ordinary differential equation using MATLAB.

Theory -

Ordinary Differential Equations (ODEs) - 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.

Simple Pendulum - The simple pendulum is an idealized mathematical model of a pendulum.This is a weight (or bob) on the end of a massless cord suspended from a pivot, without friction. When given an initial push, it will swing back and forth at a constant amplitude. Real pendulums are subject to friction and air drag, so the amplitude of their swings declines.

The point from which the bob is suspended is known as point of suspension. Weight of the bob acts vertically downwards and the tension acts vertically upwards, so, the bob comes under equillibrium position and the bob starts oscillating if even smaller external force is applied.

Governing Equations -

θθ is the angle made by the movable link with vertical axis.

LL is the length of the movable rod

mm is the mass of bob , bb is damping coefficient of the system.

1. d2θdt2+bm⋅d(θ)dt+gLsinθ=0d2θdt2+bm⋅d(θ)dt+gLsinθ=0 second order ODE -----(1)

To solve this ODE we need to reduce the order of equation into lower order

let θ=θ1θ=θ1

so this implies d(θ)dt=d(θ1)dtd(θ)dt=d(θ1)dt = θ2θ2 ---------(2)

therefore d2θdt2=d2θ1dt2=d(θ2)dtd2θdt2=d2θ1dt2=d(θ2)dt-----(3)

Putting the value of (2) and (3) in (1) , we get-

d2θ1dt+bm⋅d(θ1)dt+gL⋅sinθ1=0d2θ1dt+bm⋅d(θ1)dt+gL⋅sinθ1=0

therefore , d(θ2)dt+bm⋅(θ2)+gL⋅sinθ1=0d(θ2)dt+bm⋅(θ2)+gL⋅sinθ1=0

we got, d(θ2)dt=−bm⋅(θ2)−gLsinθ1d(θ2)dt=-bm⋅(θ2)-gLsinθ1

Objectives -

1. To solve the second order ODE with graph of Velocity v/s Displacement.

2. To create the animation of the simple pendulum.

Procedure -

Given details -

Simulate the motion between 0-20 sec

for angular displacement = 0

angular velocity=3 rad/sec

at time t = 0

L = 1 metre

M = 1 kg

b = 0.05

g = 9.81 ms2ms2

%program for simple pendulum using ode45.
clear all
close all
clc

%input variables for the secound order ode of a simple pendulum
b = 0.05;
g = 9.81;
l = 1;
m = 0.1;

theta_0 = [0; 3];

t_span = linspace(0,20,500);

%ode 45 function to solve second order ode.
[t,results] = ode45(@(t,theta) ode_func(t,theta,b,g,l,m),t_span, theta_0 );

%plotting of angular velocity and displacement
figure(1)
plot(t,results(:,1))
hold on
plot(t,results(:,2))
grid on
title('motion of pendulum')
xlabel('time in sec')
ylabel('angular velocity(rad/s)/displacement(rad)')
legend('displacement','angular velocity')

angle = results(:,1);

%count command for iterations.
ct = 1;

for i = 1 : length(angle)
    omega = angle(i);
    figure(2)
    title('motion of pendulum')
    
    %fixed corordinates of one end of simple pendulum.
    x0 = 0;
    y0 = 0;
    
    %moving coordinates of the bob.
    x1 = l*cos((3*pi/2)+omega);
    y1 = l*sin((3*pi/2)+omega);
    
    %plotting of the bob motion.
    plot([x0 x1],[y0 y1],'linewidth',3,'color','g');
    hold on;
    plot(x1,y1,'o','markers',15,'MarkerFaceColor','k');
    hold on;
    plot([-3 3],[0 0],'linewidth',2,'color','b');
    hold off
    axis([-3 3 -2 2])
    grid on 
    pause(0.1)
    M(ct) = getframe(gcf);
    ct = ct+1;
end

%to create a animation of the simple pendulum
movie(M)
videofile = VideoWriter('simple_pendulum.avi','Uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)
ode function program:

function [dtheta_dt] = ode_func(t,theta,b,g,l,m)

theta1 = theta(1);
theta2 = theta(2);
dtheta1_dt = theta(2);
dtheta2_dt = -(b/m)*theta2 - (g/l)*sin(theta1);
dtheta_dt = [dtheta1_dt; dtheta2_dt];

end
 

Code Explanation -

1.Firstly, the given details are assigned as variables in the code. Then time span of the simple pendulum is set by using linspace command then, theta value is given as a column matrix. After this we use ode45 function and corresponding 'ode_func' function call in main program.

2.Then used plot command to plot graph corresponding to time and velocity,legend command is used for identify each values in the graph.

https://youtu.be/CeoGsPBqUuY3.For getting the accurate simulation of the simple pendulum corresponding to the second order ODE with damping we need to use the looping statements, here for loop is used to create this.

4.The pendulum has fixed coordinates mentioned as x0,y0 and x1, y1 are the other coordinates.

5.Marker command is used for making bob of the pendulum. The axis is fixed using axis command. The counter command is used to make the value changes in each for loop. Then videoWriter command is used to create the video file.