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).
makes the change of variable given by
f(x,y) = 0 in all integrals occurring in exp with integration with
respect to x; y is the new variable.
(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
Function:DBLINT('F,'R,'S,a,b)
a double-integral routine which was written in
top-level macsyma 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(X,Y) must be a translated or compiled function of two
variables, and R(X) and S(X) 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! Ask
LPH (or GJC) about using these numerical aids. The file SHARE1;DBLINT
DEMO can be run in batch or demo mode to illustrate the usage on a
sample problem; the file SHARE1;DBLNT DEMO1 is an extension of the DEMO
which also makes use of other numerical aids, FLOATDEFUNK and QUANC8.
Please send all bug notes and questions to LPH
Function:DEFINT(exp, var, low, high)
DEFinite INTegration, the same as
INTEGRATE(exp,var,low,high). This uses symbolic methods, if
you wish to use a numerical method try ROMBERG(exp,var,low,high).
Function:ERF(X)
the error function, whose derivative is:
2*EXP(-X^2)/SQRT(%PI).
Variable:ERFFLAG
default: [TRUE] if FALSE prevents RISCH from introducing the
ERF function in the answer if there were none in the integrand to
begin with.
Variable:ERRINTSCE
default: [TRUE] - If a call to the INTSCE routine is not of
the form
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.
Function:ILT(exp, lvar, ovar)
takes the inverse Laplace transform of exp with
respect to lvar and parameter ovar. exp 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.
(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
Function:INTEGRATE(exp, var)
integrates exp with respect to var or returns an
integral expression (the noun form) if it cannot perform the
integration (see note 1 below). Roughly speaking three stages are
used:
(1) INTEGRATE sees if the integrand is of the form
F(G(X))*DIFF(G(X),X) by testing whether the derivative of some
subexpression (i.e. G(X) in the above case) divides the integrand. If
so it looks up F in a table of integrals and substitutes G(X) for X in
the integral of F. This may make use of gradients in taking the
derivative. (If an unknown function appears in the integrand it must
be eliminated in this stage or else INTEGRATE will return the noun
form of the integrand.)
(2) INTEGRATE tries to match the integrand to a form for which a
specific method can be used, e.g. trigonometric substitutions.
(3) If the first two stages fail it uses the Risch algorithm.
Functional relationships must be explicitly represented in order
for INTEGRATE to work properly. INTEGRATE is not affected by
DEPENDENCIES set up with the DEPENDS command.
INTEGRATE(exp, var, low, high) finds the definite integral of exp with
respect to var from low to high or returns the noun form if it cannot
perform the integration. The limits should not contain var. Several
methods are used, including direct substitution in the indefinite
integral and contour integration. Improper integrals may use the
names INF for positive infinity and MINF for negative infinity. If an
integral "form" is desired for manipulation (for example, an integral
which cannot be computed until some numbers are substituted for some
parameters), the noun form 'INTEGRATE may be used and this will
display with an integral sign. (See Note 1 below.)
The function LDEFINT uses LIMIT to evaluate the integral at the
lower and upper limits.
Sometimes during integration the user may be asked what the sign
of an expression is. Suitable responses are POS;, ZERO;, or NEG;.
(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
Variable:INTEGRATION_CONSTANT_COUNTER
- a counter which is updated each time a
constant of integration (called by MACSYMA, e.g., "INTEGRATIONCONSTANT1")
is introduced into an expression by indefinite integration of an equation.
Variable:INTEGRATE_USE_ROOTSOF
default: [false] If not false then when the denominator of
an rational function cannot be factored, we give the integral
in a form which is a sum over the roots of the denominator:
(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.
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.
Function:INTSCE(expr,var)
INTSCE LISP contains a routine, written by Richard
Bogen, for integrating products of sines,cosines and exponentials of
the form
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.
Function:LDEFINT(exp,var,ll,ul)
yields the definite integral of exp by using
LIMIT to evaluate the indefinite integral of exp with respect to var
at the upper limit ul and at the lower limit ll.
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 file SHARE1;QQ FASL (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 name,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.
Do PRINTFILE(QQ,USAGE,SHARE1) for details.
Function:QUANC8('function name,lo,hi)
An adaptive integrator, available in
SHARE1;QQ FASL. DEMO 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"); For more details do
DESCRIBE(QQ); .
Function:RESIDUE(exp, var, val)
computes the residue in the complex plane of
the expression exp when the variable var assumes the value val. The
residue is the coefficient of (var-val)**(-1) in the Laurent series
for exp.
integrates exp with respect to var 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[TRUE] - if FALSE prevents RISCH from introducing the ERF
function in the answer if there were none in the integrand to begin
with.
(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)
Function:ROMBERG(exp,var,ll,ul)
or ROMBERG(exp,ll,ul) - Romberg Integration.
You need not load in any file to use ROMBERG, it is autoloading.
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>);
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:
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.)
Variable:ROMBERGABS
default: [0.0] (0.0B0) 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
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.
Variable:ROMBERGIT
default: [11] - The accuracy of the ROMBERG integration
command is governed by the global variables ROMBERGTOL[1.E-4] and
ROMBERGIT[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.
Variable:ROMBERGMIN
default: [0] - 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.
Variable:ROMBERGTOL
default: [1.E-4] - The accuracy of the ROMBERG integration
command is governed by the global variables ROMBERGTOL[1.E-4] and
ROMBERGIT[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.