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