SOLVING SECOND ORDER ODE

SECOND ORDER ORDINARY DIFFERENTIL EQUATIONS

AIM

To write a MATLAB program that solves the second order ordinary differential equation (ODE) representing the motion of a simple pendulum. And, which creates the animation for the same.

INTRODUCTION

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 have the ability to solve these equations then we will be able to simulate any of theses systems. The ODE that we are going to sovle in this particular context also applies to a system consisting a mass of block 'm' attached with a spring and a damper.

GOVERNING EQUATION & TRANSFORMATION

A simple pendulum consists of a mass ‘m’ hanging from a string of length ‘l’ and fixed at a pivot point . When displaced to an initial angle ’θ’ and released, the pendulum will swing back and forth with periodic motion. Also, damping 'b' is provided which would eventually bring the pendulum at rest after a cetain siginificant amount of the time.

The equation is given by:

............................. (i)

For solving any ODE of any order, we first need to break it down into sets of first order ODEs. To do that, let us say,

`theta=theta_1`
`d/dt(theta)=d/dt(theta_1)=theta_2`
`d^2/dt^2(theta)=d^2/dt^2(theta_1)=d/dt(theta_2)`

Thus, we can write equation (i) as,

`d/dt(theta_2)+(b/m)theta_2+(g/l)sintheta_1=0`

Thus,

`d/dt(theta_1)=theta_2` And,
`d/dt(theta_2)=-(b/m)theta_2-(g/l)sintheta_1`

STEPS TO WRITE THE MATLAB PROGRAM

1. Start programming in MATLAB with following three commands.
   clear all: This command clears all the data stored in the workspace previously.
   close all: This command closes all the previous figures and plots.
   clc: This command clears the command window.
   
   
2. Inputs are taken as given in the problem and defined in the program with the appropriate variables a follows:
   l = 1;
   m = 1;
   b = 0.05;
   g = 9.81;
   
   
3. Initial conditions are taken as:
   Angular displacement = 0,
   Angular Velocity = 3 rad/sec, at time t = 0.
   theta = [0;3];
   And, the motion is simulated between 0-20 sec with 500 values.
   t_span = linspace(0,20,500);
   
   
4. Solving ODE: Ode solver ode45 is used in this matlab program. The results are stored in an array of row matrix of two coloumns, for the values of displacement and velocity at different time periods. ode_func is the function made to solve the given equation where the asumptions are made in the from of code.
   [t,results]=ode45(@(t,theta) ode_func(t,theta,b,g,l,m),t_span,theta);
   
   
5. Plotting figure(1)
   Here the graph of time vs angular velocity & time vs angular displacement is plotted to get the wave form of from the initial condtions.
   plot(t,results(:,1))
   hold on
   plot(t,results(:,2));
   xlabel('Time')
   ylabel('Plot')
   
   
6. Creating Loop for Oscillatig the Pendulum & Plotting figure(2)
   As the the bob moves with a velocity and displaces for the different time period, a for loop is created to capture and display all of its postion. Finally, these frames are written down in a videofile as 'pendulum.avi' thereby animating the transient motion of the simple pendulum. The code for this step can be seen in the final code as shown in the preceeding sections of this report.
   
   

FINAL CODE FOR THE MAIN PROGRAM (pendulum.m)

clear all
close all
clc

% INPUTS
l = 1;	% length
m = 1;	% mass
b = 0.05;	% damping coefficient
g = 9.81;	% gravity

% INITIAL CONDITIONS (Angular Displacement = 0, Angular Velocity = 3 rad/sec)
theta = [0;3];

% TIME POINTS (0-20) seconds with 500 values.
t_span = linspace(0,20,500);

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

% PLOTTING (Time Vs Angular Velocity & Displacement Graph)
figure(1)
plot(t,results(:,1))
hold on
plot(t,results(:,2));
xlabel('Time')
ylabel('Plot')

% LOOP FOR OSCILLATING PENDULUM
ct=1;
for i = 1:length(results(:,1))
    
    % Initial boundary conditon
    x0 = 0;
    y0 = 0;
    
    % Defining x and y coordinates
    x1 =  l*sin(results(i,1));
    y1 = -l*cos(results(i,1));
    
    figure(2)
    plot([-0.5 0.5],[0 0],'linewidth',2,'color','black')
    hold on
    plot([0 0],[-1.5 0.5],'LineStyle','--','color','black')
    plot([x0 x1],[y0 y1],'linewidth',1,'color','b')
    plot(x0,y0,'o','markersize',4,'markerfacecolor',[0.6 0.6 0.6]);
    plot(x1, y1,'o','markersize',20,'markerfacecolor','r')
    axis([-1.5 1.5 -1.5 0.5])
    pause(0.01)
    M(ct)=getframe(gcf)
    ct=ct+1;
    clf()   % Clearing the current graphic image
    
end

% Creating a Movie / Animation
movie(M);
videofile = VideoWriter('pendulum.avi','Uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)

CODE FOR THE FUNCTION (ode_func.m)

function [dtheta_dt]= ode_func(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

ERRORS FACED

1. Not using the hold on & hold off or clc() command at the proper step.
2. Determining the x and y values.

Steps To Overcome those errors:

1. Make sure the use of hold on and hold off or clc() commands after the proper plot command or the graph wont be as per requirement and to avoid the path trace to left in the graph.
2. In the motion of the simple pendulum the coridnates x, y are interchanges like: X = l*sinθ and Y = -l*cosθ;

Screenshot of the Plots when the hold on command was missing from the code.

FINAL PLOTS & RESULTS

figure(1): time vs angular velocity & time vs angular displacement plot:

Animation: transient motion of the simple pendulum: