Week 3 - Solving second order ODEs

AIM: Write a MATLAB program to solve 2nd Order Differential equation and also show an animation of motion of simple pendulum with damping.

OBJECTIVE:

1) To solve a 2nd order differential equation.

2) To create a animation video of simple pendulum under damping.

GOVERNING EQUATION:

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.

Similarly, there are hundreds of ODEs that describe key phenomena and if we solve these equations then we will be able to simulate any of theses systems.

This ODE represents the equation of motion of a simple pendulum with damping

In the above equation,

g = gravity in m/s2,

L = length of the pendulum in m,

m = mass of the ball in kg,

b=damping coefficient.

we use,

L=1 metre,

m=1 kg,

b=0.08.

g=9.81 m/s2.

PROGRAM:

1) Main program for diplacement, velocity time graph and video of simple pendulum under damping

% write a program to stimulate a second order differential equation
close all
clear all
clc


% input geometric parameters
b=0.08;
g=9.81;
l=1;
m=1;

% initial condition
% i.e at theta = 0, displacement = 0 and velocity = 5
theta_0 = [0;5];


%time interval to integrate
t_span = linspace(0,10,100);


% solve ODE
[t,result] = ode45(@(t,theta) ode_func(t,theta,b,g,l,m),t_span,theta_0);

%plotting displacement and time graph
figure(1)
plot(t,result(:,1),'color','r')
xlabel('time')
ylabel('displacement')
axis([0 10 -2 2])
hold on


%plotting velocity and time graph
figure(2)
plot(t,result(:,2))
ylabel('velocity')
xlabel('time')
axis([0 10 -5 5])


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


%for loop to counter all frames
ct=1;

for i=1:length(displacement);
    theta3=displacement(i);
    
   
    x0=0;
    y0=0;
    x=l*sind(theta3);
    y=-l*cosd(theta3);
    
    figure(3)
    
    %plotting fix end
    plot([-0.1 0.1],[0 0],'linewidth',3,'color','k')
    hold on
    
    
    %joining the fix point and end point of line
    plot([x0 x],[y0 y],'linewidth',3,'color','k')
    hold on
  
    
    %plotting the bob
    plot(x,y,'-o','markers',10,'markerfacecolor','r')
    hold off
    
    
    axis([-0.2 0.2 -1.5 0.5])
    pause(0.0000004)
    M(ct)=getframe(gcf);
    ct=ct+1;
   
end

%command to create animation with all frames
movie(M)
Videofile=VideoWriter('pendulum.avi','uncompressed AVI');
open(Videofile)
writeVideo(Videofile,M)
close(Videofile)

2) function program to divide 2nd order differential equation into 2 single order differential equation

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

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


end

OUTPUT:

1) Graph of velocity and displacement with time of under damped system

2) Motion of simple pendulum under damping