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, and the IMSL version of Romberg, DCADRE. There is also an adaptave integrator which uses the Newton-Cotes 8 panel quadrature rule, called QUANC8. 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).
(C1) 'INTEGRATE(%E**SQRT(A*Y),Y,0,4);
4
/
[ SQRT(A) SQRT(Y)
(D1) I (%E ) dY
]
/
0
(C2) CHANGEVAR(D1,Y-Z^2/A,Z,Y);
2 SQRT(A)
/
[ Z
2 I Z %E dZ
]
/
0
(D4) ---------------------
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.
(C3) SUM(A[I]*X^(I-2),I,0,INF);
INF
====
\ I - 2
(D3) > A X
/ I
====
I = 0
(C4) CHANGEVAR(%,I-2-N,N,I);
INF
====
\ N
(D4) > 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.
(C1) 'INTEGRATE(SINH(A*X)*F(T-X),X,0,T)+B*F(T)=T**2;
T
/
[ 2
(D1) I (SINH(A X) F(T - X)) dX + B F(T) = T
]
/
0
(C2) LAPLACE(%,T,S);
A LAPLACE(F(T), T, S)
(D2) ---------------------
2 2
S - A
2
+ B LAPLACE(F(T), T, S) = --
3
S
(C3) LINSOLVE([%],['LAPLACE(F(T),T,S)]);
SOLUTION
2 2
2 S - 2 A
(E3) LAPLACE(F(T), T, S) = --------------------
5 2 3
B S + (A - A B) S
(D3) [E3]
(C4) ILT(E3,S,T);
IS A B (A B - 1) POSITIVE, NEGATIVE, OR ZERO?
POS;
2
SQRT(A) SQRT(A B - B) T
2 COSH(------------------------)
B
(D4) F(T) = - --------------------------------
A
2
A T 2
+ ------- + ------------------
A B - 1 3 2 2
A B - 2 A B + A
(C1) INTEGRATE(SIN(X)**3,X);
3
COS (X)
(D1) ------- - COS(X)
3
(C2) 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
(D2) BETA(A + 1, - - A)
2
(C3) GRADEF(Q(X),SIN(X**2));
(D3) Q(X)
(C4) DIFF(LOG(Q(R(X))),X);
d 2
(-- R(X)) SIN(R (X))
dX
(D4) --------------------
Q(R(X))
(C5) INTEGRATE(%,X);
(D5) 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.
(C6) INTEGRATE(1/(X^3+X+1),X);
/
[ 1
(D6) I ---------- dX
] 3
/ X + X + 1
(C4) 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)
(D4) ----------------- - --------------- + ---------------
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.
(C5) INTEGRATE_USE_ROOTSOF:true; (D5) TRUE
(C6) 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)
(D6) ----------------------------------------------------------
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.
(C1) RESIDUE(S/(S**2+A**2),S,A*%I);
1
(D1) -
2
(C2) RESIDUE(SIN(A*X)/X**4,X,0);
3
A
(D2) - --
6
(C1) RISCH(X^2*ERF(X),X);
2 2
- X X 3 2
%E (%E SQRT(%PI) X ERF(X) + X + 1)
(D1) ------------------------------------------
3 SQRT(%PI)
(C2) DIFF(%,X),RATSIMP;
2
(D2) 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. (The IMSL version of Romberg integration is now available in Macsyma. Do DESCRIBE(DCADRE); for more information.)
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.
Go to the first, previous, next, last section, table of contents.