# Integration

## Introduction to Integration

Maxima has several routines for handling integration. The `integrate` function makes use of most of them. There is also the `antid` package, which handles an unspecified function (and its derivatives, of course). For numerical uses, there is the `romberg` function; an adaptave integrator which uses the Newton-Cotes 8 panel quadrature rule, called `quanc8`; and a set of adaptive integrators from Quadpack, named `quad_qag`, `quad_qags`, etc. Hypergeometric functions are being worked on, see `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 `g(x)` and `h(x)`.

## Definitions for Integration

Function: changevar (expr, f(x,y), y, x)
Makes the change of variable given by `f(x,y) = 0` in all integrals occurring in expr with integration with respect to x. The new variable is y.

```(%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
```

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

`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., `i = j+ ...`, not a higher degree function. E.g.,

```(%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
```

Function: dblint (f, r, s, a, b)
A double-integral routine which was written in top-level Maxima and then translated and compiled to machine code. Use `load (dblint)` to access this package. It uses the Simpson's rule method in both the x and y directions to calculate

```/b /s(x)
|  |
|  |    f(x,y) dy dx
|  |
/a /r(x)
```

The function f must be a translated or compiled function of two variables, and r and s must each be a translated or compiled function of one variable, while a and b must be floating point numbers. The routine has two global variables which determine the number of divisions of the x and y intervals: `dblint_x` and `dblint_y`, both of which are initially 10, and can be changed independently to other integer values (there are `2*dblint_x+1` points computed in the x direction, and `2*dblint_y+1` in the y direction). The routine subdivides the X axis and then for each value of X it first computes `r(x)` and `s(x)`; then the Y axis between `r(x)` and `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 `r(x)` and `s(x)` are, so if the distance `s(x) - 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 `dblint_x` and `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 f, r, and s be either translated or compiled prior to calling `dblint`. This will result in orders of magnitude speed improvement over interpreted code in many cases!

`demo (dblint)` executes a demonstration of `dblint` applied to an example problem.

Function: defint (expr, x, a, b)
Attempts to compute a definite integral. `defint` is called by `integrate` when limits of integration are specified, i.e., when `integrate` is called as `integrate (expr, x, a, b)`. Thus from the user's point of view, it is sufficient to call `integrate`.

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

Function: erf (x)
Represents the error function, whose derivative is: `2*exp(-x^2)/sqrt(%pi)`.

Option variable: erfflag
Default value: `true`

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

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

```(%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
```

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

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

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

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

`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 (`inf` or `minf`), trigonometric functions with limits of integration equal to zero and `%pi` or `2 %pi`, rational functions, integrals related to the definitions of the `beta` and `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 `ldefint`.

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

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

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

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

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

`integrate` is not, by default, declared to be linear. See `declare` and `linear`.

`integrate` attempts integration by parts only in a few special cases.

Examples:

• Elementary indefinite and definite integrals.
```(%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
```
• Use of `assume` and interactive query.
```(%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
```
• Change of variable. There are two changes of variable in this example: one using a derivative established by `gradef`, and one using the derivation `diff(r(x))` of an unspecified function `r(x)`.
```(%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)))
```
• Return value contains the `'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. `grind` shows the noun form `'integrate` in the result. See also `integrate_use_rootsof` for more on integrals of rational functions.
```(%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\$
```
• 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 `f_1` in this example contains the noun form of `integrate`. The double-single-quotes operator `'`' causes the integral to be evaluated, and the result becomes the body of `f_2`.
```(%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
```

System variable: integration_constant_counter
Default value: 0

`integration_constant_counter` is a counter which is updated each time a constant of integration (named by Maxima, e.g., `integrationconstant1`) is introduced into an expression by indefinite integration of an equation.

Option variable: integrate_use_rootsof
Default value: `false`

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

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

```(%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)
```

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:

```(%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)
```

Alternatively the user may compute the roots of the denominator separately, and then express the integrand in terms of these roots, e.g., `1/((x - a)*(x - b)*(x - c))` or `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.

Function: ldefint (expr, x, a, b)
Attempts to compute the definite integral of expr by using `limit` to evaluate the indefinite integral of expr with respect to x at the upper limit b and at the lower limit a. If it fails to compute the definite integral, `ldefint` returns an expression containing limits as noun forms.

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

The calculation makes use of the global variable `potentialzeroloc` which must be `nonlist` or of the form

```[indeterminatej=expressionj, indeterminatek=expressionk, ...]
```

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. `potentialzeroloc` is initially set to 0.

Function: qq
The package `qq` (which may be loaded with `load ("qq")`) contains a function `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 `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 `lo` to `hi` as in `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 `quanc8_relerr` (default value=1.0e-4) and `quanc8_abserr` (default value=1.0e-8) which give the relative error test:

```|integral(function) - computed value| < quanc8_relerr*|integral(function)|
```

and the absolute error test:

```|integral(function) - computed value| < quanc8_abserr
```

`printfile ("qq.usg")` yields additional information.

Function: quanc8 (expr, a, 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 `quanc8`, available in 3 or 4 arg versions. Absolute and relative error checks are used. To use it do `load ("qq")`. See also `qq`.

Function: residue (expr, z, z_0)
Computes the residue in the complex plane of the expression expr when the variable z assumes the value z_0. The residue is the coefficient of `(z - z_0)^(-1)` in the Laurent series for expr.

```(%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
```

Function: risch (expr, x)
Integrates expr with respect to 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 `integrate` can't do. `integrate` will automatically apply `risch` if given these cases.

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

```(%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)
```

Function: romberg (expr, x, a, b)
Function: romberg (expr, a, b)
Romberg integration. There are two ways to use this function. The first is an inefficient way like the definite integral version of `integrate`: `romberg (<integrand>, <variable of integration>, <lower limit>, <upper limit>)`.

Examples:

```(%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
```

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

```romberg (<function name>, <lower limit>, <upper limit>);
```

Continuing the above example, we have:

```(%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
```

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

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

`romberg` may be called recursively and thus can do double and triple integrals.

Example:

```(%i1) assume (x > 0)\$
(%i2) integrate (integrate (x*y/(x+y), y, 0, x/2), x, 1, 3)\$
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
```

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

```(%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
```

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

It is even shorter this way, and the variables do not need to be declared because they are in the context of `romberg`. Use of `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 `rombergtol` and `rombergit` switches.

Option variable: rombergabs
Default value: 0.0

Assuming that successive estimates produced by `romberg` are `y`, `y`, `y`, etc., then `romberg` will return after `n` iterations if (roughly speaking)

``` (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)
```

is `true`. (The condition on the number of iterations given by `rombergmin` must also be satisfied.) Thus if `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

```'integrate (exp(-x), x, 0, 50)
```

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

```(%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
```

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

```(%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
```

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

Option variable: rombergit
Default value: 11

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

Option variable: rombergmin
Default value: 0

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

Option variable: rombergtol
Default value: 1e-4

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

Function: tldefint (expr, x, a, b)
Equivalent to `ldefint` with `tlimswitch` set to `true`.

Function: quad_qag (f(x), x, a, b, key, epsrel, limit)
Numerically evaluate the integral

@math{integrate (f(x), x, a, b)}

The function to be integrated is f(x), with dependent variable x, and the function is to be integrated between the limits a and b. key is the integrator to be used and should be an integer between 1 and 6, inclusive. The value of key selects the order of the Gauss-Kronrod integration rule.

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

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

`quad_qag` returns a list of four elements:

• an approximation to the integral,
• the estimated absolute error of the approximation,
• the number integrand evaluations,
• an error code.

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

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

Examples:

```(%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
```

Function: quad_qags (f(x), x, a, b, epsrel, limit)

Numerically integrate the given function using adaptive quadrature with extrapolation. The function to be integrated is f(x), with dependent variable x, and the function is to be integrated between the limits a and b.

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

`quad_qags` returns a list of four elements:

• an approximation to the integral,
• the estimated absolute error of the approximation,
• the number integrand evaluations,
• an error code.

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

`0`
no problems were encountered;
`1`
too many sub-intervals were done;
`2`
excessive roundoff error is detected;
`3`
`4`
failed to converge
`5`
integral is probably divergent or slowly convergent
`6`
if the input is invalid.

Examples:

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

Note that `quad_qags` is more accurate and efficient than `quad_qag` for this integrand.

Function: quad_qagi (f(x), x, a, inftype, epsrel, limit)

Numerically evaluate one of the following integrals

@math{integrate (f(x), x, a, inf)} @math{integrate (f(x), x, minf, a)} @math{integrate (f(x), x, a, minf, inf)}

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

The parameter inftype determines the integration interval as follows:

`inf`
The interval is from a to positive infinity.
`minf`
The interval is from negative infinity to a.
`both`
The interval is the entire real line.

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

`quad_qagi` returns a list of four elements:

• an approximation to the integral,
• the estimated absolute error of the approximation,
• the number integrand evaluations,
• an error code.

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

`0`
no problems were encountered;
`1`
too many sub-intervals were done;
`2`
excessive roundoff error is detected;
`3`
`4`
failed to converge
`5`
integral is probably divergent or slowly convergent
`6`
if the input is invalid.

Examples:

```(%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
```

Function: quad_qawc (f(x), x, c, a, b, epsrel, limit)

Numerically compute the Cauchy principal value of

@math{integrate (f(x)/(x - c), x, a, b)}

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

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

`quad_qawc` returns a list of four elements:

• an approximation to the integral,
• the estimated absolute error of the approximation,
• the number integrand evaluations,
• an error code.

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

`0`
no problems were encountered;
`1`
too many sub-intervals were done;
`2`
excessive roundoff error is detected;
`3`
`6`
if the input is invalid.

Examples:

```(%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
```

Function: quad_qawf (f(x), x, a, omega, trig, epsabs, limit, maxp1, limlst)

Numerically compute the a Fourier-type integral using the Quadpack QAWF routine. The integral is

@math{integrate (f(x)*w(x), x, a, inf)}

The weight function w is selected by trig:

`cos`
w(x) = cos (omega x)
`sin`
w(x) = sin (omega x)

The optional arguments are:

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

epsabs and limit are the desired relative error and the maximum number of subintervals, respectively. epsrel defaults to 1e-8 and limit is 200.

`quad_qawf` returns a list of four elements:

• an approximation to the integral,
• the estimated absolute error of the approximation,
• the number integrand evaluations,
• an error code.

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

`0`
no problems were encountered;
`1`
too many sub-intervals were done;
`2`
excessive roundoff error is detected;
`3`
`6`
if the input is invalid.

Examples:

```(%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
```

Function: quad_qawo (f(x), x, a, b, omega, trig, epsabs, limit, maxp1, limlst)

Numerically compute the integral using the Quadpack QAWO routine:

@math{integrate (f(x)*w(x), x, a, b)}

The weight function w is selected by trig:

`cos`
w(x) = cos (omega x)
`sin`
w(x) = sin (omega x)

The optional arguments are:

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

epsabs and limit are the desired relative error and the maximum number of subintervals, respectively. epsrel defaults to 1e-8 and limit is 200.

`quad_qawo` returns a list of four elements:

• an approximation to the integral,
• the estimated absolute error of the approximation,
• the number integrand evaluations,
• an error code.

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

`0`
no problems were encountered;
`1`
too many sub-intervals were done;
`2`
excessive roundoff error is detected;
`3`
`6`
if the input is invalid.

Examples:

```(%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
```

Function: quad_qaws (f(x), x, a, b, alfa, beta, wfun, epsabs, limit)

Numerically compute the integral using the Quadpack QAWS routine:

@math{integrate (f(x)*w(x), x, a, b)}

The weight function w is selected by wfun:

`1`
w(x) = (x - a)^alfa (b - x)^beta
`2`
w(x) = (x - a)^alfa (b - x)^beta log(x - a)
`3`
w(x) = (x - a)^alfa (b - x)^beta log(b - x)
`2`
w(x) = (x - a)^alfa (b - x)^beta log(x - a) log(b - x)

The optional arguments are:

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

epsabs and limit are the desired relative error and the maximum number of subintervals, respectively. epsrel defaults to 1e-8 and limit is 200.

`quad_qaws` returns a list of four elements:

• an approximation to the integral,
• the estimated absolute error of the approximation,
• the number integrand evaluations,
• an error code.

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

`0`
no problems were encountered;
`1`
too many sub-intervals were done;
`2`
excessive roundoff error is detected;
`3`
`6`
if the input is invalid.

Examples:

```(%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
```