Info Catalog octave: Orthogonal Collocation octave: Numerical Integration

octave: Functions of Multiple Variables

 
 23.3 Functions of Multiple Variables
 ====================================
 
 Octave does not have built-in functions for computing the integral of
 functions of multiple variables directly.  It is possible, however, to
 compute the integral of a function of multiple variables using the
 existing functions for one-dimensional integrals.
 
    To illustrate how the integration can be performed, we will integrate
 the function
 
      f(x, y) = sin(pi*x*y)*sqrt(x*y)
 
    for x and y between 0 and 1.
 
    The first approach creates a function that integrates f with respect
 to x, and then integrates that function with respect to y.  Because
 ‘quad’ is written in Fortran it cannot be called recursively.  This
 means that ‘quad’ cannot integrate a function that calls ‘quad’, and
 hence cannot be used to perform the double integration.  Any of the
 other integrators, however, can be used which is what the following code
 demonstrates.
 
      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 above process can be simplified with the ‘dblquad’ and
 ‘triplequad’ functions for integrals over two and three variables.  For
 example:
 
      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: triplequad XREFtriplequad, quad XREFquad,
DONTPRINTYET       quadv XREFquadv, quadl XREFquadl, *notequadgk:
DONTPRINTYET       quadv XREFquadv, quadl XREFquadl, quadgk

      XREFquadgk, quadcc XREFquadcc, 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 ([]).
 
DONTPRINTYET       See also: dblquad XREFdblquad, quad XREFquad, *noteDONTPRINTYET       See also: dblquad XREFdblquad, quad XREFquad, 
      quadv XREFquadv, quadl XREFquadl, quadgk XREFquadgk,
      quadcc XREFquadcc, trapz XREFtrapz.
 
    The above mentioned approach works, but is fairly slow, and that
 problem increases exponentially with the dimensionality of the integral.
 Another possible solution is to use Orthogonal Collocation as described
 in the previous section (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.
 

Info Catalog octave: Orthogonal Collocation octave: Numerical Integration