Week 3 - Solving second order ODEs

Aim : To solve Second order ordinary differential equation using MATLAB.

Theory :

An ordinary differential equation consists of one independent variabe and and the derivatives of those functions .

Simple pendulum : This is a mathematical model of a pendulum . There is a weight on the end of a massless cord suspended from a pivot , without friction. When providing an initial push , it will swig back and forth at a constant amplitude . Whereas , real pendulum are subjected to friction and air drag, so the amplitude of the swing decline due to this .

Governing equation :

g = gravity in m/s2,

L = length of the pendulum in m,

m = mass of the ball in kg,

b=damping coefficient.

Theta = Angle made by movable link with vertical axis .

Second order ODE

Objectives :

1. To solve the second order ODE and represent the graph of Velocity Vs dispacement .

2. To create the animation of the simple pendulum.

Procedure :

Given detais :

Simulate the motion between 0-20 sec

for angular displacement=0

angular velocity=3 rad/sec

at time t=0 second

L=1 metre,

m=1 kg,

b=0.05

g=9.81 m/s2.

Matlab codes :

1. Function of ODE45 for main program

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

t=linspace(0,20,600);

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

2. Main Program

* Program to represent the graph between angular dispacement vs angular velocity.

% Objective : Solving 2nd order ODEs of a pendulum
clear all
close all
clc

b=0.05 % Damping coefficient
%%

g=9.81 % gravity (m/s^2)
l=1 % Length of the spring (m)
m=1 % Mass of the bob (kg)

% Initial condition
% angular dispalcement is 0 and angular velocity is 3
theta_O=[0;3];

% Time span 
% Pendulum oscillates between 0-20 seconds

t_span=linspace(0,20,600);

% Solving ODE
[t,results]=ode45(@(t,theta) ode_func(t,theta,b,g,l,m),t_span,theta_O);

% Plotting

figure(1)
plot(t,results(:,1))
hold on
plot(t,results(:,2))

title('Damped Harmonic motion of a pendulum')
legend('Angular displacement','Angular velocity')
ylabel('plot')
xlabel('Time(seconds)')

* Program to create the animation of simple pendulum.

% Objective : Solving 2nd order ODEs of a pendulum
clear all
close all
clc

b=0.05 % Damping coefficient
%%

g=9.81 % gravity (m/s^2)
l=1 % Length of the spring (m)
m=1 % Mass of the bob (kg)

% Initial condition
% angular dispalcement is 0 and angular velocity is 3
theta_O=[0;3];

% Time span 
% Pendulum oscillates between 0-20 seconds

t_span=linspace(0,20,600);

% Solving ODE
[t,results]=ode45(@(t,theta) ode_func(t,theta,b,g,l,m),t_span,theta_O);



% For loop
% Before starting the for,loop, a counter value was defined 

ct=1;
for i=1:length(results(:,1))
    x_start=0
    y_start=0
    x_end=1*sin(results(i,1))
    y_end=-1*cos(results(i,1))

    %Plotting the animation
    figure(2)
    plot([-1 1],[0 0],'linewidth',5)
    hold on
    plot([x_start x_end ],[y_start y_end],'linewidth',3)
    hold on

    plot(x_end,y_end,'o','markersize',20,'MarkerFaceColor','y')
    title('Animation of transient behaviour of a pendulum')
    axis([-1.5 1.5 -1.5 1])
    pause(0.05)
    hold off 
    M(ct)=getframe(gcf)
    ct=ct+1
end 

% Animation of the pendulum
Movie(M)
videofile=VideoWriter('Pendulum_Animation_avi', 'Uncompressed_AVI')
open (videofile)
writeVideo(videofile,M)
close(videofile)

Code explaination :

1. Firstly, The given details are assign as variable 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' to call main program .

2. Then used plot command to plot graph corresponding to time and velocity. Legend command is used for identify different parameters for comparision in graph.

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

4. The pendulum has fixed co ordinates known as x0,y0 & x1,y1 are the other co ordinates .

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 .

Results:

Reference Screeen shots:

1. Displacement vs velocity plot .

2. Animation of simple pendulum.

3. Animation file link

https://drive.google.com/file/d/14YLRendX753STpUtD3-_IT-QRh9ECE2F/view?usp=sharing

Conclusion :

1.The graph of the Displacement and angular velocity is plotted .

2. Animation of simple pendulum is created .