NDIFFQ Numerical

12.6: Definitions for Numerical

Function: FFT (real-array, imag-array)

Fast Fourier Transform. This package may be loaded by doing LOAD(FFT); There is also an IFT command, for Inverse Fourier Transform. These functions perform a (complex) fast fourier transform on either 1 or 2 dimensional FLOATING-POINT arrays, obtained by:

ARRAY(<ary>,FLOAT,<dim1>); or
ARRAY(<ary>,FLOAT,<dim1>,<dim2>);

For 1D arrays

<dim1> = 2^n-1

and for 2D arrays

<dim1>=<dim2>=2^n-1

(i.e. the array is square). (Recall that MACSYMA arrays are indexed from a 0 origin so that there will be 2^n and (2^n)^2 arrays elements in the above two cases.) This package also contains two other functions, POLARTORECT and RECTTOPOLAR. Do DESCRIBE(cmd) for details. For details on the implementation, do PRINTFILE(FFT,USAGE,SHARE); .

Variable: FORTINDENT

default: [0] - controls the left margin indentation of expressions printed out by the FORTRAN command. 0 gives normal printout (i.e. 6 spaces), and positive values will causes the expressions to be printed farther to the right.

Function: FORTMX (name,matrix)

converts a MACSYMA matrix into a sequence of FORTRAN assignment statements of the form name(i,j)=<corresponding matrix element>. This command is now obsolete. FORTMX(name,matrix); may now be done as FORTRAN(name=matrix);. (If "name" is bound, FORTRAN('name=matrix); may be necessary.) Please convert code that uses the FORTMX command as it may be flushed some day.

Function: FORTRAN (exp)

converts exp into a FORTRAN linear expression in legal FORTRAN with 6 spaces inserted at the beginning of each line, continuation lines, and ** rather than ^ for exponentiation. When the option FORTSPACES[FALSE] is TRUE, the FORTRAN command fills out to 80 columns using spaces. If FORTRAN is called on a bound symbolic atom, e.g. FORTRAN(X); where X:A*B$ has been done, then X={value of X}, e.g. X=A*B will be generated. In particular, if e.g. M:MATRIX(...); has been done, then FORTRAN(M); will generate the appropriate assignment statements of the form name(i,j)=<corresponding matrix element>. FORTINDENT[0] controls the left margin of expressions printed out, 0 is the normal margin (i.e. indented 6 spaces), increasing it will cause the expression to be printed further to the right.

Variable: FORTSPACES

default: [FALSE] - if TRUE, the FORTRAN command fills out to 80 columns using spaces.

Function: HORNER (exp, var)

will convert exp into a rearranged representation as in Horner's rule, using var as the main variable if it is specified. Var may also be omitted in which case the main variable of the CRE form of exp is used. HORNER sometimes improves stability if expr is to be numerically evaluated. It is also useful if MACSYMA is used to generate programs to be run in FORTRAN (see DESCRIBE(STRINGOUT);)

(C1) 1.0E-20*X^2-5.5*X+5.2E20;
                                2
(D1)                   1.0E-20 X  - 5.5 X + 5.2E+20
(C2) HORNER(%,X),KEEPFLOAT:TRUE;
(D2)                  X (1.0E-20 X - 5.5) + 5.2E+20
(C3) D1,X=1.0E20;
ARITHMETIC OVERFLOW
(C4) D2,X=1.0E20;
(D4)                          6.9999999E+19

Function: IFT (real-array, imag-array)

Inverse Fourier Transform. Do LOAD(FFT); to load in this package. These functions (FFT and IFT) perform a (complex) fast fourier transform on either 1 or 2 dimensional FLOATING-POINT arrays, obtained by: ARRAY(<ary>,FLOAT,<dim1>); or ARRAY(<ary>,FLOAT,<dim1>,<dim2>); For 1D arrays <dim1> must equal 2^n-1, and for 2D arrays <dim1>=<dim2>=2^n-1 (i.e. the array is square). (Recall that MACSYMA arrays are indexed from a 0 origin so that there will be 2^n and (2^n)^2 arrays elements in the above two cases.) For details on the implementation, do PRINTFILE(FFT,USAGE,SHARE); .

Function: INTERPOLATE (func,x,a,b)

finds the zero of func as x varies. The last two args give the range to look in. The function must have a different sign at each endpoint. If this condition is not met, the action of the of the function is governed by INTPOLERROR[TRUE]). If INTPOLERROR is TRUE then an error occurs, otherwise the value of INTPOLERROR is returned (thus for plotting INTPOLERROR might be set to 0.0). Otherwise (given that MACSYMA can evaluate the first argument in the specified range, and that it is continuous) INTERPOLATE is guaranteed to come up with the zero (or one of them if there is more than one zero). The accuracy of INTERPOLATE is governed by INTPOLABS[0.0] and INTPOLREL[0.0] which must be non-negative floating point numbers. INTERPOLATE will stop when the first arg evaluates to something less than or equal to INTPOLABS or if successive approximants to the root differ by no more than INTPOLREL * <one of the approximants>. The default values of INTPOLABS and INTPOLREL are 0.0 so INTERPOLATE gets as good an answer as is possible with the single precision arithmetic we have. The first arg may be an equation. The order of the last two args is irrelevant. Thus

INTERPOLATE(SIN(X)=X/2,X,%PI,.1);
   is equivalent to
INTERPOLATE(SIN(X)=X/2,X,.1,%PI);

The method used is a binary search in the range specified by the last two args. When it thinks the function is close enough to being linear, it starts using linear interpolation. An alternative syntax has been added to interpolate, this replaces the first two arguments by a function name. The function MUST be TRANSLATEd or compiled function of one argument. No checking of the result is done, so make sure the function returns a floating point number.

F(X):=(MODE_DECLARE(X,FLOAT),SIN(X)-X/2.0);
INTERPOLATE(SIN(X)-X/2,X,0.1,%PI)       time= 60 msec
INTERPOLATE(F(X),X,0.1,%PI);            time= 68 msec
TRANSLATE(F);
INTERPOLATE(F(X),X,0.1,%PI);            time= 26 msec
INTERPOLATE(F,0.1,%PI);                 time=  5 msec

There is also a Newton method interpolation routine, do DESCRIBE(NEWTON); .

Variable: INTPOLABS

default: [0.0] - The accuracy of the INTERPOLATE command is governed by INTPOLABS[0.0] and INTPOLREL[0.0] which must be non-negative floating point numbers. INTERPOLATE will stop when the first arg evaluates to something less than or equal to INTPOLABS or if successive approximants to the root differ by no more than INTPOLREL * <one of the approximants>. The default values of INTPOLABS and INTPOLREL are 0.0 so INTERPOLATE gets as good an answer as is possible with the single precision arithmetic we have.

Variable: INTPOLERROR

default: [TRUE] - Governs the behavior of INTERPOLATE. When INTERPOLATE is called, it determines whether or not the function to be interpolated satisfies the condition that the values of the function at the endpoints of the interpolation interval are opposite in sign. If they are of opposite sign, the interpolation proceeds. If they are of like sign, and INTPOLERROR is TRUE, then an error is signaled. If they are of like sign and INTPOLERROR is not TRUE, the value of INTPOLERROR is returned. Thus for plotting, INTPOLERROR might be set to 0.0.

Variable: INTPOLREL

default: [0.0] - The accuracy of the INTERPOLATE command is governed by INTPOLABS[0.0] and INTPOLREL[0.0] which must be non-negative floating point numbers. INTERPOLATE will stop when the first arg evaluates to something less than or equal to INTPOLABS or if successive approximants to the root differ by no more than INTPOLREL * <one of the approximants>. The default values of INTPOLABS and INTPOLREL are 0.0 so INTERPOLATE gets as good an answer as is possible with the single precision arithmetic we have.

Function: NEWTON (exp,var,X0,eps)

The file NEWTON 1 on the SHARE directory contains a function which will do interpolation using Newton's method. It may be accessed by LOAD(NEWTON); . The Newton method can do things that INTERPOLATE will refuse to handle, since INTERPOLATE requires that everything evaluate to a flonum. Thus NEWTON(x^2-a^2,x,a/2,a^2/100); will say that it can't tell if flonum*a^2<a^2/100. Doing ASSUME(a>0); and then doing NEWTON again works. You get x=a+<small flonum>*a which is symbolic all the way. INTERPOLATE(x^2-a^2,x,a/2,2*a); complains that .5*a is not flonum... An adaptive integrator which uses the Newton-Cotes 8 panel quadrature rule is available in SHARE1;QQ FASL. Do DESCRIBE(QQ) for details.

Function: POLARTORECT (magnitude-array, phase-array)

converts from magnitude and phase form into real and imaginary form putting the real part in the magnitude array and the imaginary part into the phase array

<real>=<magnitude>*COS(<phase>) ==>
  <imaginary>=<magnitude>*SIN(<phase>

This function is part of the FFT package. Do LOAD(FFT); to use it. Like FFT and IFT this function accepts 1 or 2 dimensional arrays. However, the array dimensions need not be a power of 2, nor need the 2D arrays be square.

Function: RECTTOPOLAR (real-array, imag-array)

undoes POLARTORECT. The phase is given in the range from -%PI to %PI. This function is part of the FFT package. Do LOAD(FFT); to use it. Like FFT and IFT this function accepts 1 or 2 dimensional arrays. However, the array dimensions need not be a power of 2, nor need the 2D arrays be square.

NDIFFQ

Numerical