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
% the length of the x-vector.
% OUTPUTS
% Integral of y dx from x(1) to x(n)
% 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 is the number of cubics
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