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.

Contents

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://numericalmethods.eng.usf.edu/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://numericalmethods.eng.usf.edu/blog/integrationdiscrete.m')
disp(' ')
disp('LAST REVISED')
disp('   April 3, 2009')
disp(' ')

ABSTRACT
   This program shows you how to integrate
   a discrete function
 
AUTHOR
   Autar K Kaw of http://autarkaw.wordpress.com
 
MFILE SOURCE
   http://numericalmethods.eng.usf.edu/blog/integrationdiscrete.m
 
LAST REVISED
   April 3, 2009
 

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(' ')

____________________________________
CASE#1
LOWER LIMIT AND UPPER LIMITS OF INTEGRATION MATCH x(1) AND x(LAST)
 
INPUTS
The x-data is

x =

    1.0000    2.0000    4.0000    6.5000

The y-data is

y =

     2     7    16    18

  Lower limit of integration, a= 1
  Upper limit of integration, b= 6.5 

THE CODE

intvalue=trapz(x,y);

DISPLAYING OUTPUTS

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

OUTPUTS
  Value of integral is = 70  
___________________________________________

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('  ')

CASE#2
LOWER LIMIT AND UPPER LIMITS OF INTEGRATION DO not MATCH x(1) AND x(LAST)
  
INPUTS
The x-data is

x =

    1.0000    2.0000    4.0000    6.5000

The y-data is

y =

     2     7    16    18

  Lower limit of integration, a= 3
  Upper limit of integration, b= 6
  Number of subdivisions used for splining = 1000  
  

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('___________________________________________')

OUTPUTS
  Value of integral is = 50.4424  
___________________________________________

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('  ')

CASE#3
LOWER LIMIT AND UPPER LIMITS OF INTEGRATION DO not MATCH x(1) AND x(LAST) 
AND X-DATA IS NOT IN ASCENDING OR DESCENDING ORDER
   
INPUTS
The x-data is

x =

    1.0000    4.0000    2.0000    6.5000

The y-data is

y =

     2    16     7    18

  Lower limit of integration, a= 3
  Upper limit of integration, b= 6
  Number of subdivisions used for splining = 1000  
  

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('  ')

OUTPUTS
  Value of integral is = 50.4424  

Published with MATLAB® 7.6