Go to the first, previous, next, last section, table of contents.
MACSYMA has several routines for handling integration. The INTEGRATE command 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, do DESCRIBE(SPECINT); for details. Generally speaking, MACSYMA 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).
(%i1) 'INTEGRATE(%E**SQRT(A*Y),Y,0,4); 4 / [ SQRT(A) SQRT(Y) (%o1) I (%E ) dY ] / 0 (%i2) CHANGEVAR(%o1,Y-Z^2/A,Z,Y); 2 SQRT(A) / [ Z 2 I Z %E dZ ] / 0 (%o4) --------------------- A
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.
(%i3) SUM(A[I]*X^(I-2),I,0,INF); INF ==== \ I - 2 (%o3) > A X / I ==== I = 0 (%i4) CHANGEVAR(%,I-2-N,N,I); INF ==== \ N (%o4) > A X / N + 2 ==== N = - 2
EXP(A*X+B)*COS(C*X)^N*SIN(C*X)
then the regular integration program will be invoked if the switch ERRINTSCE[TRUE] is TRUE. If it is FALSE then INTSCE will err out.
(%i1) 'INTEGRATE(SINH(A*X)*F(T-X),X,0,T)+B*F(T)=T**2; T / [ 2 (%o1) I (SINH(A X) F(T - X)) dX + B F(T) = T ] / 0 (%i2) LAPLACE(%,T,S); A LAPLACE(F(T), T, S) (%o2) --------------------- 2 2 S - A 2 + B LAPLACE(F(T), T, S) = -- 3 S (%i3) LINSOLVE([%],['LAPLACE(F(T),T,S)]); SOLUTION 2 2 2 S - 2 A (%t3) LAPLACE(F(T), T, S) = -------------------- 5 2 3 B S + (A - A B) S (%o3) [%t3] (%i4) ILT(%t3,S,T); IS A B (A B - 1) POSITIVE, NEGATIVE, OR ZERO? POS; 2 SQRT(A) SQRT(A B - B) T 2 COSH(------------------------) B (%o4) F(T) = - -------------------------------- A 2 A T 2 + ------- + ------------------ A B - 1 3 2 2 A B - 2 A B + A
(%i1) INTEGRATE(SIN(X)**3,X); 3 COS (X) (%o1) ------- - COS(X) 3 (%i2) INTEGRATE(X**A/(X+1)**(5/2),X,0,INF); IS A + 1 POSITIVE, NEGATIVE, OR ZERO? POS; IS 2 A - 3 POSITIVE, NEGATIVE, OR ZERO? NEG; 3 (%o2) BETA(A + 1, - - A) 2 (%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)))
(Note 1) The fact that MACSYMA does not perform certain integrals does
not always imply that the integral does not exist in closed form. In
the example below the integration call returns the noun form but the
integral can be found fairly easily. For example, one can compute the
roots of X^3+X+1 = 0
to rewrite the integrand in the form
1/((X-A)*(X-B)*(X-C))
where A, B and C are the roots. MACSYMA will integrate this equivalent form although the integral is quite complicated.
(%i6) INTEGRATE(1/(X^3+X+1),X); / [ 1 (%o6) I ---------- dX ] 3 / X + X + 1
(%i4) 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) (%o4) ----------------- - --------------- + --------------- 7 14 7 SQRT(3)
but now we set the flag to be true and the first part of the integral will undergo further simplification.
(%i5) INTEGRATE_USE_ROOTSOF:true; (%o5) TRUE
(%i6) integrate(1/(1+x+x^5),x); ==== 2 \ (%R1 - 4 %R1 + 5) LOG(x - %R1) > ------------------------------- / 2 ==== 3 %R1 - 2 %R1 3 2 %R1 in ROOTSOF(x - x + 1) (%o6) ---------------------------------------------------------- 7 2 x + 1 2 5 ATAN(-------) LOG(x + x + 1) SQRT(3) - --------------- + --------------- 14 7 SQRT(3)
Note that it may be that we want to approximate the roots in the complex plane, and then provide the function factored, since we will then be able to group the roots and their complex conjugates, so as to give a better answer.
EXP(A*X+B)*COS(C*X)^N*SIN(C*X)^M
The call is INTSCE(expr,var) expr may be any expression, but if it is not in the above form then the regular integration program will be invoked if the switch ERRINTSCE[TRUE] is TRUE. If it is FALSE then INTSCE will err out.
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.
(%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
(%i1) RISCH(X^2*ERF(X),X); 2 2 - X X 3 2 %E (%E SQRT(%PI) X ERF(X) + X + 1) (%o1) ------------------------------------------ 3 SQRT(%PI) (%i2) DIFF(%,X),RATSIMP; 2 (%o2) X ERF(X)
Examples: ROMBERG(SIN(Y),Y,1,%PI); TIME= 39 MSEC. 1.5403023 F(X):=1/(X^5+X+1); ROMBERG(F(X),X,1.5,0); TIME= 162 MSEC. - 0.75293843
The second is an efficient way that is used as follows:
ROMBERG(<function name>,<lower limit>,<upper limit>);
Example: F(X):=(MODE_DECLARE([FUNCTION(F),X],FLOAT),1/(X^5+X+1)); TRANSLATE(F); ROMBERG(F,1.5,0); TIME= 13 MSEC. - 0.75293843
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, do DESCRIBE(ROMBERGABS,ROMBERGMIN); for details. ROMBERG may be called recursively and thus can do double and triple integrals.
Example: INTEGRATE(INTEGRATE(X*Y/(X+Y),Y,0,X/2),X,1,3); 13/3 (2 LOG(2/3) + 1) %,NUMER; 0.81930233 DEFINE_VARIABLE(X,0.0,FLOAT,"Global variable in function F")$ F(Y):=(MODE_DECLARE(Y,FLOAT), X*Y/(X+Y) )$ G(X):=ROMBERG('F,0,X/2)$ ROMBERG(G,1,3); 0.8193023
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:
G1(X):=(MODE_DECLARE(X,FLOAT), ROMBERG(X*Y/(X+Y),Y,0,X/2))$ ROMBERG(G1,1,3); 0.8193023
The advantage here is shortness.
Q(A,B):=ROMBERG(ROMBERG(X*Y/(X+Y),Y,0,X/2),X,A,B)$ Q(1,3); 0.8193023
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.
Integral(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.
F(X):=(MODE_DECLARE(N,INTEGER,X,FLOAT),N:N+1,EXP(-X))$ TRANSLATE(F)$ /* First of all try doing the whole integral at once */ BLOCK([ROMBERGTOL:1.E-6,ROMBERABS:0.],N:0,ROMBERG(F,0,50)); ==> 1.00000003 N; ==> 257 /* Number of function evaluations*/
Now do the integral intelligently, by first doing Integral(exp(-x),x,0,10) and then setting ROMBERGABS to 1.E-6*(this partial integral).
BLOCK([ROMBERGTOL:1.E-6,ROMBERGABS:0.,SUM:0.], N:0,SUM:ROMBERG(F,0,10),ROMBERGABS:SUM*ROMBERGTOL,ROMBERGTOL:0., SUM+ROMBERG(F,10,50)); ==> 1.00000001 /* Same as before */ N; ==> 130
So if F(X) were a function that took a long time to compute, the second method would be about 2 times quicker.
using a simple adaptive integrator.
The function to be integrated is fun
, with dependent
variable var
, 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.
A list is returned containing the following elements:
The error code is
0
1
2
3
6
Examples:
(C1) quad_qag(x^(1/2)*log(1/x),x,0,1,3); (D1) [.4444444444492108, 3.170096850288723e-9, 961, 0] (C2) integrate(x^(1/2)*log(1/x),x,0,1); 4 (D2) - 9
Numerically integrate the given function using adaptive quadrature with
extrapolation. The function to be integrated is fun
, with
dependent variable var
, 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.
A list is returned containing the following elements:
The error code is
0
1
2
3
4
5
6
Examples:
(C1) quad_qags(x^(1/2)*log(1/x), x, 0 ,1); (D1) [.4444444444444448, 1.11022302462516e-15, 315, 0]
We that QAGS is more accurate and efficient than QAG for this integrand.
Numerically evaluate one of the following integrals
using the Quadpack QAGI routine. The function to be integrated is
fun
, with dependent variable var
, and the function is to
be integrated over an infinite range.
The parameter inftype
determines the integration interval as
follows:
inf
a
to positive infinity.
minf
a
.
both
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.
A list is returned containing the following elements:
The error code is
0
1
2
3
4
5
6
Examples:
(C1) quad_qagi(x^2*exp(-4*x),x,0,inf); (D1) [0.03125, 2.95916102995002e-11, 105, 0] (C2) integrate(x^2*exp(-4*x),x,0,inf); 1 (D2) -- 32
Numerically compute the Cauchy principal value of
using the Quadpack QAWC routine. The function to be integrated is
fun/(x-c)
, with dependent variable var
, 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.
A list is returned containing the following elements:
The error code is
0
1
2
3
6
Examples:
(C1) quad_qawc(2^(-5)*((x-1)^2+4^(-5))^(-1), x, 2, 0, 5); (D2) [- 3.130120337415925, 1.306830140249558e-8, 495, 0] (C2) integrate(2^(-alpha)*(((x-1)^2 + 4^(-alpha))*(x-2))^(-1), x, 0, 5); ALPHA ALPHA 9 4 9 4 LOG(------------- + -------------) ALPHA ALPHA 64 4 + 4 64 4 + 4 (D3) (----------------------------------------- ALPHA 2 4 + 2 3 ALPHA ------- 2 ALPHA/2 2 4 ATAN(4 4 ) - --------------------------- ALPHA 2 4 + 2 3 ALPHA ------- 2 ALPHA/2 2 4 ATAN(4 ) ALPHA - -------------------------)/2 ALPHA 2 4 + 2 (C3) ev(%,alpha=5,numer); (D3) - 3.130120337415917
Numerically compute the a Fourier-type integral using the Quadpack QAWF routine. The integral is
The weight function w(x) is selected by trig
:
cos
sin
The optional arguments are:
epsabs
limit
limit
-limlst
)/2 is the
maximum number of subintervals to use. Default is 200.
maxp1
limlst
epsabs
and limit
are the desired
relative error and the maximum number of subintervals, respectively.
epsrel
defaults to 1e-8 and limit
is 200.
A list is returned containing the following elements:
The error code is
0
1
2
3
6
Examples:
(C1) quad_qawf(2^(-5)*((x-1)^2+4^(-5))^(-1), x, 2, 0, 5); (D2) [- 3.130120337415925, 1.306830140249558e-8, 495, 0] (C2) integrate(2^(-alpha)*(((x-1)^2 + 4^(-alpha))*(x-2))^(-1), x, 0, 5); ALPHA ALPHA 9 4 9 4 LOG(------------- + -------------) ALPHA ALPHA 64 4 + 4 64 4 + 4 (D3) (----------------------------------------- ALPHA 2 4 + 2 3 ALPHA ------- 2 ALPHA/2 2 4 ATAN(4 4 ) - --------------------------- ALPHA 2 4 + 2 3 ALPHA ------- 2 ALPHA/2 2 4 ATAN(4 ) ALPHA - -------------------------)/2 ALPHA 2 4 + 2 (C3) ev(%,alpha=5,numer); (D3) - 3.130120337415917
Numerically compute the integral using the Quadpack QAWO routine:
The weight function w(x) is selected by trig
:
cos
sin
The optional arguments are:
epsabs
limit
limit
-limlst
)/2 is the
maximum number of subintervals to use. Default is 200.
maxp1
limlst
epsabs
and limit
are the desired
relative error and the maximum number of subintervals, respectively.
epsrel
defaults to 1e-8 and limit
is 200.
A list is returned containing the following elements:
The error code is
0
1
2
3
6
Examples:
(C1) quad_qawo(x^(-1/2)*exp(-2^(-2)*x), x, 1d-8, 20*2^2, 1, cos); (D2) [1.376043389877607, 4.72710759424899e-11, 765, 0] (C2) 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) (D29) ----------------------------------------------------- 2 ALPHA SQRT(2 + 1) (C3) ev(%,alpha=2,numer); (D3) 1.376043390090716
Numerically compute the integral using the Quadpack QAWS routine:
The weight function w(x) is selected by wfun
:
1
2
3
2
The optional arguments are:
epsabs
limit
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.
A list is returned containing the following elements:
The error code is
0
1
2
3
6
Examples:
(C1) quad_qaws(1/(x+1+2^(-4)), x, -1, 1, -0.5, -0.5, 1); (D1) [8.750097361672832, 1.24321522715422e-10, 170, 0] (C2) 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) (D2) ------------------------------- ALPHA 4 2 + 2 (C3) ev(%,alpha=4,numer); (D3) 8.750097361672829
Go to the first, previous, next, last section, table of contents.