@menu
* Introduction to Integration::  
* Definitions for Integration::  
* Introduction to QUADPACK::
* Definitions for QUADPACK::
@end menu

@node Introduction to Integration, Definitions for Integration, Integration, Integration
@section Introduction to Integration

Maxima has several routines for handling integration.
The @code{integrate} function makes use of most of them.  There is also the
@code{antid} package, which handles an unspecified function (and its
derivatives, of course).  For numerical uses, there is the @code{romberg}
function; an
adaptave integrator which uses the Newton-Cotes 8 panel quadrature
rule, called @code{quanc8}; and a set of adaptive integrators from Quadpack,
named @code{quad_qag}, @code{quad_qags}, etc., which are described under the heading @code{QUADPACK}.
Hypergeometric functions are being worked on,
see @code{specint} for details.
Generally speaking, Maxima only handles integrals which are
integrable in terms of the "elementary functions" (rational functions,
trigonometrics, logs, exponentials, radicals, etc.) and a few
extensions (error function, dilogarithm).  It does not handle
integrals in terms of unknown functions such as @code{g(x)} and @code{h(x)}.

@c end concepts Integration
@node Definitions for Integration, Introduction to QUADPACK, Introduction to Integration, Integration
@section Definitions for Integration

@c NEEDS WORK
@deffn {Function} changevar (@var{expr}, @var{f(x,y)}, @var{y}, @var{x})
Makes the change of variable given by
@code{@var{f(x,y)} = 0} in all integrals occurring in @var{expr} with integration with
respect to @var{x}.
The new variable is @var{y}.

@c HMM, THIS EXAMPLE YIELDS A CORRECT BUT SLIGHTLY STRANGE RESULT...
@example
(%i1) assume(a > 0)$
(%i2) 'integrate (%e**sqrt(a*y), y, 0, 4);
                      4
                     /
                     [    sqrt(a) sqrt(y)
(%o2)                I  %e                dy
                     ]
                     /
                      0
(%i3) changevar (%, y-z^2/a, z, y);
                      0
                     /
                     [                abs(z)
                   2 I            z %e       dz
                     ]
                     /
                      - 2 sqrt(a)
(%o3)            - ----------------------------
                                a
@end example

An expression containing a noun form, such as the instances of @code{'integrate} above,
may be evaluated by @code{ev} with the @code{nouns} flag.
For example, the expression returned by @code{changevar} above may be evaluated
by @code{ev (%o3, nouns)}.

@code{changevar} may also be used to changes in the indices of a sum or
product.  However, it must be realized that when a change is made in a
sum or product, this change must be a shift, i.e., @code{i = j+ ...}, not a
higher degree function.  E.g.,

@example
(%i4) sum (a[i]*x^(i-2), i, 0, inf);
                         inf
                         ====
                         \         i - 2
(%o4)                     >    a  x
                         /      i
                         ====
                         i = 0
(%i5) changevar (%, i-2-n, n, i);
                        inf
                        ====
                        \               n
(%o5)                    >      a      x
                        /        n + 2
                        ====
                        n = - 2
@end example

@end deffn

@c THIS ITEM IS A MESS, BUT DON'T BOTHER TO CLEAN IT UP:
@c THE GAUSS-KRONROD FUNCTIONS (QUADPACK) MAKE THIS OBSOLETE
@deffn {Function} dblint (@var{f}, @var{r}, @var{s}, @var{a}, @var{b})
A double-integral routine which was written in
top-level Maxima and then translated and compiled to machine code.
Use @code{load (dblint)} to access this package.  It uses the Simpson's rule
method in both the x and y directions to calculate

@example
/b /s(x)
|  |
|  |    f(x,y) dy dx
|  |
/a /r(x)
@end example

The function @var{f} must be a translated or compiled function of two
variables, and @var{r} and @var{s} must each be a translated or compiled
function of one variable, while @var{a} and @var{b} must be floating point
numbers. The routine has two global variables which determine the
number of divisions of the x and y intervals: @code{dblint_x} and @code{dblint_y},
both of which are initially 10, and can be changed independently to
other integer values (there are @code{2*dblint_x+1} points computed in the x
direction, and @code{2*dblint_y+1} in the y direction).
The routine subdivides the X axis and then for each value of X it
first computes @code{@var{r}(x)} and @code{@var{s}(x)}; then the Y axis between @code{@var{r}(x)} and @code{@var{s}(x)} is
subdivided and the integral along the Y axis is performed using
Simpson's rule; then the integral along the X axis is done using
Simpson's rule with the function values being the Y-integrals. This
procedure may be numerically unstable for a great variety of reasons,
but is reasonably fast: avoid using it on highly oscillatory functions
and functions with singularities (poles or branch points in the
region).  The Y integrals depend on how far apart @code{@var{r}(x)} and @code{@var{s}(x)} are,
so if the distance @code{@var{s}(x) - @var{r}(x)} varies rapidly with X, there may be
substantial errors arising from truncation with different step-sizes
in the various Y integrals. One can increase @code{dblint_x} and @code{dblint_y} in
an effort to improve the coverage of the region, at the expense of
computation time. The function values are not saved, so if the
function is very time-consuming, you will have to wait for
re-computation if you change anything (sorry).
It is required that the functions @var{f}, @var{r}, and @var{s} be either translated or
compiled prior to calling @code{dblint}. This will result in orders of
magnitude speed improvement over interpreted code in many cases!

@code{demo (dblint)} executes a demonstration of @code{dblint} applied to an example problem.
@c demo (dblint_1) FAILS WITH Could not find `fltdfnk.mc' -- DON'T BOTHER TO MENTION IT. !!!
@c @code{demo (dblint_1)} executes another demonstration.

@end deffn

@deffn {Function} defint (@var{expr}, @var{x}, @var{a}, @var{b})
Attempts to compute a definite integral.
@code{defint} is called by @code{integrate} when limits of integration are specified,
i.e., when @code{integrate} is called as @code{integrate (@var{expr}, @var{x}, @var{a}, @var{b})}.
Thus from the user's point of view, it is sufficient to call @code{integrate}.
@c SHOULD WE BOTHER TO DOCUMENT defint ??? NO FUNCTIONALITY HERE THAT IS NOT ALREADY PRESENT IN integrate !!!

@code{defint} returns a symbolic expression,
either the computed integral or the noun form of the integral.
See @code{quad_qag} and related functions for numerical approximation of definite integrals.

@end deffn

@c NEEDS EXPANSION AND EXAMPLES
@deffn {Function} erf (@var{x})
Represents the error function, whose derivative is:
@code{2*exp(-x^2)/sqrt(%pi)}.

@end deffn

@defvr {Option variable} erfflag
Default value: @code{true}

When @code{erfflag} is @code{false}, prevents @code{risch} from introducing the
@code{erf} function in the answer if there were none in the integrand to
begin with.

@end defvr

@c NEEDS WORK
@deffn {Function} ilt (@var{expr}, @var{t}, @var{s})
Computes the inverse Laplace transform of @var{expr} with
respect to @var{t} and parameter @var{s}.  @var{expr} must be a ratio of
polynomials whose denominator has only linear and quadratic factors.
By using the functions @code{laplace} and @code{ilt} together with the @code{solve} or
@code{linsolve} functions the user can solve a single differential or
convolution integral equation or a set of them.

@example
(%i1) 'integrate (sinh(a*x)*f(t-x), x, 0, t) + b*f(t) = t**2;
              t
             /
             [                                    2
(%o1)        I  f(t - x) sinh(a x) dx + b f(t) = t
             ]
             /
              0
(%i2) laplace (%, t, s);
                               a laplace(f(t), t, s)   2
(%o2)  b laplace(f(t), t, s) + --------------------- = --
                                       2    2           3
                                      s  - a           s
(%i3) linsolve ([%], ['laplace(f(t), t, s)]);
                                        2      2
                                     2 s  - 2 a
(%o3)     [laplace(f(t), t, s) = --------------------]
                                    5         2     3
                                 b s  + (a - a  b) s
(%i4) ilt (rhs (first (%)), s, t);
Is  a b (a b - 1)  positive, negative, or zero?

pos;
               sqrt(a b (a b - 1)) t
        2 cosh(---------------------)       2
                         b               a t
(%o4) - ----------------------------- + -------
              3  2      2               a b - 1
             a  b  - 2 a  b + a

                                                       2
                                             + ------------------
                                                3  2      2
                                               a  b  - 2 a  b + a
@end example

@end deffn

@deffn {Function} integrate (@var{expr}, @var{x})
@deffnx {Function} integrate (@var{expr}, @var{x}, @var{a}, @var{b})
Attempts to symbolically compute the integral of @var{expr} with respect to @var{x}.
@code{integrate (@var{expr}, @var{x})} is an indefinite integral,
while @code{integrate (@var{expr}, @var{x}, @var{a}, @var{b})} is a definite integral,
with limits of integration @var{a} and @var{b}.
The limits should not contain @var{x}, although @code{integrate} does not enforce this restriction.
@var{a} need not be less than @var{b}.
If @var{b} is equal to @var{a}, @code{integrate} returns zero.

See @code{quad_qag} and related functions for numerical approximation of definite integrals.
See @code{residue} for computation of residues (complex integration).
See @code{antid} for an alternative means of computing indefinite integrals.

The integral (an expression free of @code{integrate}) is returned if @code{integrate} succeeds.
Otherwise the return value is
the noun form of the integral (the quoted operator @code{'integrate})
or an expression containing one or more noun forms.
The noun form of @code{integrate} is displayed with an integral sign.

In some circumstances it is useful to construct a noun form by hand,
by quoting @code{integrate} with a single quote, e.g., @code{'integrate (@var{expr}, @var{x})}.
For example, the integral may depend on some parameters which are not yet computed.
The noun may be applied to its arguments by @code{ev (@var{i}, nouns)}
where @var{i} is the noun form of interest.

@c BEGIN EXPOSITION ON HEURISTICS
@code{integrate} handles definite integrals separately from indefinite,
and employs a range of heuristics to handle each case.
Special cases of definite integrals include limits of integration equal to
zero or infinity (@code{inf} or @code{minf}),
trigonometric functions with limits of integration equal to zero and @code{%pi} or @code{2 %pi},
rational functions,
integrals related to the definitions of the @code{beta} and @code{psi} functions,
and some logarithmic and trigonometric integrals.
Processing rational functions may include computation of residues.
If an applicable special case is not found,
an attempt will be made to compute the indefinite integral and evaluate it at the limits of integration.
This may include taking a limit as a limit of integration goes to infinity or negative infinity;
see also @code{ldefint}.

Special cases of indefinite integrals include trigonometric functions,
exponential and logarithmic functions,
and rational functions.
@code{integrate} may also make use of a short table of elementary integrals.

@code{integrate} may carry out a change of variable
if the integrand has the form @code{f(g(x)) * diff(g(x), x)}.
@code{integrate} attempts to find a subexpression @code{g(x)} such that
the derivative of @code{g(x)} divides the integrand.
This search may make use of derivatives defined by the @code{gradef} function.
See also @code{changevar} and @code{antid}.

If none of the preceding heuristics find the indefinite integral,
the Risch algorithm is executed.
The flag @code{risch} may be set as an @code{evflag},
in a call to @code{ev} or on the command line,
e.g., @code{ev (integrate (@var{expr}, @var{x}), risch)} or @code{integrate (@var{expr}, @var{x}), risch}.
If @code{risch} is present, @code{integrate} calls the @code{risch} function
without attempting heuristics first. See also @code{risch}.
@c END EXPOSITION ON HEURISTICS

@code{integrate} works only with functional relations represented explicitly with the @code{f(x)} notation.
@code{integrate} does not respect implicit dependencies established by the @code{depends} function.

@code{integrate} may need to know some property of a parameter in the integrand.
@code{integrate} will first consult the @code{assume} database,
and, if the variable of interest is not there,
@code{integrate} will ask the user.
Depending on the question,
suitable responses are @code{yes;} or @code{no;},
or @code{pos;}, @code{zero;}, or @code{neg;}.

@code{integrate} is not, by default, declared to be linear. See @code{declare} and @code{linear}.

@code{integrate} attempts integration by parts only in a few special cases.

Examples:

@itemize @bullet
@item
Elementary indefinite and definite integrals.

@example
(%i1) integrate (sin(x)^3, x);
                           3
                        cos (x)
(%o1)                   ------- - cos(x)
                           3
(%i2) integrate (x/ sqrt (b^2 - x^2), x);
                                 2    2
(%o2)                    - sqrt(b  - x )
(%i3) integrate (cos(x)^2 * exp(x), x, 0, %pi);
                               %pi
                           3 %e      3
(%o3)                      ------- - -
                              5      5
(%i4) integrate (x^2 * exp(-x^2), x, minf, inf);
                            sqrt(%pi)
(%o4)                       ---------
                                2
@end example

@item
Use of @code{assume} and interactive query.

@example
(%i1) assume (a > 1)$
(%i2) integrate (x**a/(x+1)**(5/2), x, 0, inf);
    2 a + 2
Is  -------  an integer?
       5

no;
Is  2 a - 3  positive, negative, or zero?

neg;
                                   3
(%o2)                  beta(a + 1, - - a)
                                   2
@end example

@item
Change of variable. There are two changes of variable in this example:
one using a derivative established by @code{gradef},
and one using the derivation @code{diff(r(x))} of an unspecified function @code{r(x)}.

@example
(%i3) gradef (q(x), sin(x**2));
(%o3)                         q(x)
(%i4) diff (log (q (r (x))), x);
                      d               2
                     (-- (r(x))) sin(r (x))
                      dx
(%o4)                ----------------------
                            q(r(x))
(%i5) integrate (%, x);
(%o5)                     log(q(r(x)))
@end example

@item
Return value contains the @code{'integrate} noun form.
In this example, Maxima can extract one factor of the denominator
of a rational function, but cannot factor the remainder or otherwise find its integral.
@code{grind} shows the noun form @code{'integrate} in the result.
See also @code{integrate_use_rootsof} for more on integrals of rational functions.

@example
(%i1) expand ((x-4) * (x^3+2*x+1));
                    4      3      2
(%o1)              x  - 4 x  + 2 x  - 7 x - 4
(%i2) integrate (1/%, x);
                              /  2
                              [ x  + 4 x + 18
                              I ------------- dx
                              ]  3
                 log(x - 4)   / x  + 2 x + 1
(%o2)            ---------- - ------------------
                     73               73
(%i3) grind (%);
log(x-4)/73-('integrate((x^2+4*x+18)/(x^3+2*x+1),x))/73$
@end example

@item
Defining a function in terms of an integral.
The body of a function is not evaluated when the function is defined.
Thus the body of @code{f_1} in this example contains the noun form of @code{integrate}.
The quote-quote operator @code{'@w{}'} causes the integral to be evaluated,
and the result becomes the body of @code{f_2}.

@example
(%i1) f_1 (a) := integrate (x^3, x, 1, a);
                                     3
(%o1)           f_1(a) := integrate(x , x, 1, a)
(%i2) ev (f_1 (7), nouns);
(%o2)                          600
(%i3) /* Note parentheses around integrate(...) here */
      f_2 (a) := ''(integrate (x^3, x, 1, a));
                                   4
                                  a    1
(%o3)                   f_2(a) := -- - -
                                  4    4
(%i4) f_2 (7);
(%o4)                          600
@end example
@end itemize

@end deffn

@defvr {System variable} integration_constant_counter
Default value: 0

@c WHEN DOES integrationconstant1 SHOW UP IN THE OUTPUT OF integrate ???
@c integrate (a, x) YIELDS "a x", NOT "a x + integrationconstant1" !!!
@code{integration_constant_counter} is a counter which is updated each time a
constant of integration (named by Maxima, e.g., @code{integrationconstant1})
is introduced into an expression by indefinite integration of an equation.

@end defvr

@defvr {Option variable} integrate_use_rootsof
Default value: @code{false}

When @code{integrate_use_rootsof} is @code{true} and the denominator of
a rational function cannot be factored, @code{integrate} returns the integral
in a form which is a sum over the roots (not yet known) of the denominator.

For example, with @code{integrate_use_rootsof} set to @code{false},
@code{integrate} returns an unsolved integral of a rational function in noun form:

@example
(%i1) integrate_use_rootsof: false$
(%i2) integrate (1/(1+x+x^5), x);
        /  2
        [ x  - 4 x + 5
        I ------------ dx                            2 x + 1
        ]  3    2                2            5 atan(-------)
        / x  - x  + 1       log(x  + x + 1)          sqrt(3)
(%o2)   ----------------- - --------------- + ---------------
                7                 14             7 sqrt(3)
@end example

Now we set the flag to be true and the unsolved part of the
integral will be expressed as a summation over the roots of the denominator of the rational function:

@example
(%i3) integrate_use_rootsof: true$
(%i4) integrate (1/(1+x+x^5), x);
      ====        2
      \       (%r4  - 4 %r4 + 5) log(x - %r4)
       >      -------------------------------
      /                    2
      ====            3 %r4  - 2 %r4
                      3    2
      %r4 in rootsof(x  - x  + 1)
(%o4) ----------------------------------------------------------
               7

                                                             2 x + 1
                                         2            5 atan(-------)
                                    log(x  + x + 1)          sqrt(3)
                                  - --------------- + ---------------
                                          14             7 sqrt(3)
@end example

Alternatively the user may compute the roots of the denominator separately,
and then express the integrand in terms of these roots,
e.g., @code{1/((x - a)*(x - b)*(x - c))} or @code{1/((x^2 - (a+b)*x + a*b)*(x - c))}
if the denominator is a cubic polynomial.
Sometimes this will help Maxima obtain a more useful result.

@end defvr

@c NEEDS EXAMPLES
@deffn {Function} ldefint (@var{expr}, @var{x}, @var{a}, @var{b})
Attempts to compute the definite integral of @var{expr} by using
@code{limit} to evaluate the indefinite integral of @var{expr} with respect to @var{x}
at the upper limit @var{b} and at the lower limit @var{a}.
If it fails to compute the definite integral,
@code{ldefint} returns an expression containing limits as noun forms.

@code{ldefint} is not called from @code{integrate},
so executing @code{ldefint (@var{expr}, @var{x}, @var{a}, @var{b})} may yield a different result than
@code{integrate (@var{expr}, @var{x}, @var{a}, @var{b})}.
@code{ldefint} always uses the same method to evaluate the definite integral,
while @code{integrate} may employ various heuristics and may recognize some special cases.

@end deffn

@c UMM, IS THERE SOME TEXT MISSING HERE ???
@c WHAT IS THIS ABOUT EXACTLY ??
@deffn {Function} potential (@var{givengradient})
The calculation makes use of the global variable @code{potentialzeroloc[0]}
which must be @code{nonlist} or of the form

@example
[indeterminatej=expressionj, indeterminatek=expressionk, ...]
@end example

the
former being equivalent to the nonlist expression for all right-hand
sides in the latter.  The indicated right-hand sides are used as the
lower limit of integration.  The success of the integrations may
depend upon their values and order. @code{potentialzeroloc} is initially set
to 0.

@end deffn

@c THIS ITEM IS A MESS BUT DON'T BOTHER TO FIX IT:
@c THE GAUSS-KRONROD FUNCTIONS (QUADPACK) MAKE THIS OBSOLETE
@deffn {Function} qq
The package @code{qq} (which may be loaded with @code{load ("qq")})
contains a function @code{quanc8} which can take either 3 or 4 arguments. The
3 arg version computes the integral of the function specified as the
first argument over the interval from lo to hi as in
@code{quanc8 ('function, lo, hi)}.
The function name should be quoted.  The 4 arg version will compute
the integral of the function or expression (first arg) with respect to
the variable (second arg) over the interval from @code{lo} to @code{hi} as in
@code{quanc8(<f(x) or expression in x>, x, lo, hi)}.
The method used is the Newton-Cotes 8th order polynomial quadrature,
and the routine is adaptive. It will thus spend time dividing the
interval only when necessary to achieve the error conditions specified
by the global variables @code{quanc8_relerr} (default value=1.0e-4) and
@code{quanc8_abserr} (default value=1.0e-8) which give the relative error
test:

@example
|integral(function) - computed value| < quanc8_relerr*|integral(function)|
@end example

and the absolute error test:

@example
|integral(function) - computed value| < quanc8_abserr
@end example

@code{printfile ("qq.usg")} yields additional information.

@end deffn

@deffn {Function} quanc8 (@var{expr}, @var{a}, @var{b})
An adaptive integrator.
Demonstration and usage files are provided.  The method is to
use Newton-Cotes 8-panel quadrature rule, hence the function name
@code{quanc8}, available in 3 or 4 arg versions.  Absolute and relative error
checks are used.  To use it do @code{load ("qq")}.  See also @code{qq}.

@end deffn

@deffn {Function} residue (@var{expr}, @var{z}, @var{z_0})
Computes the residue in the complex plane of
the expression @var{expr} when the variable @var{z} assumes the value @var{z_0}.  The
residue is the coefficient of @code{(@var{z} - @var{z_0})^(-1)} in the Laurent series
for @var{expr}.

@example
(%i1) residue (s/(s**2+a**2), s, a*%i);
                                1
(%o1)                           -
                                2
(%i2) residue (sin(a*x)/x**4, x, 0);
                                 3
                                a
(%o2)                         - --
                                6
@end example

@end deffn

@deffn {Function} risch (@var{expr}, @var{x})
Integrates @var{expr} with respect to @var{x} using the
transcendental case of the Risch algorithm.  (The algebraic case of
the Risch algorithm has not been implemented.)  This currently
handles the cases of nested exponentials and logarithms which the main
part of @code{integrate} can't do.  @code{integrate} will automatically apply @code{risch}
if given these cases.

@code{erfflag}, if @code{false}, prevents @code{risch} from introducing the @code{erf}
function in the answer if there were none in the integrand to begin
with.

@example
(%i1) risch (x^2*erf(x), x);
                                                        2
             3                      2                - x
        %pi x  erf(x) + (sqrt(%pi) x  + sqrt(%pi)) %e
(%o1)   -------------------------------------------------
                              3 %pi
(%i2) diff(%, x), ratsimp;
                             2
(%o2)                       x  erf(x)
@end example

@end deffn

@c NEEDS WORK BUT DON'T BOTHER TO CLEAN IT UP:
@c THE GAUSS-KRONROD FUNCTIONS (QUADPACK) MAKE THIS OBSOLETE
@deffn {Function} romberg (@var{expr}, @var{x}, @var{a}, @var{b})
@deffnx {Function} romberg (@var{expr}, @var{a}, @var{b})
Romberg integration.
There are two ways to use this function.  The first is an inefficient
way like the definite integral version of @code{integrate}:
@code{romberg (<integrand>, <variable of integration>, <lower limit>, <upper limit>)}.

Examples:

@example
(%i1) showtime: true$
(%i2) romberg (sin(y), y, 0, %pi);
Evaluation took 0.00 seconds (0.01 elapsed) using 25.293 KB.
(%o2)                   2.000000016288042
(%i3) 1/((x-1)^2+1/100) + 1/((x-2)^2+1/1000) + 1/((x-3)^2+1/200)$
(%i4) f(x) := ''%$
(%i5) rombergtol: 1e-6$
(%i6) rombergit: 15$
(%i7) romberg (f(x), x, -5, 5);
Evaluation took 11.97 seconds (12.21 elapsed) using 12.423 MB.
(%o7)                   173.6730736617464
@c INCLUDE THIS COMPARISON TO EXACT RESULT ??? YIELDS A LOT OF "RAT replaced" MESSAGES !!!
@c integrate (f(x), x, -5, 5) - %, numer;
@end example

The second is an efficient way that is used as follows:

@example
romberg (<function name>, <lower limit>, <upper limit>);
@end example

Continuing the above example, we have:
@example
(%i8) f(x) := (mode_declare ([function(f), x], float), ''(%th(5)))$
(%i9) translate(f);
(%o9)                          [f]
(%i10) romberg (f, -5, 5);
Evaluation took 3.51 seconds (3.86 elapsed) using 6.641 MB.
(%o10)                  173.6730736617464
@end example

The first argument must be a translated or compiled function.  (If it
is compiled it must be declared to return a @code{flonum}.)  If the first
argument is not already translated, @code{romberg} will not attempt to
translate it but will give an error.

The accuracy of the integration is governed by the global variables
@code{rombergtol} (default value 1.E-4) and @code{rombergit} (default value 11).
@code{romberg} will return a result if the relative difference in successive
approximations is less than @code{rombergtol}.  It will try halving the
stepsize @code{rombergit} times before it gives up.  The number of iterations
and function evaluations which @code{romberg} will do is governed by
@code{rombergabs} and @code{rombergmin}.

@code{romberg} may be called recursively and thus can do double and triple
integrals.

Example:
@example
(%i1) assume (x > 0)$
(%i2) integrate (integrate (x*y/(x+y), y, 0, x/2), x, 1, 3)$
(%i3) radcan (%);
                    26 log(3) - 26 log(2) - 13
(%o3)             - --------------------------
                                3
(%i4) %,numer;
(%o4)                   .8193023963959073
(%i5) define_variable (x, 0.0, float, "Global variable in function F")$
(%i6) f(y) := (mode_declare (y, float), x*y/(x+y))$
(%i7) g(x) := romberg ('f, 0, x/2)$  
(%i8) romberg (g, 1, 3);
(%o8)                   .8193022864324522
@end example

The advantage with this way is that the function @code{f} can be used for other 
purposes, like plotting. The disadvantage is that you have to think up 
a name for both the function @code{f} and its free variable @code{x}.
Or, without the global:

@example
(%i1) g_1(x) := (mode_declare (x, float), romberg (x*y/(x+y), y, 0, x/2))$
(%i2) romberg (g_1, 1, 3);
(%o2)                   .8193022864324522
@end example

The advantage here is shortness.

@example
(%i3) q (a, b) := romberg (romberg (x*y/(x+y), y, 0, x/2), x, a, b)$
(%i4) q (1, 3);
(%o4)                   .8193022864324522
@end example

It is even shorter this way, and the variables do not need to be declared 
because they are in the context of @code{romberg}.
Use of @code{romberg} for multiple integrals can have great disadvantages,
though.  The amount of extra calculation needed because of the
geometric information thrown away by expressing multiple integrals
this way can be incredible.  The user should be sure to understand and
use the @code{rombergtol} and @code{rombergit} switches.

@end deffn

@defvr {Option variable} rombergabs
Default value: 0.0

Assuming that successive estimates
produced by @code{romberg} are @code{y[0]}, @code{y[1]}, @code{y[2]}, etc., then @code{romberg} will
return after @code{n} iterations if (roughly speaking)

@example
 (abs(y[n]-y[n-1]) <= rombergabs or
 abs(y[n]-y[n-1])/(if y[n]=0.0 then 1.0 else y[n]) <= rombergtol)
@end example

is @code{true}.  (The condition on the number of iterations given by
@code{rombergmin} must also be satisfied.)
Thus if @code{rombergabs} is 0.0 (the default) you just get the relative
error test.  The usefulness of the additional variable comes when you
want to perform an integral, where the dominant contribution comes
from a small region.  Then you can do the integral over the small
dominant region first, using the relative accuracy check, followed by
the integral over the rest of the region using the absolute accuracy
check.

Example:  Suppose you want to compute

@example
'integrate (exp(-x), x, 0, 50)
@end example

(numerically) with a relative accuracy of  1 part in 10000000.
Define the function.  @code{n} is a counter, so we can see how many
function evaluations were needed.
First of all try doing the whole integral at once.

@example
(%i1) f(x) := (mode_declare (n, integer, x, float), n:n+1, exp(-x))$
(%i2) translate(f)$
Warning-> n is an undefined global variable.
(%i3) block ([rombergtol: 1.e-6, romberabs: 0.0], n:0, romberg (f, 0, 50));
(%o3)                   1.000000000488271
(%i4) n;
(%o4)                          257
@end example

That approach required 257 function evaluations.
 Now do the integral intelligently, by first doing
     @code{'integrate (exp(-x), x, 0, 10)} and then setting @code{rombergabs} to 1.E-6 times (this
     partial integral).
This approach takes only 130 function evaluations.

@example
(%i5) block ([rombergtol: 1.e-6, rombergabs:0.0, sum:0.0],
  n: 0, sum: romberg (f, 0, 10), rombergabs: sum*rombergtol, rombergtol:0.0,
      sum + romberg (f, 10, 50));
(%o5)                   1.000000001234793
(%i6) n;
(%o6)                          130
@end example

So if @code{f(x)} were a function that took a long time to compute, the
second method would be about 2 times quicker.

@end defvr

@defvr {Option variable} rombergit
Default value: 11

The accuracy of the @code{romberg} integration
command is governed by the global variables @code{rombergtol} and
@code{rombergit}.  @code{romberg} will return a result if the relative
difference in successive approximations is less than @code{rombergtol}.  It
will try halving the stepsize @code{rombergit} times before it gives up.

@end defvr

@defvr {Option variable} rombergmin
Default value: 0

@code{rombergmin} governs the minimum number of function
evaluations that @code{romberg} will make.  @code{romberg} will evaluate its first
arg. at least @code{2^(rombergmin+2)+1} times.  This is useful for
integrating oscillatory functions, when the normal converge test might
sometimes wrongly pass.

@end defvr

@defvr {Option variable} rombergtol
Default value: 1e-4

The accuracy of the @code{romberg} integration
command is governed by the global variables @code{rombergtol} and
@code{rombergit}.  @code{romberg} will return a result if the relative
difference in successive approximations is less than @code{rombergtol}.  It
will try halving the stepsize @code{rombergit} times before it gives up.

@end defvr

@c NEEDS EXPANSION, CLARIFICATION, AND EXAMPLES
@deffn {Function} tldefint (@var{expr}, @var{x}, @var{a}, @var{b})
Equivalent to @code{ldefint} with @code{tlimswitch} set to @code{true}.

@end deffn

@footnotestyle end

@node Introduction to QUADPACK, Definitions for QUADPACK, Definitions for Integration, Integration
@section Introduction to QUADPACK


@c FOLLOWING TEXT ADAPTED WITH HEAVY MODIFICATION FROM http://www.netlib.org/slatec/src/qpdoc.f

QUADPACK is a collection of functions for the numerical
computation of one-dimensional definite integrals.
It originated from a joint project of
R. Piessens @footnote{Applied Mathematics and Programming Division, K.U. Leuven},
E. de Doncker @footnote{Applied Mathematics and Programming Division, K.U. Leuven},
C. Ueberhuber @footnote{Institut f@"ur Mathematik, T.U. Wien},
and D. Kahaner @footnote{National Bureau of Standards, Washington, D.C., U.S.A}.

The QUADPACK library included in Maxima is an automatic translation
(via the program @code{f2cl}) of the Fortran source code of QUADPACK as it appears in
the SLATEC Common Mathematical Library, Version 4.1 @footnote{http://www.netlib.org/slatec}.
The SLATEC library is dated July 1993, but the QUADPACK functions
were written some years before.
There is another version of QUADPACK at Netlib @footnote{http://www.netlib.org/quadpack};
it is not clear how that version differs from the SLATEC version.

The QUADPACK functions included in Maxima are all automatic,
in the sense that these functions attempt to compute a result to a specified accuracy,
requiring an unspecified number of function evaluations.
Maxima's Lisp translation of QUADPACK also includes some non-automatic functions,
but they are not exposed at the Maxima level.

Further information about QUADPACK can be found in the QUADPACK book
@footnote{R. Piessens, E. de Doncker-Kapenga, C.W. Uberhuber, and D.K. Kahaner.
@i{QUADPACK: A Subroutine Package for Automatic Integration.}
Berlin: Springer-Verlag, 1983, ISBN 0387125531.}.

@subsection Overview

@table @code
@item quad_qag
Integration of a general function over a finite interval.
@code{quad_qag} implements a simple globally adaptive integrator using the strategy of Aind (Piessens, 1973).
The caller may choose among 6 pairs of Gauss-Kronrod quadrature
formulae for the rule evaluation component.
The high-degree rules are suitable for strongly oscillating integrands.

@item quad_qags
Integration of a general function over a finite interval.
@code{quad_qags} implements globally adaptive interval subdivision with extrapolation
(de Doncker, 1978) by the Epsilon algorithm (Wynn, 1956).

@item quad_qagi
Integration of a general function over an infinite or semi-infinite interval.
The interval is mapped onto a finite interval and
then the same strategy as in @code{quad_qags} is applied.

@item quad_qawo
Integration of @math{cos(omega x) f(x)} or @math{sin(omega x) f(x)} over a finite interval,
where @math{omega} is a constant.
The rule evaluation component is based on the modified Clenshaw-Curtis technique.
@code{quad_qawo} applies adaptive subdivision with extrapolation, similar to @code{quad_qags}.

@item quad_qawf
Calculates a Fourier cosine or Fourier sine transform on a semi-infinite interval.
The same approach as in @code{quad_qawo} is applied on successive finite intervals,
and convergence acceleration by means of the Epsilon algorithm (Wynn, 1956) 
is applied to the series of the integral contributions.

@item quad_qaws
Integration of @math{w(x) f(x)} over a finite interval @math{[a, b]},
where @math{w} is a function of the form @math{(x - a)^alpha (b - x)^beta v(x)}
and @math{v(x)} is 1 or @math{log(x - a)} or @math{log(b - x)} or @math{log(x - a) log(b - x)},
and @math{alpha > -1} and @math{beta > -1}.
A globally adaptive subdivision strategy is applied,
with modified Clenshaw-Curtis integration on the subintervals which contain @math{a} or @math{b}.

@item quad_qawc
Computes the Cauchy principal value of @math{f(x)/(x - c)} over a finite interval @math{(a, b)}
and specified @math{c}.
The strategy is globally adaptive, and modified
Clenshaw-Curtis integration is used on the subranges
which contain the point @math{x = c}.
@end table

@node Definitions for QUADPACK, , Introduction to QUADPACK, Integration
@section Definitions for QUADPACK

@c THERE ARE OPTIONAL ARGUMENTS WHICH MAKES LISTING THE VARIANTS A LITTLE TEDIOUS
@c NEED A MORE CONVENIENT (AND NONAMBIGUOUS) NOTATION FOR OPTIONAL ARGUMENTS
@deffn {Function} quad_qag (@var{f(x)}, @var{x}, @var{a}, @var{b}, @var{key}, @var{epsrel}, @var{limit})
@deffnx {Function} quad_qag (@var{f}, @var{x}, @var{a}, @var{b}, @var{key}, @var{epsrel}, @var{limit})

Integration of a general function over a finite interval.
@code{quad_qag} implements a simple globally adaptive integrator using the strategy of Aind (Piessens, 1973).
The caller may choose among 6 pairs of Gauss-Kronrod quadrature
formulae for the rule evaluation component.
The high-degree rules are suitable for strongly oscillating integrands.

@code{quad_qag} computes the integral

@ifhtml
@math{integrate (f(x), x, a, b)}
@end ifhtml
@ifinfo
@math{integrate (f(x), x, a, b)}
@end ifinfo
@tex
$$\int_a^b {f(x) dx}$$
@end tex

The function to be integrated is @var{f(x)}, with dependent
variable @var{x}, and the function is to be integrated between the
limits @var{a} and @var{b}.  @var{key} is the integrator to be used
and should be an integer between 1 and 6, inclusive.  The value of
@var{key} selects the order of the Gauss-Kronrod integration rule.
High-order rules are suitable for strongly oscillating integrands.

The integrand may be specified as the name of a Maxima or Lisp function or operator,
a Maxima lambda expression, or a general Maxima expression.

The numerical integration is done adaptively by subdividing the
integration region into sub-intervals until the desired accuracy is
achieved.

The optional arguments @var{epsrel} and @var{limit} are the desired
relative error and the maximum number of subintervals, respectively.
@var{epsrel} defaults to 1e-8 and @var{limit} is 200.

@code{quad_qag} returns a list of four elements:

@itemize
@item
an approximation to the integral,
@item
the estimated absolute error of the approximation, 
@item
the number integrand evaluations,
@item
an error code.
@end itemize

The error code (fourth element of the return value) can have the values:

@table @code
@item 0
if no problems were encountered; 
@item 1
if too many sub-intervals were done;
@item 2
if excessive roundoff error is detected;
@item 3
if extremely bad integrand behavior occurs;
@item 6
if the input is invalid.

@end table

@c NEED CROSS REFS HERE -- EITHER CROSS REF A QUADPACK OVERVIEW, OR CROSS REF EACH OF THE quad_* FUNCTIONS

Examples:

@example
(%i1) quad_qag (x^(1/2)*log(1/x), x, 0, 1, 3);
(%o1)    [.4444444444492108, 3.1700968502883E-9, 961, 0]
(%i2) integrate (x^(1/2)*log(1/x), x, 0, 1);
                                4
(%o2)                           -
                                9
@end example
@end deffn

@c THERE ARE OPTIONAL ARGUMENTS WHICH MAKES LISTING THE VARIANTS A LITTLE TEDIOUS
@c NEED A MORE CONVENIENT (AND NONAMBIGUOUS) NOTATION FOR OPTIONAL ARGUMENTS
@deffn {Function} quad_qags (@var{f(x)}, @var{x}, @var{a}, @var{b}, @var{epsrel}, @var{limit})
@deffnx {Function} quad_qags (@var{f}, @var{x}, @var{a}, @var{b}, @var{epsrel}, @var{limit})

Integration of a general function over a finite interval.
@code{quad_qags} implements globally adaptive interval subdivision with extrapolation
(de Doncker, 1978) by the Epsilon algorithm (Wynn, 1956).

@code{quad_qags} computes the integral

@ifhtml
@math{integrate (f(x), x, a, b)}
@end ifhtml
@ifinfo
@math{integrate (f(x), x, a, b)}
@end ifinfo
@tex
$$\int_a^b {f(x) dx}$$
@end tex

The function to be integrated is @var{f(x)}, with
dependent variable @var{x}, and the function is to be integrated
between the limits @var{a} and @var{b}.

The integrand may be specified as the name of a Maxima or Lisp function or operator,
a Maxima lambda expression, or a general Maxima expression.

The optional arguments @var{epsrel} and @var{limit} are the desired
relative error and the maximum number of subintervals, respectively.
@var{epsrel} defaults to 1e-8 and @var{limit} is 200.

@code{quad_qags} returns a list of four elements:

@itemize
@item
an approximation to the integral,
@item
the estimated absolute error of the approximation, 
@item
the number integrand evaluations,
@item
an error code.
@end itemize

The error code (fourth element of the return value) can have the values:

@table @code
@item 0
no problems were encountered; 
@item 1
too many sub-intervals were done;
@item 2
excessive roundoff error is detected;
@item 3
extremely bad integrand behavior occurs;
@item 4
failed to converge
@item 5
integral is probably divergent or slowly convergent
@item 6
if the input is invalid.
@end table

@c NEED CROSS REFS HERE -- EITHER CROSS REF A QUADPACK OVERVIEW, OR CROSS REF EACH OF THE quad_* FUNCTIONS

Examples:

@example
(%i1) quad_qags (x^(1/2)*log(1/x), x, 0 ,1);
(%o1)   [.4444444444444448, 1.11022302462516E-15, 315, 0]
@end example

Note that @code{quad_qags} is more accurate and efficient than @code{quad_qag} for this integrand.

@end deffn

@c THERE ARE OPTIONAL ARGUMENTS WHICH MAKES LISTING THE VARIANTS A LITTLE TEDIOUS
@c NEED A MORE CONVENIENT (AND NONAMBIGUOUS) NOTATION FOR OPTIONAL ARGUMENTS
@deffn {Function} quad_qagi (@var{f(x)}, @var{x}, @var{a}, @var{inftype}, @var{epsrel}, @var{limit})
@deffnx {Function} quad_qagi (@var{f}, @var{x}, @var{a}, @var{inftype}, @var{epsrel}, @var{limit})

Integration of a general function over an infinite or semi-infinite interval.
The interval is mapped onto a finite interval and
then the same strategy as in @code{quad_qags} is applied.

@code{quad_qagi} evaluates one of the following integrals

@ifhtml
@math{integrate (f(x), x, a, inf)}
@end ifhtml
@ifinfo
@math{integrate (f(x), x, a, inf)}
@end ifinfo
@tex
$$\int_a^\infty {f(x) dx}$$
@end tex

@ifhtml
@math{integrate (f(x), x, minf, a)}
@end ifhtml
@ifinfo
@math{integrate (f(x), x, minf, a)}
@end ifinfo
@tex
$$\int_\infty^a {f(x) dx}$$
@end tex

@ifhtml
@math{integrate (f(x), x, minf, inf)}
@end ifhtml
@ifinfo
@math{integrate (f(x), x, minf, inf)}
@end ifinfo
@tex
$$\int_{-\infty}^\infty {f(x) dx}$$
@end tex

using the Quadpack QAGI routine.  The function to be integrated is
@var{f(x)}, with dependent variable @var{x}, and the function is to
be integrated over an infinite range.

The integrand may be specified as the name of a Maxima or Lisp function or operator,
a Maxima lambda expression, or a general Maxima expression.

The parameter @var{inftype} determines the integration interval as follows:

@table @code
@item inf
The interval is from @var{a} to positive infinity.
@item minf
The interval is from negative infinity to @var{a}.
@item both
The interval is the entire real line.
@end table

The optional arguments @var{epsrel} and @var{limit} are the desired
relative error and the maximum number of subintervals, respectively.
@var{epsrel} defaults to 1e-8 and @var{limit} is 200.

@code{quad_qagi} returns a list of four elements:

@itemize
@item
an approximation to the integral,
@item
the estimated absolute error of the approximation, 
@item
the number integrand evaluations,
@item
an error code.
@end itemize

The error code (fourth element of the return value) can have the values:

@table @code
@item 0
no problems were encountered; 
@item 1
too many sub-intervals were done;
@item 2
excessive roundoff error is detected;
@item 3
extremely bad integrand behavior occurs;
@item 4
failed to converge
@item 5
integral is probably divergent or slowly convergent
@item 6
if the input is invalid.

@end table

@c NEED CROSS REFS HERE -- EITHER CROSS REF A QUADPACK OVERVIEW, OR CROSS REF EACH OF THE quad_* FUNCTIONS

Examples:

@example
(%i1) quad_qagi (x^2*exp(-4*x), x, 0, inf);
(%o1)        [0.03125, 2.95916102995002E-11, 105, 0]
(%i2) integrate (x^2*exp(-4*x), x, 0, inf);
                               1
(%o2)                          --
                               32
@end example

@end deffn

@c THERE ARE OPTIONAL ARGUMENTS WHICH MAKES LISTING THE VARIANTS A LITTLE TEDIOUS
@c NEED A MORE CONVENIENT (AND NONAMBIGUOUS) NOTATION FOR OPTIONAL ARGUMENTS
@deffn {Function} quad_qawc (@var{f(x)}, @var{x}, @var{c}, @var{a}, @var{b}, @var{epsrel}, @var{limit})
@deffnx {Function} quad_qawc (@var{f}, @var{x}, @var{c}, @var{a}, @var{b}, @var{epsrel}, @var{limit})

Computes the Cauchy principal value of @math{f(x)/(x - c)} over a finite interval.
The strategy is globally adaptive, and modified
Clenshaw-Curtis integration is used on the subranges
which contain the point @math{x = c}.

@code{quad_qawc} computes the Cauchy principal value of

@ifhtml
@math{integrate (f(x)/(x - c), x, a, b)}
@end ifhtml
@ifinfo
@math{integrate (f(x)/(x - c), x, a, b)}
@end ifinfo
@tex
$$\int_{a}^{b}{{{f\left(x\right)}\over{x-c}}\>dx}$$
@end tex

using the Quadpack QAWC routine.  The function to be integrated is
@code{@var{f(x)}/(@var{x} - @var{c})}, with dependent variable @var{x}, and the function
is to be integrated over the interval @var{a} to @var{b}.

The integrand may be specified as the name of a Maxima or Lisp function or operator,
a Maxima lambda expression, or a general Maxima expression.

The optional arguments @var{epsrel} and @var{limit} are the desired
relative error and the maximum number of subintervals, respectively.
@var{epsrel} defaults to 1e-8 and @var{limit} is 200.

@code{quad_qawc} returns a list of four elements:

@itemize
@item
an approximation to the integral,
@item
the estimated absolute error of the approximation, 
@item
the number integrand evaluations,
@item
an error code.
@end itemize

The error code (fourth element of the return value) can have the values:

@table @code
@item 0
no problems were encountered; 
@item 1
too many sub-intervals were done;
@item 2
excessive roundoff error is detected;
@item 3
extremely bad integrand behavior occurs;
@item 6
if the input is invalid.

@end table

Examples:

@example
(%i1) quad_qawc (2^(-5)*((x-1)^2+4^(-5))^(-1), x, 2, 0, 5);
(%o1)    [- 3.130120337415925, 1.306830140249558E-8, 495, 0]
(%i2) integrate (2^(-alpha)*(((x-1)^2 + 4^(-alpha))*(x-2))^(-1), x, 0, 5);
Principal Value
                       alpha
        alpha       9 4                 9
       4      log(------------- + -------------)
                      alpha           alpha
                  64 4      + 4   64 4      + 4
(%o2) (-----------------------------------------
                        alpha
                     2 4      + 2

         3 alpha                       3 alpha
         -------                       -------
            2            alpha/2          2          alpha/2
      2 4        atan(4 4       )   2 4        atan(4       )   alpha
    - --------------------------- - -------------------------)/2
                alpha                        alpha
             2 4      + 2                 2 4      + 2
(%i3) ev (%, alpha=5, numer);
(%o3)                    - 3.130120337415917
@end example

@end deffn

@c THERE ARE OPTIONAL ARGUMENTS WHICH MAKES LISTING THE VARIANTS A LITTLE TEDIOUS
@c NEED A MORE CONVENIENT (AND NONAMBIGUOUS) NOTATION FOR OPTIONAL ARGUMENTS
@deffn {Function} quad_qawf (@var{f(x)}, @var{x}, @var{a}, @var{omega}, @var{trig}, @var{epsabs}, @var{limit}, @var{maxp1}, @var{limlst})
@deffnx {Function} quad_qawf (@var{f}, @var{x}, @var{a}, @var{omega}, @var{trig}, @var{epsabs}, @var{limit}, @var{maxp1}, @var{limlst})

Calculates a Fourier cosine or Fourier sine transform on a semi-infinite interval
using the Quadpack QAWF function.
The same approach as in @code{quad_qawo} is applied on successive finite intervals,
and convergence acceleration by means of the Epsilon algorithm (Wynn, 1956) 
is applied to the series of the integral contributions.

@code{quad_qawf} computes the integral

@ifhtml
@math{integrate (f(x)*w(x), x, a, inf)}
@end ifhtml
@ifinfo
@math{integrate (f(x)*w(x), x, a, inf)}
@end ifinfo
@tex
$$\int_a^\infty f(x) w(x) dx$$
@end tex

The weight function @math{w} is selected by @var{trig}:

@table @code
@item cos
@math{w(x) = cos (omega x)}
@item sin
@math{w(x) = sin (omega x)}
@end table

The integrand may be specified as the name of a Maxima or Lisp function or operator,
a Maxima lambda expression, or a general Maxima expression.

The optional arguments are:

@table @var
@item epsabs
Desired absolute error of approximation.  Default is 1d-10.
@item limit
Size of internal work array.  (@var{limit} - @var{limlst})/2 is the
maximum number of subintervals to use.  Default is 200.
@item maxp1
Maximum number of Chebyshev moments.  Must be greater than 0.  Default
is 100.
@item limlst
Upper bound on the number of cycles.  Must be greater than or equal to
3.  Default is 10.
@end table

@c MERGE THESE LINES INTO PRECEDING TABLE
@var{epsabs} and @var{limit} are the desired
relative error and the maximum number of subintervals, respectively.
@var{epsrel} defaults to 1e-8 and @var{limit} is 200.

@code{quad_qawf} returns a list of four elements:

@itemize
@item
an approximation to the integral,
@item
the estimated absolute error of the approximation, 
@item
the number integrand evaluations,
@item
an error code.
@end itemize

The error code (fourth element of the return value) can have the values:

@table @code
@item 0
no problems were encountered; 
@item 1
too many sub-intervals were done;
@item 2
excessive roundoff error is detected;
@item 3
extremely bad integrand behavior occurs;
@item 6
if the input is invalid.

@end table

Examples:

@example
(%i1) quad_qawf (exp(-x^2), x, 0, 1, 'cos);
(%o1)   [.6901942235215714, 2.84846300257552E-11, 215, 0]
(%i2) integrate (exp(-x^2)*cos(x), x, 0, inf);
                          - 1/4
                        %e      sqrt(%pi)
(%o2)                   -----------------
                                2
(%i3) ev (%, numer);
(%o3)                   .6901942235215714
@end example

@end deffn

@c THERE ARE OPTIONAL ARGUMENTS WHICH MAKES LISTING THE VARIANTS A LITTLE TEDIOUS
@c NEED A MORE CONVENIENT (AND NONAMBIGUOUS) NOTATION FOR OPTIONAL ARGUMENTS
@deffn {Function} quad_qawo (@var{f(x)}, @var{x}, @var{a}, @var{b}, @var{omega}, @var{trig}, @var{epsabs}, @var{limit}, @var{maxp1}, @var{limlst})
@deffnx {Function} quad_qawo (@var{f}, @var{x}, @var{a}, @var{b}, @var{omega}, @var{trig}, @var{epsabs}, @var{limit}, @var{maxp1}, @var{limlst})

Integration of @math{cos(omega x) f(x)} or @math{sin(omega x) f(x)} over a finite interval,
where @math{omega} is a constant.
The rule evaluation component is based on the modified Clenshaw-Curtis technique.
@code{quad_qawo} applies adaptive subdivision with extrapolation, similar to @code{quad_qags}.

@code{quad_qawo} computes the integral using the Quadpack QAWO
routine:

@ifhtml
@math{integrate (f(x)*w(x), x, a, b)}
@end ifhtml
@ifinfo
@math{integrate (f(x)*w(x), x, a, b)}
@end ifinfo
@tex
$$\int_a^b f(x) w(x) dx$$
@end tex

The weight function @math{w} is selected by @var{trig}:

@table @code
@item cos
@math{w(x) = cos (omega x)}
@item sin
@math{w(x) = sin (omega x)}
@end table

The integrand may be specified as the name of a Maxima or Lisp function or operator,
a Maxima lambda expression, or a general Maxima expression.

The optional arguments are:

@table @var
@item epsabs
Desired absolute error of approximation.  Default is 1d-10.
@item limit
Size of internal work array.  (@var{limit} - @var{limlst})/2 is the
maximum number of subintervals to use.  Default is 200.
@item maxp1
Maximum number of Chebyshev moments.  Must be greater than 0.  Default
is 100.
@item limlst
Upper bound on the number of cycles.  Must be greater than or equal to
3.  Default is 10.
@end table

@c MERGE THESE LINES INTO PRECEDING TABLE
@var{epsabs} and @var{limit} are the desired
relative error and the maximum number of subintervals, respectively.
@var{epsrel} defaults to 1e-8 and @var{limit} is 200.

@code{quad_qawo} returns a list of four elements:

@itemize
@item
an approximation to the integral,
@item
the estimated absolute error of the approximation, 
@item
the number integrand evaluations,
@item
an error code.
@end itemize

The error code (fourth element of the return value) can have the values:

@table @code
@item 0
no problems were encountered; 
@item 1
too many sub-intervals were done;
@item 2
excessive roundoff error is detected;
@item 3
extremely bad integrand behavior occurs;
@item 6
if the input is invalid.

@end table

Examples:

@example
(%i1) quad_qawo (x^(-1/2)*exp(-2^(-2)*x), x, 1d-8, 20*2^2, 1, cos);
(%o1)     [1.376043389877692, 4.72710759424899E-11, 765, 0]
(%i2) rectform (integrate (x^(-1/2)*exp(-2^(-alpha)*x) * cos(x), x, 0, inf));
                   alpha/2 - 1/2            2 alpha
        sqrt(%pi) 2              sqrt(sqrt(2        + 1) + 1)
(%o2)   -----------------------------------------------------
                               2 alpha
                         sqrt(2        + 1)
(%i3) ev (%, alpha=2, numer);
(%o3)                     1.376043390090716
@end example

@end deffn

@c THERE ARE OPTIONAL ARGUMENTS WHICH MAKES LISTING THE VARIANTS A LITTLE TEDIOUS
@c NEED A MORE CONVENIENT (AND NONAMBIGUOUS) NOTATION FOR OPTIONAL ARGUMENTS
@deffn {Function} quad_qaws (@var{f(x)}, @var{x}, @var{a}, @var{b}, @var{alpha}, @var{beta}, @var{wfun}, @var{epsabs}, @var{limit})
@deffnx {Function} quad_qaws (@var{f}, @var{x}, @var{a}, @var{b}, @var{alpha}, @var{beta}, @var{wfun}, @var{epsabs}, @var{limit})

Integration of @math{w(x) f(x)} over a finite interval,
where @math{w(x)} is a certain algebraic or logarithmic function.
A globally adaptive subdivision strategy is applied,
with modified Clenshaw-Curtis integration on the subintervals which contain the endpoints
of the interval of integration.

@code{quad_qaws} computes the integral using the Quadpack QAWS
routine:

@ifhtml
@math{integrate (f(x)*w(x), x, a, b)}
@end ifhtml
@ifinfo
@math{integrate (f(x)*w(x), x, a, b)}
@end ifinfo
@tex
$$\int_a^b f(x) w(x) dx$$
@end tex

The weight function @math{w} is selected by @var{wfun}:

@table @code
@item 1
@math{w(x) = (x - a)^alpha (b - x)^beta}
@item 2
@math{w(x) = (x - a)^alpha (b - x)^beta log(x - a)}
@item 3
@math{w(x) = (x - a)^alpha (b - x)^beta log(b - x)}
@item 4
@math{w(x) = (x - a)^alpha (b - x)^beta log(x - a) log(b - x)}
@end table

The integrand may be specified as the name of a Maxima or Lisp function or operator,
a Maxima lambda expression, or a general Maxima expression.

The optional arguments are:

@table @var
@item epsabs
Desired absolute error of approximation.  Default is 1d-10.
@item limit
Size of internal work array.  (@var{limit} - @var{limlst})/2 is the
maximum number of subintervals to use.  Default is 200.
@end table

@c MERGE THESE LINES INTO PRECEDING TABLE
@var{epsabs} and @var{limit} are the desired
relative error and the maximum number of subintervals, respectively.
@var{epsrel} defaults to 1e-8 and @var{limit} is 200.

@code{quad_qaws} returns a list of four elements:

@itemize
@item
an approximation to the integral,
@item
the estimated absolute error of the approximation, 
@item
the number integrand evaluations,
@item
an error code.
@end itemize

The error code (fourth element of the return value) can have the values:

@table @code
@item 0
no problems were encountered; 
@item 1
too many sub-intervals were done;
@item 2
excessive roundoff error is detected;
@item 3
extremely bad integrand behavior occurs;
@item 6
if the input is invalid.

@end table

Examples:

@example
(%i1) quad_qaws (1/(x+1+2^(-4)), x, -1, 1, -0.5, -0.5, 1);
(%o1)     [8.750097361672832, 1.24321522715422E-10, 170, 0]
(%i2) integrate ((1-x*x)^(-1/2)/(x+1+2^(-alpha)), x, -1, 1);
       alpha
Is  4 2      - 1  positive, negative, or zero?

pos;
                          alpha         alpha
                   2 %pi 2      sqrt(2 2      + 1)
(%o2)              -------------------------------
                               alpha
                            4 2      + 2
(%i3) ev (%, alpha=4, numer);
(%o3)                     8.750097361672829
@end example

@end deffn