(octave.info)Functions of Multiple Variables
23.3 Functions of Multiple Variables
====================================
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
-- : dblquad (F, XA, XB, YA, YB)
-- : dblquad (F, XA, XB, YA, YB, TOL)
-- : dblquad (F, XA, XB, YA, YB, TOL, QUADF)
-- : 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: Note: integral2, *note integral3:
XREFintegral3, Note: triplequad, *note quad:
XREFquad, Note: quadv, Note: quadl, Note:
quadgk, Note: quadcc, *note trapz:
XREFtrapz.
-- : triplequad (F, XA, XB, YA, YB, ZA, ZB)
-- : triplequad (F, XA, XB, YA, YB, ZA, ZB, TOL)
-- : triplequad (F, XA, XB, YA, YB, ZA, ZB, TOL, QUADF)
-- : 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: Note: integral3, *note integral2:
XREFintegral2, Note: dblquad, Note: quad,
Note: quadv, Note: quadl, *note quadgk:
XREFquadgk, Note: quadcc, Note: 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) where X is a vector and Y is a scalar. It should
return a vector of the same length and orientation as X.
Additional optional parameters can be specified using "PROPERTY",
VALUE pairs. Valid properties are:
‘AbsTol’
Define the absolute error tolerance for the quadrature. The
default value is 1e-10 (1e-5 for single).
‘RelTol’
Define the relative error tolerance for the quadrature. The
default value is 1e-6 (1e-4 for single).
‘MaxFunEvals’
The maximum number of function calls to the vectorized
function F. The default value is 5000.
‘Singular’
Enable/disable transforms to weaken singularities on the edge
of the integration domain. The default value is TRUE.
‘Vectorized’
Option to disable vectorized integration, forcing Octave to
use only scalar inputs when calling the integrand. The
default value is FALSE.
‘FailurePlot’
If ‘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, is
the Triangle Area x Height or ‘1/2 * BASE * WIDTH * HEIGHT’.
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: Note: integral2, *note dblquad:
XREFdblquad, Note: integral, Note: quad,
Note: quadgk, Note: quadv, *note quadl:
XREFquadl, Note: quadcc, Note: trapz, Note:
integral3, Note: 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) where X is a vector and Y is a scalar. It should
return a vector of the same length and orientation as X.
Additional optional parameters can be specified using "PROPERTY",
VALUE pairs. Valid properties are:
‘AbsTol’
Define the absolute error tolerance for the quadrature. The
default value is 1e-10 (1e-5 for single).
‘RelTol’
Define the relative error tolerance for the quadrature. The
default value is 1e-6 (1e-4 for single).
‘Method’
Specify 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.
‘Vectorized’
Enable 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 and only accept X and Y as scalars to be used. The
default value is ‘true’.
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, is
the Triangle Area x Height or ‘1/2 * BASE * WIDTH * HEIGHT’.
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: Note: quad2d, Note: dblquad,
Note: integral, Note: quad, *note quadgk:
XREFquadgk, Note: quadv, Note: quadl, Note:
quadcc, Note: trapz, *note integral3:
XREFintegral3, Note: 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) where X is a vector and Y is a scalar. It should
return a vector of the same length and orientation as X.
Additional optional parameters can be specified using "PROPERTY",
VALUE pairs. Valid properties are:
‘AbsTol’
Define the absolute error tolerance for the quadrature. The
default value is 1e-10 (1e-5 for single).
‘RelTol’
Define the relative error tolerance for the quadrature. The
default value is 1e-6 (1e-4 for single).
‘Method’
Specify 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.
‘Vectorized’
Enable 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 and only accept X and Y as scalars to be used. The
default value is ‘true’.
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
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) sqrt (1 - X.^2 - Y.^2);
Q = integral3 (F, 0, 1, 0, YMAX)
⇒ 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’.
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: Note: triplequad, *note integral:
XREFintegral, Note: quad, Note: quadgk, Note:
quadv, Note: quadl, Note: quadcc,
Note: trapz, Note: integral2, Note:
quad2d, Note: 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 (Note: 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.
automatically generated by info2www version 1.2.2.9