@c Copyright (C) 1996, 1997 John W. Eaton
@c This is part of the Octave manual.
@c For copying conditions, see the file gpl.texi.

@node Numerical Integration
@chapter Numerical Integration

Octave comes with several built-in functions for computing the integral
of a function numerically.  These functions all solve 1-dimensional
integration problems.

@menu
* Functions of One Variable:: 
* Functions of Multiple Variables:: 
* Orthogonal Collocation::      
@end menu

@node Functions of One Variable
@section Functions of One Variable

Octave supports three different algorithms for computing the integral
@iftex
@tex
$$
 \int_a^b f(x) d x
$$
@end tex
@end iftex
of a function @math{f} over the interval from @math{a} to @math{b}.
These are

@table @code
@item quad
Numerical integration based on Gaussian quadrature.

@item quadl
Numerical integration using an adaptive Lobatto rule.

@item trapz
Numerical integration using the trapezodial method.
@end table

@noindent
Besides these functions Octave also allows you to perform cumulative
numerical integration using the trapezodial method through the
@code{cumtrapz} function.

@DOCSTRING(quad)

@DOCSTRING(quad_options)

Here is an example of using @code{quad} to integrate the function
@iftex
@tex
$$
 f(x) = x \sin (1/x) \sqrt {|1 - x|}
$$
from $x = 0$ to $x = 3$.
@end tex
@end iftex
@ifnottex

@example
  @var{f}(@var{x}) = @var{x} * sin (1/@var{x}) * sqrt (abs (1 - @var{x}))
@end example

@noindent
from @var{x} = 0 to @var{x} = 3.
@end ifnottex

This is a fairly difficult integration (plot the function over the range
of integration to see why).

The first step is to define the function:

@example
@group
function y = f (x)
  y = x .* sin (1 ./ x) .* sqrt (abs (1 - x));
endfunction
@end group
@end example

Note the use of the `dot' forms of the operators.  This is not necessary
for the call to @code{quad}, but it makes it much easier to generate a
set of points for plotting (because it makes it possible to call the
function with a vector argument to produce a vector result).

Then we simply call quad:

@example
@group
[v, ier, nfun, err] = quad ("f", 0, 3)
     @result{} 1.9819
     @result{} 1
     @result{} 5061
     @result{} 1.1522e-07
@end group
@end example

Although @code{quad} returns a nonzero value for @var{ier}, the result
is reasonably accurate (to see why, examine what happens to the result
if you move the lower bound to 0.1, then 0.01, then 0.001, etc.).

@DOCSTRING(quadl)

@DOCSTRING(trapz)

@DOCSTRING(cumtrapz)

@node Orthogonal Collocation
@section Orthogonal Collocation

@DOCSTRING(colloc)

Here is an example of using @code{colloc} to generate weight matrices
for solving the second order differential equation
@iftex
@tex
$u^\prime - \alpha u^{\prime\prime} = 0$ with the boundary conditions
$u(0) = 0$ and $u(1) = 1$.
@end tex
@end iftex
@ifnottex
@var{u}' - @var{alpha} * @var{u}'' = 0 with the boundary conditions
@var{u}(0) = 0 and @var{u}(1) = 1.
@end ifnottex

First, we can generate the weight matrices for @var{n} points (including
the endpoints of the interval), and incorporate the boundary conditions
in the right hand side (for a specific value of
@iftex
@tex
$\alpha$).
@end tex
@end iftex
@ifnottex
@var{alpha}).
@end ifnottex

@example
@group
n = 7;
alpha = 0.1;
[r, a, b] = colloc (n-2, "left", "right");
at = a(2:n-1,2:n-1);
bt = b(2:n-1,2:n-1);
rhs = alpha * b(2:n-1,n) - a(2:n-1,n);
@end group
@end example

Then the solution at the roots @var{r} is

@example
u = [ 0; (at - alpha * bt) \ rhs; 1]
     @result{} [ 0.00; 0.004; 0.01 0.00; 0.12; 0.62; 1.00 ]
@end example

@node Functions of Multiple Variables
@section Functions of Multiple Variables

Octave does not have built-in functions for computing the integral
of functions of multiple variables.  It is however possible to compute
the integral of a function of multiple variables using the functions
for one-dimensional integrals.

To illustrate how the integration can be performed, we will integrate
the function
@iftex
@tex
$$
  f(x, y) = \sin(\pi x y)\sqrt{x y}
$$
@end tex
@end iftex
@ifnottex
@example
f(x, y) = sin(pi*x*y)*sqrt(x*y)
@end example
@end ifnottex
for @math{x} and @math{y} between 0 and 1.

The first approach creates a function that integrates @math{f} with
respect to @math{x}, and then integrates that function with respect to
@math{y}.  Since @code{quad} is written in Fortran it cannot be called
recursively.  This means that @code{quad} cannot integrate a function
that calls @code{quad}, and hence cannot be used to perform the double
integration.  It is however possible with @code{quadl}, which is what
the following code does.

@example
function I = g(y)
  I = ones(1, length(y));
  for i = 1:length(y)
    f = @@(x) sin(pi.*x.*y(i)).*sqrt(x.*y(i));
    I(i) = quadl(f, 0, 1);
  endfor
endfunction

I = quadl("g", 0, 1)
      @result{} 0.30022
@end example

The above mentioned approach works but is fairly slow, and that problem
increases exponentially with the dimensionality the problem.  Another
possible solution is to use Orthogonal Collocation as described in the
previous section.  The integral of a function @math{f(x,y)} for
@math{x} and @math{y} between 0 and 1 can be approximated using @math{n}
points by
@iftex
@tex
$$
 \int_0^1 \int_0^1 f(x,y) d x d y \approx \sum_{i=1}^n \sum_{j=1}^n q_i q_j f(r_i, r_j),
$$
@end tex
@end iftex
@ifnottex
the sum over @code{i=1:n} and @code{j=1:n} of @code{q(i)*q(j)*f(r(i),r(j))},
@end ifnottex
where @math{q} and @math{r} is as returned by @code{colloc(n)}.  The
generalisation to more than two variables is straight forward.  The
following code computes the studied integral using @math{n=7} points.

@example
f = @@(x,y) sin(pi*x*y').*sqrt(x*y');
n = 7;
[t, A, B, q] = colloc(n);
I = q'*f(t,t)*q;
      @result{} 0.30022
@end example

@noindent
It should be noted that the number of points determines the quality
of the approximation.  If the integration needs to be performed between
@math{a} and @math{b} instead of 0 and 1, a change of variables is needed.