Octave includes several functions for computing the integral of functions of multiple variables. This procedure can generally be performed by creating a function that integrates f with respect to x, and then integrates that function with respect to y. This procedure can be performed manually using the following example which integrates the function:
f(x, y) = sin(pi*x*y) * sqrt(x*y)
for x and y between 0 and 1.
Using quadgk in the example below, a double integration can be
performed. (Note that any of the 1-D quadrature functions can be used in this
fashion except for quad since it is written in Fortran and cannot be
called recursively.)
function q = g(y)
q = ones (size (y));
for i = 1:length (y)
f = @(x) sin (pi*x.*y(i)) .* sqrt (x.*y(i));
q(i) = quadgk (f, 0, 1);
endfor
endfunction
I = quadgk ("g", 0, 1)
⇒ 0.30022
The algorithm above is implemented in the function dblquad for integrals
over two variables. The 3-D equivalent of this process is implemented in
triplequad for integrals over three variables. As an example, the
result above can be replicated with a call to dblquad as shown below.
I = dblquad (@(x, y) sin (pi*x.*y) .* sqrt (x.*y), 0, 1, 0, 1)
⇒ 0.30022
q = dblquad (f, xa, xb, ya, yb) ¶q = dblquad (f, xa, xb, ya, yb, tol) ¶q = dblquad (f, xa, xb, ya, yb, tol, quadf) ¶q = dblquad (f, xa, xb, ya, yb, tol, quadf, …) ¶Numerically evaluate the double integral of f.
f is a function handle, inline function, or string containing the name of the function to evaluate. The function f must have the form z = f(x,y) where x is a vector and y is a scalar. It should return a vector of the same length and orientation as x.
xa, ya and xb, yb are the lower and upper limits of integration for x and y respectively. The underlying integrator determines whether infinite bounds are accepted.
The optional argument tol defines the absolute tolerance used to integrate each sub-integral. The default value is 1e-6.
The optional argument quadf specifies which underlying integrator
function to use. Any choice but quad is available and the default
is quadcc.
Additional arguments, are passed directly to f. To use the default
value for tol or quadf one may pass ':' or an empty
matrix ([]).
See also: integral2, integral3, triplequad, quad, quadv, quadl, quadgk, quadcc, trapz.
q = triplequad (f, xa, xb, ya, yb, za, zb) ¶q = triplequad (f, xa, xb, ya, yb, za, zb, tol) ¶q = triplequad (f, xa, xb, ya, yb, za, zb, tol, quadf) ¶q = triplequad (f, xa, xb, ya, yb, za, zb, tol, quadf, …) ¶Numerically evaluate the triple integral of f.
f is a function handle, inline function, or string containing the name of the function to evaluate. The function f must have the form w = f(x,y,z) where either x or y is a vector and the remaining inputs are scalars. It should return a vector of the same length and orientation as x or y.
xa, ya, za and xb, yb, zb are the lower and upper limits of integration for x, y, and z respectively. The underlying integrator determines whether infinite bounds are accepted.
The optional argument tol defines the absolute tolerance used to integrate each sub-integral. The default value is 1e-6.
The optional argument quadf specifies which underlying integrator
function to use. Any choice but quad is available and the default
is quadcc.
Additional arguments, are passed directly to f. To use the default
value for tol or quadf one may pass ':' or an empty
matrix ([]).
See also: integral3, integral2, dblquad, quad, quadv, quadl, quadgk, quadcc, trapz.
The recursive algorithm for quadrature presented above is referred to as
"iterated". A separate 2-D integration method is implemented in the
function quad2d. This function performs a "tiled" integration
by subdividing the integration domain into rectangular regions and performing
separate integrations over those domains. The domains are further subdivided
in areas requiring refinement to reach the desired numerical accuracy. For
certain functions this method can be faster than the 2-D iteration used in the
other functions above.
q = quad2d (f, xa, xb, ya, yb) ¶q = quad2d (f, xa, xb, ya, yb, prop, val, …) ¶[q, err, iter] = quad2d (…) ¶Numerically evaluate the two-dimensional integral of f using adaptive quadrature over the two-dimensional domain defined by xa, xb, ya, yb using tiled integration. Additionally, ya and yb may be scalar functions of x, allowing for the integration over non-rectangular domains.
f is a function handle, inline function, or string containing the
name of the function to evaluate. The function f must be of the form
z = f(x,y), and all operations must be vectorized such that x
and y accept array inputs and return array outputs of the same size.
(It can be assumed that x and y will either be same-size arrays
or one will be a scalar.) The underlying integrators will input arrays of
integration points into f and/or use internal vector expansions to
speed computation that can produce unpredictable results if f is not
restricted to elementwise operations. For integrands where this is
unavoidable, the ("Vectorized") option described below may produce
more reliable results.
Additional optional parameters can be specified using
"property", value pairs. Valid properties are:
AbsTolDefine the absolute error tolerance for the quadrature. The default value is 1e-10 (1e-5 for single).
RelTolDefine the relative error tolerance for the quadrature. The default value is 1e-6 (1e-4 for single).
MaxFunEvalsThe maximum number of function calls to the vectorized function f. The default value is 5000.
SingularEnable/disable transforms to weaken singularities on the edge of the integration domain. The default value is true.
VectorizedEnable or disable vectorized integration. A value of false forces
Octave to use only scalar inputs when calling the integrand, which enables
integrands f(x,y) that have not been vectorized or only accept
scalar values of x or y. The default value is true.
Note that this is achieved by wrapping f(x,y) with the function
arrayfun, which may significantly decrease computation speed.
FailurePlotIf quad2d fails to converge to the desired error tolerance before
MaxFunEvals is reached, a plot of the areas that still need refinement
is created. The default value is false.
Adaptive quadrature is used to minimize the estimate of error until the following is satisfied:
error <= max (AbsTol, RelTol*|q|)
The optional output err is an approximate bound on the error in the
integral abs (q - I), where I is the exact value
of the integral. The optional output iter is the number of vectorized
function calls to the function f that were used.
Example 1 : integrate a rectangular region in x-y plane
f = @(x,y) 2*ones (size (x)); q = quad2d (f, 0, 1, 0, 1) ⇒ q = 2
The result is a volume, which for this constant-value integrand, is just
Length * Width * Height.
Example 2 : integrate a triangular region in x-y plane
f = @(x,y) 2*ones (size (x)); ymax = @(x) 1 - x; q = quad2d (f, 0, 1, 0, ymax) ⇒ q = 1
The result is a volume, which for this constant-value integrand
f = 2, is the Triangle Area x Height or
1/2 * Base * Width * Height.
Example 3 : integrate a non-vectorized function over a square region
f = @(x,y) sinc (x) * sinc (y)); q = quad2d (f, -1, 1, -1, 1) ⇒ q = 12.328 (incorrect) q = quad2d (f, -1, 1, -1, 1, "Vectorized", false) ⇒ q = 1.390 (correct) f = @(x,y) sinc (x) .* sinc (y); q = quad2d (f, -1, 1, -1, 1) ⇒ q = 1.390 (correct)
The first result is incorrect as the non-elementwise operator between the
sinc functions in f create unintended matrix multiplications between
the internal integration arrays used by quad2d. In the second
result, setting "Vectorized" to false forces quad2d to
perform scalar internal operations to compute the integral, resulting in
the correct numerical result at the cost of about a 20x increase in
computation time. In the third result, vectorizing the integrand f
using the elementwise multiplication operator gets the correct result
without increasing computation time.
Programming Notes: If there are singularities within the integration region it is best to split the integral and place the singularities on the boundary.
Known MATLAB incompatibility: If tolerances are left unspecified, and
any integration limits are of type single, then Octave’s integral
functions automatically reduce the default absolute and relative error
tolerances as specified above. If tighter tolerances are desired they
must be specified. MATLAB leaves the tighter tolerances appropriate
for double inputs in place regardless of the class of the
integration limits.
Reference: L.F. Shampine, MATLAB program for quadrature in 2D, Applied Mathematics and Computation, pp. 266–274, Vol 1, 2008.
See also: integral2, dblquad, integral, quad, quadgk, quadv, quadl, quadcc, trapz, integral3, triplequad.
Finally, the functions integral2 and integral3 are provided
as general 2-D and 3-D integration functions. They will auto-select between
iterated and tiled integration methods and, unlike dblquad and
triplequad, will work with non-rectangular integration domains.
q = integral2 (f, xa, xb, ya, yb) ¶q = integral2 (f, xa, xb, ya, yb, prop, val, …) ¶[q, err] = integral2 (…) ¶Numerically evaluate the two-dimensional integral of f using adaptive quadrature over the two-dimensional domain defined by xa, xb, ya, yb (scalars may be finite or infinite). Additionally, ya and yb may be scalar functions of x, allowing for integration over non-rectangular domains.
f is a function handle, inline function, or string containing the
name of the function to evaluate. The function f must be of the form
z = f(x,y), and all operations must be vectorized such that x
and y accept array inputs and return array outputs of the same size.
(It can be assumed that x and y will either be same-size arrays
or one will be a scalar.) The underlying integrators will input arrays of
integration points into f and/or use internal vector expansions to
speed computation that can produce unpredictable results if f is not
restricted to elementwise operations. For integrands where this is
unavoidable, the ("Vectorized") option described below may produce
more reliable results.
Additional optional parameters can be specified using
"property", value pairs. Valid properties are:
AbsTolDefine the absolute error tolerance for the quadrature. The default value is 1e-10 (1e-5 for single).
RelTolDefine the relative error tolerance for the quadrature. The default value is 1e-6 (1e-4 for single).
MethodSpecify the two-dimensional integration method to be used, with valid
options being "auto" (default), "tiled", or
"iterated". When using "auto", Octave will choose the
"tiled" method unless any of the integration limits are infinite.
VectorizedEnable or disable vectorized integration. A value of false forces
Octave to use only scalar inputs when calling the integrand, which enables
integrands f(x,y) that have not been vectorized or only accept
scalar values of x or y. The default value is true.
Note that this is achieved by wrapping f(x,y) with the function
arrayfun, which may significantly decrease computation speed.
Adaptive quadrature is used to minimize the estimate of error until the following is satisfied:
error <= max (AbsTol, RelTol*|q|)
err is an approximate bound on the error in the integral
abs (q - I), where I is the exact value of the
integral.
Example 1 : integrate a rectangular region in x-y plane
f = @(x,y) 2*ones (size (x)); q = integral2 (f, 0, 1, 0, 1) ⇒ q = 2
The result is a volume, which for this constant-value integrand, is just
Length * Width * Height.
Example 2 : integrate a triangular region in x-y plane
f = @(x,y) 2*ones (size (x)); ymax = @(x) 1 - x; q = integral2 (f, 0, 1, 0, ymax) ⇒ q = 1
The result is a volume, which for this constant-value integrand
f = 2, is the Triangle Area x Height or
1/2 * Base * Width * Height.
Example 3 : integrate a non-vectorized function over a square region
f = @(x,y) sinc (x) * sinc (y)); q = integral2 (f, -1, 1, -1, 1) ⇒ q = 12.328 (incorrect) q = integral2 (f, -1, 1, -1, 1, "Vectorized", false) ⇒ q = 1.390 (correct) f = @(x,y) sinc (x) .* sinc (y); q = integral2 (f, -1, 1, -1, 1) ⇒ q = 1.390 (correct)
The first result is incorrect as the non-elementwise operator between the
sinc functions in f create unintended matrix multiplications between
the internal integration arrays used by integral2. In the second
result, setting "Vectorized" to false forces integral2 to
perform scalar internal operations to compute the integral, resulting in
the correct numerical result at the cost of about a 20x increase in
computation time. In the third result, vectorizing the integrand f
using the elementwise multiplication operator gets the correct result
without increasing computation time.
Programming Notes: If there are singularities within the integration region it is best to split the integral and place the singularities on the boundary.
Known MATLAB incompatibility: If tolerances are left unspecified, and
any integration limits are of type single, then Octave’s integral
functions automatically reduce the default absolute and relative error
tolerances as specified above. If tighter tolerances are desired they
must be specified. MATLAB leaves the tighter tolerances appropriate
for double inputs in place regardless of the class of the
integration limits.
Reference: L.F. Shampine, MATLAB program for quadrature in 2D, Applied Mathematics and Computation, pp. 266–274, Vol 1, 2008.
See also: quad2d, dblquad, integral, quad, quadgk, quadv, quadl, quadcc, trapz, integral3, triplequad.
q = integral3 (f, xa, xb, ya, yb, za, zb) ¶q = integral3 (f, xa, xb, ya, yb, za, zb, prop, val, …) ¶Numerically evaluate the three-dimensional integral of f using adaptive quadrature over the three-dimensional domain defined by xa, xb, ya, yb, za, zb (scalars may be finite or infinite). Additionally, ya and yb may be scalar functions of x and za, and zb maybe be scalar functions of x and y, allowing for integration over non-rectangular domains.
f is a function handle, inline function, or string containing the
name of the function to evaluate. The function f must be of the form
z = f(x,y,z), and all operations must be vectorized such that
x, y, and z accept array inputs and return array outputs
of the same size. (It can be assumed that x, y, and z
will either be same-size arrays or scalars.) The underlying integrators
will input arrays of integration points into f and/or use internal
vector expansions to speed computation that can produce unpredictable
results if f is not restricted to elementwise operations. For
integrands where this is unavoidable, the ("Vectorized") option
described below may produce more reliable results.
Additional optional parameters can be specified using
"property", value pairs. Valid properties are:
AbsTolDefine the absolute error tolerance for the quadrature. The default value is 1e-10 (1e-5 for single).
RelTolDefine the relative error tolerance for the quadrature. The default value is 1e-6 (1e-4 for single).
MethodSpecify the two-dimensional integration method to be used, with valid
options being "auto" (default), "tiled", or
"iterated". When using "auto", Octave will choose the
"tiled" method unless any of the integration limits are infinite.
VectorizedEnable or disable vectorized integration. A value of false forces
Octave to use only scalar inputs when calling the integrand, which enables
integrands f(x,y,z) that have not been vectorized or only accept
scalar values of x, y, or z. The default value is
true. Note that this is achieved by wrapping f(x,y,z) with
the function arrayfun, which may significantly decrease computation
speed.
Adaptive quadrature is used to minimize the estimate of error until the following is satisfied:
error <= max (AbsTol, RelTol*|q|)
err is an approximate bound on the error in the integral
abs (q - I), where I is the exact value of the
integral.
Example 1 : integrate over a rectangular volume
f = @(x,y,z) ones (size (x)); q = integral3 (f, 0, 1, 0, 1, 0, 1) ⇒ q = 1.00000
For this constant-value integrand, the result is a volume which is just
Length * Width * Height.
Example 2 : integrate over a spherical volume
f = @(x,y) ones (size (x)); ymax = @(x) sqrt (1 - x.^2); zmax = @(x,y) sqrt (1 - x.^2 - y.^2); q = integral3 (f, 0, 1, 0, ymax, 0, zmax) ⇒ q = 0.52360
For this constant-value integrand, the result is a volume which is 1/8th
of a unit sphere or 1/8 * 4/3 * pi.
Example 3 : integrate a non-vectorized function over a cubic volume
f = @(x,y) sinc (x) * sinc (y), * sinc (z); q = integral3 (f, -1, 1, -1, 1, -1, 1) ⇒ q = 14.535 (incorrect) q = integral3 (f, -1, 1, -1, 1, -1, 1, "Vectorized", false) ⇒ q = 1.6388 (correct) f = @(x,y,z) sinc (x) .* sinc (y), .* sinc (z); q = integral3 (f, -1, 1, -1, 1, -1, 1) ⇒ q = 1.6388 (correct)
The first result is incorrect as the non-elementwise operator between the
sinc functions in f create unintended matrix multiplications between
the internal integration arrays used by integral3. In the second
result, setting "Vectorized" to false forces integral3 to
perform scalar internal operations to compute the integral, resulting in
the correct numerical result at the cost of about a 30x increase in
computation time. In the third result, vectorizing the integrand f
using the elementwise multiplication operator gets the correct result
without increasing computation time.
Programming Notes: If there are singularities within the integration region it is best to split the integral and place the singularities on the boundary.
Known MATLAB incompatibility: If tolerances are left unspecified, and
any integration limits are of type single, then Octave’s integral
functions automatically reduce the default absolute and relative error
tolerances as specified above. If tighter tolerances are desired they
must be specified. MATLAB leaves the tighter tolerances appropriate
for double inputs in place regardless of the class of the
integration limits.
Reference: L.F. Shampine, MATLAB program for quadrature in 2D, Applied Mathematics and Computation, pp. 266–274, Vol 1, 2008.
See also: triplequad, integral, quad, quadgk, quadv, quadl, quadcc, trapz, integral2, quad2d, dblquad.
The above integrations can be fairly slow, and that problem increases
exponentially with the dimensionality of the integral. Another possible
solution for 2-D integration is to use Orthogonal Collocation as described in
the previous section (see Orthogonal Collocation). The integral of a
function f(x,y) for x and y between 0 and 1 can be
approximated using n points by
the sum over i=1:n and j=1:n of q(i)*q(j)*f(r(i),r(j)),
where q and r is as returned by colloc (n). The
generalization to more than two variables is straight forward. The
following code computes the studied integral using n=8 points.
f = @(x,y) sin (pi*x*y') .* sqrt (x*y');
n = 8;
[t, ~, ~, q] = colloc (n);
I = q'*f(t,t)*q;
⇒ 0.30022
It should be noted that the number of points determines the quality of the approximation. If the integration needs to be performed between a and b, instead of 0 and 1, then a change of variables is needed.