Week 3 - Solving second order ODEs

We have to solve given second order ordinary differential equation (ODE) in matlab. ODE is a differential equation containing one or more functions of independent variable and the derivatives of that functions.  Our objective is to write a program that solves the ODE which represents the equation of motion of a simple pendulum with damping. Simple pendulum is an idealized mathematical model which is nothing but the weight or bob on the end of massless cord suspended from a pivot without friction.

Given details are -

g (gravity) = 9.81 m/s^2,

l (length of the pendulum) = 1 m,

m (mass of the ball) = 1 kg,

b (damping coefficient) = 0.05,

Simulate the motion between 0-20 sec,

for angular displacement=0,

angular velocity=3 rad/sec at time t=0.

Now lets write the programme for ODE which represents the equation of motion of a simple pendulum with damping.

1. Function code of ode45 for main programme -

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

2. Main programme for ODE solving -

b = 0.05;
g = 9.81;
l = 1;
m = 1;

% initial conditions are
theta_0= [0;3];

% time points 
t_span = linspace(0,10,500);

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

figure(1)
plot(t,results(:,1))
hold on
plot(t,results(:,2))
title('motion of pendulum')
xlabel('time(s)')
ylabel('angular velocity')
grid on
hold off


disp = results(:,1);
velocity = results(:,2);
ct = 1;
for i= 1:length(disp)
    
    x0=0; 
    y0=0;
    x1= 1-l*sin(disp(i));
    y1= 1-l*cos(disp(i));
    
    figure(2)
    plot([0 2],[1 1],'linewidth',3,'color','b');
    hold on
    plot([1 x1],[1 y1],[x1 x1],[y1 y1],'-o','markersize',20,'markerfacecolor','y');
    plot(x0,y0,'-o','markersize',3,'markerfacecolor','b');
    hold off
    axis([0 2 -0.5 2])
    grid on 
    
    K(ct) = getframe(gcf);
    ct =ct+1;
end

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

The given details are assigned as variables in the above programme. We have create the ode_function which is call in main programme to used in ode45 command.

Code explanation -

• linspace - To create an array of time span values.
• ode45 - To solve differential equation using ode_function.
• plot - Plot the graph of corrensponding time and velocity.
• legend - To add legend to axes.
• axis - To fix the x axis coordinate and y axis coordinate. 
• xlabel/ylabel - To label the axis.
• title - To give title to the plot.
• hold on/off - To retains plots in the current axes so that new plots added to the axes do not delete existing plots.
• grid on/off - To display or hide axes grid lines.
• for loop - For accurate simulation of simple pendulum corresponding to second oreder ODE we use looping statement.
• counter - To meke value change in each for loop.
• getframe - To capture the figures as movie frame.
• movie - To play recorded movie frames.
• VideoWriter - Creat video objects to write video with name and video profile.
• writeVideo - To write video data with mention file.
• open and close - To open and close file for writing video data

output -

Graph of velocity vs time shows that displacement and angular velocity is decreasing by the increasing time due to damper intorduction.

Animation for ODE represents the equation of motion of a simple pendulum with damping is created.

Animation link -

https://youtu.be/NCvfngwW2Vw