Week 3 - Solving second order ODEs

Solving second order ODE's for a Simple Pendulum

AIM:

TO STUDY TRANSIENT BEHAVIOR OF A PENDULUM BY SOLVING SECOND ORDER ORDINARY DIFFERENTIAL EQUATION AND ALSO SIMULATE THE PENDULUM MOTION WITH THE HELP OF MATLAB PROGRAM.

THEORY:

A simple pendulum consists of a mass m hanging from a string of length L and fixed at a pivot point P. When displaced to an initial angle and released, the pendulum will swing back and forth with periodic motion.

GOVERNING EQUATIONS:

The way the Simple 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.

`(d^2 θ)/(dt^2 )+ b/m dθ/dt+ g/L sinθ=0` -------(1)

In the above equation,

g = gravity in m/s²,

L= length of the pendulum in m,

m = mass of the ball in kg,

b=damping coefficient.

To solve this 2^nd order ODE we have to reduce it into lower order ODE

Let `θ=θ_1` -----------(2)

So we can say that,`dθ/dt=〖dθ〗_1/dt=θ_2` -----------(3)

Therefore,

`(d^2 θ)/(dt^2 )=(d^2 θ_1)/(dt^2 )=d/dt ((dθ_1)/dt)=d/dt(θ_2)` -------------(4)

By substituting (2),(3),(4) in (1) we get

`(dθ_2)/dt+ b/m θ_2+g/L sinθ_1=0`

After rearranging this equation we get

`(dθ_2)/dt= -(b/m θ_2 )-g/L sinθ_1`

`〖dθ〗_1/dt=θ_2`

We can also write these two differential equation into matrix form,

`d/dt [■(θ_1@θ_2 )]=[■(θ_2@-b/m θ_2-(g )/L sinθ_1 )]`

OBJECTIVES:

1. To study the transient behavior of the simple pendulum system with damping by solving second order ODE.
2. Animate the motion of the simple pendulum.

INPUT VALUES:

L = 1 m

m = 1 kg

b = 0.05

g = 9.81 m/s²

FUNCTION PROGRAM:

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

theta1=theta(1)
theta2=theta(2)

%second order ODE reduced into single order ODE
dtheta1_dt=theta2;
dtheta2_dt=-(b/m)*theta2-(g/l)*sin(theta1);

%To make an Array of two single order ODE
dtheta_dt=[dtheta1_dt;dtheta2_dt];
end

MAIN PROGRAM:

clear all 
close all
clc

%Inputs
l=1;    %in m
b=0.05;
g=9.81; % in m/s^2
m=1;    % in kg

% Initial conditions 
theta_0=[0,3];
%time
t_span=linspace(0,20,500);
%solve ode
[t,results]=ode45(@(t,theta)pendulum_damped_ode(t,theta,b,g,l,m),t_span,theta_0);
figure(1)
plot(t,results(:,1),'color','b','linewidth',2);
hold on
plot(t,results(:,2),'color','r','linewidth',2);
hold off

xlabel('Time');
ylabel('Plot');

ct=1
for i=1:length(results(:,1));
    
    % Starting point coordinates
    x0=0;
    y0=0;
    
    % End point coordinates
    x1=l*sin(results(i,1));
    y1=-l*cos(results(i,1));
    
    figure(2)
    %To plot ceiling from which pendulum hangs
    plot([-1 1],[0 0],'linewidth',4,'color','k');
    hold on
    %To plot massless thread
    line([x0 x1],[y0 y1],'linewidth',2,'color','b');
    hold on
    
    % To represent Bob  
    plot(x1,y1,'-o','markers',20,'markerfacecolor','r');
    grid on
    %Limit of x-axis and y-axis
    axis([-1.5 1.5 -1.5 0.5])
    pause(0.003);
    hold off
    % To get figure window and save it as an array
    q(ct)=getframe(gcf);
    
    ct=ct+1;
end
% creating Movie of all these frame
movie(q)
% To dump out all these frame and saving the movie in ".avi" format 
videofile=VideoWriter('pendulum_damped.avi','uncompressed AVI');
open(videofile)
% Writing the movie in the video file
writeVideo(videofile,q)
close(videofile)

EXPLANATION OF PROGRAM:

1. Function Program

• In this program we have create a function covering all the parameters and variables. Remember that the name of the program should be saved with the function name.
• Reduce second order ODE into single order ODE.
• Now make an Array of two single order ODE.
• And at last end the function program by using end command.

2. Main Program

• First we give the initial value for the system.
• Then specify time span over which the pendulum’s transient behavior has to be observed .
• Now, the ode45 solver is used to calculate the results by calling function inside it.
• Further, we plot the Displacement v/s Time and Velocity v/s Time graphs on the same figure to study the respective behaviours by using plot command
• Now Plot ceiling , massless thread and by using markers plot bob inside for loop and get the figure window and store it in an array ‘q’.
• And then terminate the loop.
• movie command is used to create movie/animation by combining all the figures saved in ‘q’ file.
• And lastly save this animation as .avi.

OUTPUT :

OUTPUT PLOT:

ANIMATION:

The link for the animation has been mentioned below.

Motion of a Simple Pendulum with damping

CONCLUSION:

We have successfully solved the second order ODE and studied the transient behavior of 2^nd order ODE .

And also animate the motion of simple pendulum.