[ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
20.1 Introduction to Integration | ||
20.2 Definitions for Integration |
[ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
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)
.
[ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
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 |
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.
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.
Represents the error function, whose derivative is:
2*exp(-x^2)/sqrt(%pi)
.
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.
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 |
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:
(%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 |
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 |
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))) |
'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$ |
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 |
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.
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.
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[0]
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.
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.
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
.
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 |
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) |
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)$ (%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 |
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 |
The advantage here is shortness.
(%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.
Default value: 0.0
Assuming that successive estimates
produced by romberg
are y[0]
, y[1]
, y[2]
, 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.
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.
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.
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.
Equivalent to ldefint
with tlimswitch
set to true
.
Numerically evaluate the integral
integrate (f(x), x, a, b)
using a simple adaptive integrator.
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:
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 |
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:
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
extremely bad integrand behavior occurs;
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.
Numerically evaluate one of the following integrals
integrate (f(x), x, a, inf)
integrate (f(x), x, minf, a)
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:
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
extremely bad integrand behavior occurs;
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 |
Numerically compute the Cauchy principal value of
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:
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
extremely bad integrand behavior occurs;
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 |
Numerically compute the a Fourier-type integral using the Quadpack QAWF routine. The integral is
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:
Desired absolute error of approximation. Default is 1d-10.
Size of internal work array. (limit - limlst)/2 is the maximum number of subintervals to use. Default is 200.
Maximum number of Chebyshev moments. Must be greater than 0. Default is 100.
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:
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
extremely bad integrand behavior occurs;
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 |
Numerically compute the integral using the Quadpack QAWO routine:
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:
Desired absolute error of approximation. Default is 1d-10.
Size of internal work array. (limit - limlst)/2 is the maximum number of subintervals to use. Default is 200.
Maximum number of Chebyshev moments. Must be greater than 0. Default is 100.
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:
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
extremely bad integrand behavior occurs;
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 |
Numerically compute the integral using the Quadpack QAWS routine:
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:
Desired absolute error of approximation. Default is 1d-10.
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:
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
extremely bad integrand behavior occurs;
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 |
[ << ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
This document was generated on March, 19 2006 using texi2html 1.76.