 Numerical Integration - Maple Programming Help

Numerical Integration

Calling Sequence

 evalf(Int(f, x=a..b, ...)) $\mathrm{evalf}\left({{\int }}_{a}^{b}f\phantom{\rule[-0.0ex]{0.3em}{0.0ex}}{ⅆ}x\right)$ evalf(Int(f, a..b, ...)) evalf(Int(f, list-of-equations, ...)) evalf(Int(f, list-of-ranges, ...)) evalf(int(f, x=a..b)) $\mathrm{evalf}\left({\int }_{a}^{b}f\phantom{\rule[-0.0ex]{0.3em}{0.0ex}}ⅆx\right)$

Parameters

 f - algebraic expression or procedure; the integrand x - name; the variable of integration a, b - endpoints of the interval of integration list-of-equations - list of equations [x1=a1..b1, ..., xn=an..bn] list-of-ranges - list of ranges [a1..b1, ..., an..bn] ... - (optional) zero or more options, as described below

Description

 • The most common command for numerical integration is evalf(Int(f, x=a..b)) where the integration command is expressed in inert form to avoid first invoking the symbolic integration routines. It is also possible to invoke evalf on an unevaluated integral returned by the symbolic int command, as in evalf(int(f, x=a..b)), if it happens that symbolic int fails (returns an unevaluated integral).
 • All numerical integration calling sequences can also be accessed directly from the int command by using the numeric option.
 • You can enter the command evalf/Int using either the 1-D or 2-D calling sequence. For example, evalf(Int(1/(x^2+1), x=0..infinity)) is equivalent to $\mathrm{evalf}\left({{\int }}_{0}^{\mathrm{\infty }}\frac{1}{{x}^{2}+1}\phantom{\rule[-0.0ex]{0.3em}{0.0ex}}{ⅆ}x\right)$.
 • The integrand f may be another unevaluated integral, that is, multiple integrals are supported. A special list syntax (see below) can be used to specify multiple integrals, rather than using nested integrals. Integrals expressed in the standard non-list notation are referred to as 1-D (one-dimensional) integrals including the case of nested 1-D integrals.
 • If the integrand f is specified as a procedure or a Maple operator, then the second argument must be a range a..b and not an equation, that is, a variable of integration must not be specified.
 • Various levels of user information are displayed during the computation if infolevel[evalf/int] is assigned values between 1 and 4.

Optional Arguments

 • Additional options may be specified as equations. (For backward compatibility some options are accepted as values rather than equations, as specified below.) An option is one of the following forms:

 method =  or   digits =  or   epsilon = methodoptions = maxintervals =

 • The specification method =  (or simply ) indicates a particular numerical integration method to be applied. The methods that can be specified are described below. By default, a hybrid symbolic-numeric strategy is applied.
 • The specification digits =  (or simply ) indicates the number of digits of precision for the computation. Some additional guard digits are carried during the computation to attempt to achieve a result with  correct digits (although a larger tolerance can be specified by using the 'epsilon' option). By default, the Maple environment variable Digits specifies the precision for the computation.
 • The specification epsilon =  specifies the relative error tolerance for the computed result. The routines attempt to achieve a final result with a relative error less than this value. By default, the relative error tolerance which the routines attempt to achieve for the final result is

$\mathrm{eps}=0.5{10}^{1-\mathrm{digits}}$

 where digits is the precision specified for the computation. In attempting to achieve this accuracy, the working value of Digits is increased as deemed necessary. It is an error to specify 'epsilon' smaller than the default value above, and for any value larger than 1e-3 the value 1e-3 is used instead if the method in use is deterministic (i.e. not the MonteCarlo or Cuba methods).
 Note: For some integrands, the numerical accuracy attained when computing values of the integrand may be insufficient to allow the value of the integral to be computed to the default tolerance $\mathrm{eps}$ (even though the computation is using some number of guard digits). In such cases, specifying a larger tolerance (relative to the setting of digits) via the 'epsilon' option may be helpful. Alternatively, increasing Digits and fixing 'epsilon' may provide the desired answer (see the end of the examples section).
 • The specification methodoptions =  specifies a list of zero or more options that are specific to a method selected with the method option. In particular, if the _d01ajc or _d01akc method is selected, one can supply an option of the form methodoptions=[maxintervals = ] to specify a maximal number of subintervals that can be used internally by those methods.
 • For backward compatibility, the option maxintervals =  for the _d01ajc and _d01akc methods can also be specified as a separate option, as an argument to Int directly, rather than in the methodoptions option.

Outline of the Numerical Integration Polyalgorithm (1-D Integrals)

 • In the default case (no particular method specified), the problem is first passed to NAG integration routines if Digits is not too large (that is, if Digits <= evalhf(Digits)). The NAG routines are in a compiled C library and hence operate at hardware floating-point speed. If the NAG routines cannot perform the integration, then some singularity handling may be performed and control may pass back to the NAG routines with a modified problem. Native Maple routines are invoked if the NAG routines cannot solve the problem (for example, if Digits is too large or if the integrand involves functions for which hardware floating-point evaluation is not supported).
 • The native Maple hybrid symbolic-numeric solution strategy is as follows. The default numerical method applied is Clenshaw-Curtis quadrature (_CCquad). If slow convergence is detected then there must be singularities in or near the interval of integration (perhaps in the complex plane). Some techniques of symbolic analysis are used to deal with the singularities. For problems with non-removable endpoint singularities, an adaptive double-exponential quadrature method (_Dexp) is applied.
 • If singularities interior to the interval are suspected, then an attempt is made to locate the singularities in order to break up the interval of integration. Finally, if still unsuccessful, then the interval is subdivided and the _Dexp method is applied, or if the method was already _Dexp or _Sinc then an adaptive Gaussian quadrature method (_Gquad) is applied.
 • For the limits of integration, the values infinity and/or -infinity are valid, and a symbolic-numeric strategy attempts to deal with singularities. Techniques employed include variable transformations, subtracting out the singularity, and integration of a truncated generalized series near the singularity.
 • No singularity handling is attempted in the case where the integrand f is specified as a procedure or a Maple operator.

Special (List) Syntax for Multiple Integrals

 • A numerical multiple integration problem may be specified in a natural way using nested one-dimensional integrals, for example:

 evalf( Int(...(Int(Int(f, x1=a1..b1), x2=a2..b2), ...), xn=an..bn) )

 where the integrand f depends on x1, x2, ..., xn. Such a problem may also be specified using the following special multiple integration notation with a list as the second argument:

 evalf( Int(f, [x1=a1..b1, x2=a2..b2, ..., xn=an..bn]) ) .

 • Additional optional arguments may be stated just as in the case of 1-D integration. Also as in 1-D integration, the integrand f may be specified as a procedure in which case the second argument must be a list of ranges: [a1..b1, a2..b2, ..., an..bn].
 • Whether a multiple integration problem is stated using nested integrals or using the list notation, the arguments will be extracted so as to invoke the same numerical multiple integration routines.

The Method Names

 • The optional argument method =  (or simply ) accepts the following method names.

 method = _DEFAULT -- equivalent to not specifying a method; the solution strategy outlined above is applied for 1-D integrals; for multiple integrals, the problem is passed to the _cuhre method and if it fails, then the problem is treated via nested 1-D integration. method = _NoNAG -- indicates to avoid calling NAG routines; otherwise follow the _DEFAULT strategy. method = _NoMultiple -- indicates to avoid calling numerical multiple integration routines; compute multiple integrals via nested 1-D integration.

Maple Methods

 • Specifying a method indicates to try only that method (in particular, no NAG methods and no singularity handling).

 method = _CCquad -- Clenshaw-Curtis quadrature method. method = _Dexp -- adaptive double-exponential method. method = _Gquad -- adaptive Gaussian quadrature method. method = _Sinc -- adaptive sinc quadrature method. method = _NCrule -- adaptive Newton-Cotes method "quanc8". Note that in contrast to the other Maple methods listed here, "quanc8" (method = _NCrule) is a fixed-order method and hence it is not recommended for very high precisions (e.g. Digits > 15).

NAG Methods

 • Specifying a method indicates to try only that method (in particular, no singularity handling and no Maple methods).

 method = _d01ajc -- for finite interval of integration; allows for badly behaved integrands; uses adaptive Gauss 10-point and Kronrod 21-point rules. method = _d01akc -- for finite interval of integration, oscillating integrands; uses adaptive Gauss 30-point and Kronrod 61-point rules. method = _d01amc -- for semi-infinite/infinite interval of integration.

Multiple Integration Methods

 • These methods are for multiple integrals over a hyperrectangle, that is, the limits of integration are finite constants.
 • Specifying a method indicates to try only that method (in particular, do not revert to nested 1-D integration).

 method = _cuhre -- dimensions 2 to 15; ACM TOMS Algorithm 698. method = _MonteCarlo -- Monte Carlo method; for low accuracy only (less than 5 digits of accuracy); NAG routine 'd01gbc'. method = _CubaVegas -- Vegas method; for low accuracy. For details and method-specific options, see evalf/Int/cuba. method = _CubaSuave -- Suave method; for low accuracy. For details and method-specific options, see evalf/Int/cuba. method = _CubaDivonne -- Divonne method; for low accuracy. For details and method-specific options, see evalf/Int/cuba. method = _CubaCuhre -- Cuhre method. For details and method-specific options, see evalf/Int/cuba.

Examples

 > $\mathrm{evalf}\left(\mathrm{Int}\left(\frac{\mathrm{exp}\left(-{x}^{3}\right)}{{x}^{2}+1},x=0..1\right)\right)$
 ${0.6649369431}$ (1)
 > $\mathrm{evalf}\left(\mathrm{Int}\left(\frac{1}{{x}^{2}+1},x=0..\mathrm{\infty }\right)\right)$
 ${1.570796327}$ (2)
 > $\mathrm{evalf}\left(\mathrm{Int}\left(\mathrm{sin}\left(x\right)\mathrm{ln}\left(x\right)\mathrm{exp}\left(-{x}^{3}\right),x=0..\mathrm{\infty }\right)\right)$
 ${-0.1957885158}$ (3)

The following integrals are computed to higher precision.

 > $\mathrm{e1}≔\frac{1}{\mathrm{\Gamma }\left(x\right)}:$
 > $\mathrm{Int}\left(\mathrm{e1},x=0..2\right)=\mathrm{evalf}\left(\mathrm{Int}\left(\mathrm{e1},x=0..2,\mathrm{digits}=20,\mathrm{method}=\mathrm{_Dexp}\right)\right)$
 ${{\int }}_{{0}}^{{2}}\frac{{1}}{{\mathrm{\Gamma }}{}\left({x}\right)}\phantom{\rule[-0.0ex]{0.3em}{0.0ex}}{ⅆ}{x}{=}{1.6263783986861406145}$ (4)
 > $\mathrm{e2}≔\frac{\mathrm{exp}\left(v-\frac{{v}^{2}}{2}\right)}{1+\frac{1}{2}\mathrm{exp}\left(v\right)}:$
 > $\mathrm{Int}\left(\mathrm{e2},v=0..\mathrm{\infty }\right)=\mathrm{evalf}\left[20\right]\left(\mathrm{Int}\left(\mathrm{e2},v=0..\mathrm{\infty }\right)\right)$
 ${{\int }}_{{0}}^{{\mathrm{\infty }}}\frac{{{ⅇ}}^{{v}{-}\frac{{1}}{{2}}{}{{v}}^{{2}}}}{{1}{+}\frac{{{ⅇ}}^{{v}}}{{2}}}\phantom{\rule[-0.0ex]{0.3em}{0.0ex}}{ⅆ}{v}{=}{1.3055168991185060654}$ (5)
 > $\mathrm{e3}≔\frac{1}{1+\mathrm{ln}\left(1+x\right)}:$
 > $\mathrm{Int}\left(\mathrm{e3},x=0..1\right)=\mathrm{evalf}\left[32\right]\left(\mathrm{Int}\left(\mathrm{e3},x=0..1,\mathrm{method}=\mathrm{_Gquad}\right)\right)$
 ${{\int }}_{{0}}^{{1}}\frac{{1}}{{1}{+}{\mathrm{ln}}{}\left({1}{+}{x}\right)}\phantom{\rule[-0.0ex]{0.3em}{0.0ex}}{ⅆ}{x}{=}{0.73716070962368003213791626905536}$ (6)
 > $r≔\mathrm{int}\left(\mathrm{sech}\left(x\right)\mathrm{exp}\left(-{x}^{2}\right),x=-\mathrm{\infty }..\mathrm{\infty }\right)$
 ${r}{≔}{{\int }}_{{-}{\mathrm{\infty }}}^{{\mathrm{\infty }}}{\mathrm{sech}}{}\left({x}\right){}{{ⅇ}}^{{-}{{x}}^{{2}}}\phantom{\rule[-0.0ex]{0.3em}{0.0ex}}{ⅆ}{x}$ (7)
 > $\mathrm{evalf}\left(r\right)$
 ${1.479061171}$ (8)
 > $\mathrm{evalf}\left[25\right]\left(r\right)$
 ${1.479061171449575890854454}$ (9)

The following command returns an error because procedure $f$ is invoked with argument $x$, a symbolic name.

 > f := proc(x) if x < 2 then 2*x else x^2 end if; end proc;
 ${f}{:=}{\mathbf{proc}}\left({x}\right)\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathbf{if}}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{x}{<}{2}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathbf{then}}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{2}{*}{x}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathbf{else}}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{x}{^}{2}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathbf{end if}}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathbf{end proc}}$ (10)
 > $\mathrm{evalf}\left(\mathrm{Int}\left(f\left(x\right),x=0..3\right)\right)$

When the integrand $f$ is a procedure, the following syntax should be used.

 > $\mathrm{evalf}\left(\mathrm{Int}\left(f,0..3\right)\right)$
 ${10.33333333}$ (11)

Note that the following command also works by delaying the evaluation of $f\left(x\right)$ via unevaluation quotes.

 > $\mathrm{evalf}\left(\mathrm{Int}\left('f'\left(x\right),x=0..3\right)\right)$
 ${10.33333333}$ (12)

Multiple integrals may be expressed as nested one-dimensional integrals.

 > $\mathrm{Int}\left(\mathrm{Int}\left(\mathrm{Int}\left(\mathrm{/}\left(\mathrm{exp}\left(\mathrm{+}\left(\mathrm{+}\left(x,y\right),z\right)\right),\mathrm{*}\left(\mathrm{*}\left(\mathrm{+}\left(\mathrm{*}\left(5,x\right),1\right),\mathrm{+}\left(\mathrm{*}\left(10,y\right),2\right)\right),\mathrm{+}\left(\mathrm{*}\left(15,z\right),3\right)\right)\right),\mathrm{=}\left(x,0..4\right)\right),\mathrm{=}\left(y,0..3\right)\right),\mathrm{=}\left(z,0..\mathrm{sqrt}\left(2\right)\right)\right)$
 ${?}$ (13)
 > $\mathrm{evalf}\left(\right)$
 ${0.9331611325}$ (14)

Numerical multiple integration may also be invoked using a list syntax.

 > $d≔1-{w}^{2}{x}^{2}{y}^{2}{z}^{2}:$
 > $g≔d\mathrm{cos}\left(wxyz\right)-dwxyz\mathrm{sin}\left(wxyz\right)$
 ${g}{≔}\left({-}{{w}}^{{2}}{}{{x}}^{{2}}{}{{y}}^{{2}}{}{{z}}^{{2}}{+}{1}\right){}{\mathrm{cos}}{}\left({w}{}{x}{}{y}{}{z}\right){-}\left({-}{{w}}^{{2}}{}{{x}}^{{2}}{}{{y}}^{{2}}{}{{z}}^{{2}}{+}{1}\right){}{w}{}{x}{}{y}{}{z}{}{\mathrm{sin}}{}\left({w}{}{x}{}{y}{}{z}\right)$ (15)
 > $\mathrm{evalf}\left(\mathrm{Int}\left(g,\left[w=0..1,x=0..1,y=0..1,z=0..1\right]\right)\right)$
 ${0.9717798177}$ (16)

When low accuracy is sufficient, the Monte Carlo method may be used.

 > $h≔\frac{1}{2+\mathrm{sin}\left(\mathrm{\pi }\mathrm{sqrt}\left(87\right)\left(\mathrm{x1}+\mathrm{x2}+\mathrm{x3}+\mathrm{x4}+\mathrm{x5}+\mathrm{x6}\right)\right)}$
 ${h}{≔}\frac{{1}}{{2}{+}{\mathrm{sin}}{}\left({\mathrm{\pi }}{}\sqrt{{87}}{}\left({\mathrm{x1}}{+}{\mathrm{x2}}{+}{\mathrm{x3}}{+}{\mathrm{x4}}{+}{\mathrm{x5}}{+}{\mathrm{x6}}\right)\right)}$ (17)
 > $\mathrm{evalf}\left(\mathrm{Int}\left(h,\left[\mathrm{x1}=-1..1,\mathrm{x2}=-1..1,\mathrm{x3}=-1..1,\mathrm{x4}=-1..1,\mathrm{x5}=-1..1,\mathrm{x6}=-1..1\right],\mathrm{method}=\mathrm{_MonteCarlo},\mathrm{\epsilon }=0.005\right)\right)$
 ${36.91495206}$ (18)

Only trust about 3 digits when epsilon = 0.5e-2.

 > $\mathrm{evalf}\left[3\right]\left(\right)$
 ${36.9}$ (19)

The following integrand has a region near x=0.5 where evaluation incurs catastrophic cancellation to the extent that the function cannot even be evaluated to 1 significant Digit at standard precision.

 > $\mathrm{igrand}≔\frac{1}{2-\mathrm{sin}\left(\mathrm{\pi }x\right)-\mathrm{sin}\left(\frac{355}{113}x\right)}$
 ${\mathrm{igrand}}{≔}\frac{{1}}{{2}{-}{\mathrm{sin}}{}\left({\mathrm{\pi }}{}{x}\right){-}{\mathrm{sin}}{}\left(\frac{{355}{}{x}}{{113}}\right)}$ (20)
 > $\mathrm{evalf}\left(\mathrm{eval}\left(\mathrm{igrand},x=0.5\right)\right)$
 ${Float}{}\left({\mathrm{\infty }}\right)$ (21)

Note that evalf fails to compute this integral with default settings.

 > $\mathrm{evalf}\left(\mathrm{Int}\left(\mathrm{igrand},x=0..1\right)\right)$
 ${{\int }}_{{0.}}^{{1.}}\frac{{1}}{{2.}{-}{1.}{}{\mathrm{sin}}{}\left({3.141592654}{}{x}\right){-}{1.}{}{\mathrm{sin}}{}\left({3.141592920}{}{x}\right)}\phantom{\rule[-0.0ex]{0.3em}{0.0ex}}{ⅆ}{x}$ (22)

So to compute the value of this integral to 10 digits, we need to add significant guard digits:

 > $\mathrm{evalf}\left[15\right]\left(\mathrm{eval}\left(\mathrm{igrand},x=0.5\right)\right)$
 ${1.11111111111111}{}{{10}}^{{14}}$ (23)
 > $\mathrm{evalf}\left[20\right]\left(\mathrm{eval}\left(\mathrm{igrand},x=0.5\right)\right)$
 ${1.1241778044582643369}{}{{10}}^{{14}}$ (24)
 > $\mathrm{evalf}\left[25\right]\left(\mathrm{eval}\left(\mathrm{igrand},x=0.5\right)\right)$
 ${1.124177605944406775029070}{}{{10}}^{{14}}$ (25)

So we add 15 digits to assure we get the answer to 10 digits:

 > $\mathrm{evalf}\left[25\right]\left(\mathrm{Int}\left(\mathrm{igrand},x=0..1,\mathrm{\epsilon }=1.{10}^{-10}\right)\right)$
 ${1.499451605234141071490295}{}{{10}}^{{7}}$ (26)

In the following example, the default setting for the maximum number of subintervals, $500$ in this case, is not enough for successful integration using the NAG method for oscillatory integrands, $\mathrm{_d01akc}$, which would be suitable for this integrand. By supplying a higher upper bound, we can get successful completion with this method.

 > $\mathrm{igrand}≔\frac{\mathrm{sin}\left(\mathrm{exp}\left(\mathrm{abs}\left(x\right)\right)\right)}{1+{x}^{2}}$
 ${\mathrm{igrand}}{≔}\frac{{\mathrm{sin}}{}\left({{ⅇ}}^{\left|{x}\right|}\right)}{{{x}}^{{2}}{+}{1}}$ (27)
 > $\mathrm{evalf}\left(\mathrm{Int}\left(\mathrm{igrand},x=-10..10,\mathrm{method}=\mathrm{_d01akc}\right)\right)$
 > $\mathrm{evalf}\left(\mathrm{Int}\left(\mathrm{igrand},x=-10..10,\mathrm{method}=\mathrm{_d01akc},\mathrm{methodoptions}=\left[\mathrm{maxintervals}=1000\right]\right)\right)$
 ${1.235076653}$ (28)