Solving second order ODEs

AIM:

To write a program to solve the ODE and simulate the motion of a simple pendulum with damping in MATLAB

THEORY:

In Engineering, ODE (ordinary differential equation) is used to describe the transient behavior of a system. A simple example is a pendulum. The way the pendulum moves depends on the Newton’s second law. When this law is written down, we get a second order ODE that describes the position of the ball with respect to time.

The formula for Second order ODE

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

Where, g = gravity in `m/s^2` ; l = length of the pendulum in m ; m = mass of the ball in kg ; b = damping coefficient

PROCEDURE:

FOR MAIN PROGRAM

• Start programming in MATLAB with these three commands, as all the previous calculations are erased.

• clear: 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.

• Now start with the inputs such as the damping coefficient, acceleration due to gravity, length of the pendulum, and mass of the ball.
• Now define the initial conditions of angular displacement and the angular velocity in a variable “theta_0”.
• Now with the variable “t_span” time span of the motion of the pendulum is defined using the “linspace” command which divides the range of 0 to 20 in 500 equal intervals.

ODE SOLVER

• Now the ODE is solved using the ODE solver that is inbuilt in Matlab i.e. “ode45”. The syntax of “ode45” is first we define the independent and the dependent variable “(t,theta)” then the function is called that contains the ODE, this function is integrated for getting the results “ode_function(t,theta,b,g,l,m)” and then the time span is given for which the function is integrated “t_span” and then lastly the initial condition is given from which the calculation starts “theta_0”. Any ODE to be integrated by “ode45” must be provided as a set of first order ODEs. All the results are stored in an array “[t, results]”.

PLOTTING DISPLACEMENT AND VELOCITY Vs TIME

• Now for plotting the results, I have used the “subplot” command which divides the figure window in the way we want. Consider the figure window as a matrix and we can select the size of the matrix and allot the position of the plot in that matrix, here it is set as 5*3 where 5 is the number of rows and 3 sis the number of the columns. A 5*3 matrix means 15 boxes are present therefore boxes 10 to 15 is assigned to this plot.
• Now using the “plot” command first graph is plotted between “t” and “results” in which data of only first column is used “results(:,1)” and by using the “hold on” command second graph between “t” and “results” containing the data of second column “results(:,2)” is plotted on the same graph and then “linewidth”’ command is used to set the width of the line and then through “color” color of the line is assigned.
• Now the range of the axis of the plot is defined, the x ranges first and then the y ranges “axis([0 2 -5.5 5.5])”
• Now labelling of x axis and y axis is done respectively, “xlabel(‘Time’)” and “ylabel(‘Plot’)” and then the legend and grid is switched on “grid on”

ANIMATION CREATING

• Now for creating an animation “for” loop is used through which only one data at a time is stored. Now the initial and final position of the pendulum is defined. For plotting the animation “subplot” command is used and boxes 1 to 9 is assigned to it.
• Now with “plot” command base at which the pendulum is attached is plotted. And then through “line” command string and then again through “plot” command the ball of the pendulum is plotted. Now the grid is switched on “grid on” and also the range of the axis is defined. The compiling of the code is paused for 0.1.
• Now for creating an animation a variable “M” is defined which stores all the frames of an iteration by using the command “getframe” and then for the increment of the frame “c+1” is done.
• Now for creating a movie “movie” command is used. Then for writing the video, giving a name and format to the video “VideoWriter” is used.
• Now “open(videofile)” opens the file and then through “writeVideo(videofile,M)” video is written in it and it is closed through “close(videofile)”.
• Now run the program.

FOR FUCNTION PROGRAM

• The name of the function and the name of the file should be same otherwise the program will not run. Here it is “ode_function”.
• Now the 2^nd order ODE is broken into two 1^st order ODEs, first one to be “dtheta1_dt” and the second to be “dtheta2_dt”
• Now the two results are combined in an array “dtheta_dt” which contains values in two rows.

CODE:

FOR MAIN PROGRAM

% to create an animation of pendulam movement by solving a 2nd order ODE
clear
close all
clc

%inputs
b=0.05; %damping coefficient
g=9.81; %acceleration due to gravity in 'm/s^2
l=1; %length of the pendulum in 'm'
m=0.1; %mass of the pendulum ball in 'kg'

%initial condition - angular displacement '0' and angular velocity '3 rad/s' 
theta_0 = [0;3];

%time span of the motion of pendulum 
t_span = linspace(0,20,500);

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

%plotting the values using subplot
subplot(5,3,[10:15])
%figure(1)
plot(t,results(:,1),'linewidth',2,'color','r')
hold on 
plot(t,results(:,2),'linewidth',2,'color','b')
axis([0 20.2 -5.5 5.5])
xlabel('Time')
ylabel('Plot')
legend('Angular Displacement','Angular Velocity')
grid on
hold off

%animation of pendulum movement
c=1;
for i=1:length(results(:,1))
    x0=0; %initial position
    y0=0; %initial position
    x1=l*sin(results(i,1)); %final position
    y1=-l*cos(results(i,1)); %final position
    
    %plotting the animation 
    subplot(5,3,[1:9])
    %figure(2)
    plot([-l l],[0 0],'linewidth',4,'color','k') %fixed base where the pendulum is attached
    hold on
    line([x0 x1],[y0 y1],'linewidth',2,'color','b') %string of the pendulum
    plot(x1,y1,'o','markersize',20,'markerfacecolor','r') %ball of the pendulum
    grid on
    axis ([-1.2 1.2 -1.2 0.2])
    pause (0.1)
    hold off
    M(c)=getframe(gcf);
    c=c+1; %increment of frame   
end
movie(M)
videofile=VideoWriter('Pendulum_Motion.avi','Uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)

FOR FUNCTION

function [dtheta_dt] = ode_function(t,theta,b,g,l,m)
%second order ODE is broken into two firstorder ODEs
theta1 = theta(1);
theta2 = theta(2); 
dtheta1_dt = theta2; %first ODE 
dtheta2_dt = -(b/m)*theta2 - (g/l)*sin(theta1); %second ODE
dtheta_dt = [dtheta1_dt; dtheta2_dt]; % combining the two results in an array in two rows
end

RESULTS:

ANIMATION: Pendulum motion

CONCLUSION:

The second order ODE is solved using “ode45” command and the motion of the pendulum is also simulated.