How do I integrate a discrete function in MATLAB?

Many students ask me how do I do this or that in MATLAB.  So I thought why not have a small series of my next few blogs do that.  In this blog, I show you how to integrate a discrete function.

The MATLAB program link is here.

The HTML version of the MATLAB program is here.

_____________________________________________________

%% HOW DO I DO THAT IN MATLAB SERIES?
% In this series, I am answering questions that students have asked
% me about MATLAB.  Most of the questions relate to a mathematical
% procedure.

%% TOPIC
% How do I integrate a discrete function?  Three cases of data are
% discussed.

%% SUMMARY

% Language : MATLAB 2008a;
% Authors : Autar Kaw;
% Mfile available at
% http://nm.mathforcollege.com/blog/integrationdiscrete.m;
% Last Revised : April 3, 2009;
% Abstract: This program shows you how to integrate a given discrete function.

clc
clear all

%% INTRODUCTION

disp(‘ABSTRACT’)
disp(‘   This program shows you how to integrate’)
disp(‘   a discrete function’)
disp(‘ ‘)
disp(‘AUTHOR’)
disp(‘   Autar K Kaw of http://autarkaw.wordpress.com’)
disp(‘ ‘)
disp(‘MFILE SOURCE’)
disp(‘   http://nm.mathforcollege.com/blog/integrationdiscrete.m’)
disp(‘ ‘)
disp(‘LAST REVISED’)
disp(‘   April 3, 2009’)
disp(‘ ‘)

%% CASE 1

%% INPUTS

% Integrate the discrete function y from x=1 to 6.5
% with y vs x data given as (1,2), (2,7), (4,16), (6.5,18)
% Defining the x-array
x=[1  2  4  6.5];
% Defining the y-array
y=[2  7  16  18];

%% DISPLAYING INPUTS
disp(‘____________________________________’)
disp(‘CASE#1’)
disp(‘LOWER LIMIT AND UPPER LIMITS OF INTEGRATION MATCH x(1) AND x(LAST)’)
disp(‘ ‘)
disp(‘INPUTS’)
disp(‘The x-data is’)
x
disp(‘The y-data is’)
y
fprintf(‘  Lower limit of integration, a= %g’,x(1))
fprintf(‘\n  Upper limit of integration, b= %g’,x(length(x)))
disp(‘ ‘)

%% THE CODE

intvalue=trapz(x,y);

%% DISPLAYING OUTPUTS

disp(‘OUTPUTS’)
fprintf(‘  Value of integral is = %g’,intvalue)
disp(‘  ‘)
disp(‘___________________________________________’)

%% CASE 2

%% INPUTS

% Integrate the discrete function y from x=3 to 6
% with y vs x data given as (1,2), (2,7), (4,16), (6.5,18)
% Defining the x-array
x=[1  2  4  6.5];
% Defining the y-array
y=[2  7  16  18];
% Lower limit of integration, a
a=3;
% Upper limit of integration, b
b=6;
%% DISPLAYING INPUTS

disp(‘CASE#2’)
disp(‘LOWER LIMIT AND UPPER LIMITS OF INTEGRATION DO not MATCH x(1) AND x(LAST)’)
disp(‘  ‘)
disp(‘INPUTS’)
disp(‘The x-data is’)
x
disp(‘The y-data is’)
y
fprintf(‘  Lower limit of integration, a= %g’,a)
fprintf(‘\n  Upper limit of integration, b= %g’,b)
% Choose how many divisions you want for splining from a to b
n=1000;
fprintf(‘\n  Number of subdivisions used for splining = %g’,n)
disp(‘  ‘)
disp(‘  ‘)

%% THE CODE

xx=a:(b-a)/n:b;
% Using spline to approximate the curve from x(1) to x(last)
yy=spline(x,y,xx);
intvalue=trapz(xx,yy);

%% DISPLAYING OUTPUTS

disp(‘OUTPUTS’)
fprintf(‘  Value of integral is = %g’,intvalue)
disp(‘  ‘)
disp(‘___________________________________________’)
%% CASE 3

%% INPUTS

% Integrate the discrete function y from x=1 to 6.5
% with y vs x data given as (1,2), (4,16), (2,7), (6.5,18)
% The x-data is not in ascending order
% Defining the x-array
x=[1  4   2 6.5];
% Defining the y-array
y=[2  16  7 18];
% Lower limit of integration, a
a=3;
% Upper limit of integration, b
b=6;
%% DISPLAYING INPUTS

disp(‘CASE#3’)
disp(‘LOWER LIMIT AND UPPER LIMITS OF INTEGRATION DO not MATCH x(1) AND x(LAST) ‘)
disp(‘AND X-DATA IS NOT IN ASCENDING OR DESCENDING ORDER’)
disp(‘   ‘)
disp(‘INPUTS’)
disp(‘The x-data is’)
x
disp(‘The y-data is’)
y
fprintf(‘  Lower limit of integration, a= %g’,a)
fprintf(‘\n  Upper limit of integration, b= %g’,b)
% Choose how many divisions you want for splining from a to b
n=1000;
fprintf(‘\n  Number of subdivisions used for splining = %g’,n)
disp(‘  ‘)
disp(‘  ‘)

%% THE CODE
[x,so] = sort(x); % so is the sort order
y = y(so); % y data is now in same order as x data
xx=a:(b-a)/n:b;
% Using spline to approximate the curve from x(1) to x(last)
yy=spline(x,y,xx);
intvalue=trapz(xx,yy);

%% DISPLAYING OUTPUTS

disp(‘OUTPUTS’)
fprintf(‘  Value of integral is = %g’,intvalue)
disp(‘  ‘)

____________________________________________________________

This post is brought to you by Holistic Numerical Methods: Numerical Methods for the STEM undergraduate at http://nm.mathforcollege.com, the textbook on Numerical Methods with Applications available from the lulu storefront, and the YouTube video lectures available at http://nm.mathforcollege.com/videos and http://www.youtube.com/numericalmethodsguy

Subscribe to the blog via a reader or email to stay updated with this blog. Let the information follow you.