% test_integration

% test_integration.m

% This MATLAB program tests integration routines for the

% integral of exp(-x^2) over [0,1].

% K. Beers.

% MIT ChE. 8/7/2002

function iflag_main = test_integration();

iflag_main = 0;

% Set the limits of integration.

a = 0;

b = 1;

% Set a character string to define the function

% to be integrated. We use .^ so that the function

% is overloaded - we can hand it multiple values of

% x at a time.

fun = 'exp(-x.^2)';

% We set the "true" value, known from the error function.

integral_true = sqrt(pi)/2*erf(1);

disp(['Integral (true value) = ', ...

num2str(integral_true)]);

% trapz() ---------------------------------------------------------

% First, we demonstrate the use of the composite trapezoid rule.

% Make a vector of grid points at which to evaluate the function.

% In this example, we use evenly-spaced grid points; however, if

% you know that your function changes most rapidly in a particular

% region of the domain [a,b], you make the grid more fine in this

% region to provide higher accuracy with the same total number

% of grid points.

resolution = 10; % number of grid points

x_val = linspace(a,b,resolution); % x-values at grid points

% Compute the value of f from the function by making an inline function

F = inline(fun); % we declare an inline function F

f_val = F(x_val); % calculate values of f(x) at grid points

% We plot the function over [a,b].

figure;

plot(x_val,f_val);

xlabel('x');

ylabel('f(x)');

title(['f(x) = ', fun]);

% Then, we estimate the value of the integral using the

% composite trapezoid rule. To increase the accuracy of

% the estimate, we increase the resolution.

integral_trap = trapz(x_val,f_val);

disp(['Integral (composite trapezoid) = ', ...

num2str(integral_trap)]);

% quad() ------------------------------------------------------------

% Next, we use the built-in function quad() to estimate the integral.

% This routine requires the user to supply a function that returns the

% value of f(x) given an input value of x. For more detail on the

% operation of this function, consult the help page.

% quad() uses the following method to calculate the integral

% - The domain [a,b] is divided into ever-smaller sub-intervalues

% and the integrals over each sub-interval are added to return

% the value of the integral over the entire domain [a,b].

% - On each sub-interval, f(x) is approximated by using low-order

% polynomial interpolation (linear, quadratic, etc.), and for

% each interpolation polynomial, the value of the integral can

% be calculated analytically.

% - The estimated integral values for each order local polynomial

% approximation are then combined in such a way that the low-order

% (in Dx = the sub-interval size) contributions to the error

% cancel out (Romberg integration)

% - The routine iteratively refines the mesh until the estimated

% error is below some specified tolerance.

% As a rule, quad() is recommended for functions that are not very smooth,

% as only low-order local approximations are used on each sub-interval.

% First, we demonstrate the use of quad where we supply the function

% explicitly through a character string.

tol = 1e-7; % specify tolerance (default = 1e-6)

integral_quad_1 = quad(fun,a,b,tol);

disp(['Integral (quad w/ explicit function) = ', ...

num2str(integral_quad_1)]);

% We can also use quad with an inline function.

integral_quad_2 = quad(F,a,b,tol);

disp(['Integral (quad w/ inline function) = ', ...

num2str(integral_quad_2)]);

% For the development of programs, it is useful to be able to

% provide quad() with the name of a function that returns f(x)

% for an input x. As an example, we will use a routine that

% return the value of f(x) = exp(-c*x.^2), where for our

% particular problem we will use a value of c=1. The value of c

% for our problem is passed to the routine through the argument

% list of quad(). The name of the function to be evaluated is

% contained in the MATLAB file exp_sqr.m

% The contents of exp_sqr.m are

% % ======================================================

% % function f = exp_sqr(x,c);

% %

% % This routine returns, for input arguments of x (overloaded)

% % and c, the values of exp(-c*x.^2).

% function f = exp_sqr(x,c);

% f = exp(-c*x.^2);

% return;

fun_directory = pwd; % set directory name containing function

addpath(fun_directory); % add to path if not done so already

fun_name = 'exp_sqr'; % name of M-file containing function

c = 1; % value of constant c for our problem

trace = []; % if non-empty, quad() prints results of each iteration

integral_quad_3 = quad(fun_name,a,b,tol,trace,c);

disp(['Integral (quad w/ external routine) = ', ...

num2str(integral_quad_3)]);

% We can also call quad() for a routine containined in the same

% m-file as our calling program when we know the name of the

% function. As an example, we define exp_sqr_2() as a routine

% that performs exactly the same job as exp_sqr() but now add

% it to the calling M-file below.

integral_quad_4 = quad(@exp_sqr_2,a,b,tol,trace,c);

disp(['Integral (quad w/ named routine) = ', ...

num2str(integral_quad_4)]);

% quadl() ----------------------------------------------------------

% In addition to quad(), MATLAB also includes the built-in function

% quadl() that uses Gauss-Lobatto quadrature (a modification of

% Gaussian quadrature) to evaluate the integral over each sub-interval.

% In Gauss-Lobatto quadrature, one places grid points at the zeros

% of the othogonal polynomial (that each satisfy a < xj < b, AND

% at one or both end points (a and/or b).

% quadl() gives higher-accuracy when the integrand f(x) varies smoothly

% and so can be approximated well locally by a low-order polynomial.

% For functions that are not smooth, quad() is recommended, since the

% higher-order local polynomial approximations used in quadl() may oscillate

% over sub-intervals containing rapid changes in f(x).

% We demonstrate the use of quadl with the routine exp_sqr_2.

integral_quadl_1 = quadl(@exp_sqr_2,a,b,tol,trace,c);

disp(['Integral (quadl w/ named routine) = ', ...

num2str(integral_quadl_1)]);

% Now, let us say that we wish to integrate f(x), but only know

% the value of the function at a discrete set of grid points.

% We could use trapz(); however, this routine uses only a local

% linear interpolation, so that the method is not very accurate.

% We can use quadl() to obtain higher-order accuracy by using

% higher-order local approximations for f(x). In the example

% below, we use the MATLAB built-in function interp1() to

% return the estimated value of f(x) by using the tabulated values

% of f and x and piecewise cubic spline interpolation. Since

% this interpolation approximation for f(x) is everywhere a low-order

% polynomial locally, we expect quadl() to provide accurate results.

integral_spline = quadl(@cubic_interpolate,a,b,tol,trace, ...

x_val,f_val);

disp(['Integral (quadl w/ spline interpolation) = ', ...

num2str(integral_spline)]);

% ode15s() ---------------------------------------------------------

% Finally, we note that if we define the function y(x) as the integral

% from t = a to t = x of the function f(t), then we can obtain the

% value of the desired integral by solving the initial value problem

% y'(t) = f(t) with the initial contiion y(a) = 0.

% This method works if the function f(t) does not undergo singularities

% within [a,b], and can be quite efficient with an adaptive step-size

% algorithm for problems in which f(t) changes rapidly (but smoothly)

% over some region of [a,b].

% The built-in MATLAB function for solving initial value ODE problems,

% ode15s(), is used below to calculate the value of this integral.

t_span = [a b];

options = odeset('RelTol',1e-10,'AbsTol',1e-10);

[t_out,y_out] = ode15s(@exp_sqr_ode,t_span,0,options,c);

integral_ode15s = y_out(length(t_out));

disp(['Integral (ode15s) = ', ...

num2str(integral_ode15s)]);

% gaussq() ----------------------------------------------------------

% While quad() and quadl() are standard functions in MATLAB, one can

% write MATLAB functions that perform Guassian quadrature in a more

% explicit format. The MATLAB central file exchange at the URL :

% www.mathworks.com/matlabcentral/fileexchange/index.jsp

% includes a routine gaussq (written by Per Brodtkorb) that performs

% Gaussian quadrature using the most-common orthogonal polynomials.

% Finally, we calculate the errors for each method.

disp(' ');

disp('Error magnitudes for each method');

% trapz()

error_trap = abs(integral_trap - integral_true);

disp(['Error (trapz) = ', ...

num2str(error_trap)]);

% quad()

integral_quad_avg = ...

mean([integral_quad_1,integral_quad_2,integral_quad_3]);

error_quad = abs(integral_quad_avg - integral_true);

disp(['Error (quad) = ', ...

num2str(error_quad)]);

% quadl()

error_quadl = abs(integral_quadl_1 - integral_true);

disp(['Error (quadl) = ', ...

num2str(error_quadl)]);

% spline interpolation of grid values

error_spline = abs(integral_spline - integral_true);

disp(['Error (quadl w/ spline interpolation) = ', ...

num2str(error_spline)]);

% ode15s

error_ode15s = abs(integral_ode15s - integral_true);

disp(['Error (ode15s) = ', ...

num2str(error_ode15s)]);

iflag_main = 1;

return;

% =========================================================

% function f = exp_sqr_2(x,c);

% This routine returns, for input arguments of x (overloaded) and c,

% the values of exp(-c*x.^2).

function f = exp_sqr_2(x,c);

f = exp(-c*x.^2);

return;

% ===========================================================

% function f = cubic_interpolate(x,x_val,f_val);

% This MATLAB function uses piecewise cubic spline interpolation

% to return a values of f(x) from tabulated values of f(x) at the

% grid points x_val that are stored in f_val.

% K. Beers. MIT ChE. 8/7/2002

function f = cubic_interpolate(x,x_val,f_val);

f = interp1(x_val,f_val,x,'spline');

return;

% =========================================================

% function f = exp_sqr_ode(x,y,c);

% This routine returns, for input arguments of x (overloaded) and c,

% the values of exp(-c*x.^2). The value of y is the cumulative

% integral up to this point, and is not used to determine the value

% of f(x).

function f = exp_sqr_ode(x,y,c);

f = exp(-c*x.^2);

return;