Integrating functions given at discrete points via MATLAB

When integrating functions given at discrete data points in MATLAB, trapz is the function of choice to go for. But we can get more accurate results by interpolating via cubic splines with the MATLAB spline function.  Since the spline function is made of piecewise cubics, Simpson’s 1/3rd rule can be used to exactly integrate them.   Here is a test program and a function written for you.  You can download the mfile here and a published version here for more readability if you wish.

clc
clear all
% Author: Autar Kaw, AutarKaw.com
% https://creativecommons.org/licenses/by-nc-sa/4.0/ 
% Testing the program with data given at discrete points
% for y=x^6
xx=[1  1.5    2   2.5     3     3.5    4   4.5   5];
yy=[1  1.5^6  2^6 2.5^6   3^6   3.5^6  4^6 4.5^6 5^6];
n=length(xx);
splineintegval=splineintegral(xx,yy);
fprintf('Value of integral using spline =%g',splineintegval)
% Exact value of integral if function was given continuously 
syms x
exact=vpaintegral(x^6,x,xx(1),xx(n));
fprintf('\n Value of integral using exact integration =%g',exact)

%% Function to integrate via spline interpolation
function splineval=splineintegral(x,y)
% This function integrates functions given at discrete data points
% INPUTS
% The x-values are given in ascending order 
% The limits of integration are x(1) to x(n), where n is teh length of
% the x-vector.
% OUTPUTS
% Integral of y dx from x(1) to x(n)
% Author: Autar Kaw, AutarKaw.com 
% https://creativecommons.org/licenses/by-nc-sa/4.0/ 
% The function finds the mid-point value of y between the given 
% x-values at the mid-point. Then since the spline is made of cubics,
% it uses the Simpson's 1/3rd rule to integrate the cubics exactly
n=length(x);
m=n-1;
% Calculating mid-points
    for i=1:1:m
        xmid(i)=(x(i)+x(i+1))*0.5;
    end
% Calculating value of y at the midpoints
polyvalmid=spline(x,y,xmid);
% Using Simpson's 1/3rd rule of integration to integrate cubics
splineval=0;
    for i=1:1:m
        splineval=splineval+(y(i)+y(i+1)+4*polyvalmid(i))*(x(i+1)-x(i))/6;
    end
end

____________________________

This post is brought to you by

Leave a Reply