[ < ] [ > ]   [ << ] [ Up ] [ >> ]         [Top] [Contents] [Index] [ ? ]

23. Numerical


[ < ] [ > ]   [ << ] [ Up ] [ >> ]         [Top] [Contents] [Index] [ ? ]

23.1 Introduction to Numerical


[ < ] [ > ]   [ << ] [ Up ] [ >> ]         [Top] [Contents] [Index] [ ? ]

23.2 Fourier packages

The fft package comprises functions for the numerical (not symbolic) computation of the fast Fourier transform. load ("fft") loads this package. See fft.

The fourie package comprises functions for the symbolic computation of Fourier series. load ("fourie") loads this package. There are functions in the fourie package to calculate Fourier integral coefficients and some functions for manipulation of expressions. See Definitions for Fourier Series.


[ < ] [ > ]   [ << ] [ Up ] [ >> ]         [Top] [Contents] [Index] [ ? ]

23.3 Definitions for Numerical

Function: polartorect (magnitude_array, phase_array)

Translates complex values of the form r %e^(%i t) to the form a + b %i. load ("fft") loads this function into Maxima. See also fft.

The magnitude and phase, r and t, are taken from magnitude_array and phase_array, respectively. The original values of the input arrays are replaced by the real and imaginary parts, a and b, on return. The outputs are calculated as

 
a: r cos (t)
b: r sin (t)

The input arrays must be the same size and 1-dimensional. The array size need not be a power of 2.

polartorect is the inverse function of recttopolar.

Function: recttopolar (real_array, imaginary_array)

Translates complex values of the form a + b %i to the form r %e^(%i t). load ("fft") loads this function into Maxima. See also fft.

The real and imaginary parts, a and b, are taken from real_array and imaginary_array, respectively. The original values of the input arrays are replaced by the magnitude and angle, r and t, on return. The outputs are calculated as

 
r: sqrt (a^2 + b^2)
t: atan2 (b, a)

The computed angle is in the range -%pi to %pi.

The input arrays must be the same size and 1-dimensional. The array size need not be a power of 2.

recttopolar is the inverse function of polartorect.

Function: ift (real_array, imaginary_array)

Fast inverse discrete Fourier transform. load ("fft") loads this function into Maxima.

ift carries out the inverse complex fast Fourier transform on 1-dimensional floating point arrays. The inverse transform is defined as

 
x[j]: sum (y[j] exp (+2 %i %pi j k / n), k, 0, n-1)

See fft for more details.

Function: fft (real_array, imaginary_array)
Function: ift (real_array, imaginary_array)
Function: recttopolar (real_array, imaginary_array)
Function: polartorect (magnitude_array, phase_array)

Fast Fourier transform and related functions. load ("fft") loads these functions into Maxima.

fft and ift carry out the complex fast Fourier transform and inverse transform, respectively, on 1-dimensional floating point arrays. The size of imaginary_array must equal the size of real_array.

fft and ift operate in-place. That is, on return from fft or ift, the original content of the input arrays is replaced by the output. The fillarray function can make a copy of an array, should it be necessary.

The discrete Fourier transform and inverse transform are defined as follows. Let x be the original data, with

 
x[i]: real_array[i] + %i imaginary_array[i]

Let y be the transformed data. The forward and inverse transforms are

 
y[k]: (1/n) sum (x[j] exp (-2 %i %pi j k / n), j, 0, n-1)

x[j]:       sum (y[j] exp (+2 %i %pi j k / n), k, 0, n-1)

Suitable arrays can be allocated by the array function. For example:

 
array (my_array, float, n-1)$

declares a 1-dimensional array with n elements, indexed from 0 through n-1 inclusive. The number of elements n must be equal to 2^m for some m.

fft can be applied to real data (imaginary array all zeros) to obtain sine and cosine coefficients. After calling fft, the sine and cosine coefficients, say a and b, can be calculated as

 
a[0]: real_array[0]
b[0]: 0

and

 
a[j]: real_array[j] + real_array[n-j]
b[j]: imaginary_array[j] - imaginary_array[n-j]

for j equal to 1 through n/2-1, and

 
a[n/2]: real_array[n/2]
b[n/2]: 0

recttopolar translates complex values of the form a + b %i to the form r %e^(%i t). See recttopolar.

polartorect translates complex values of the form r %e^(%i t) to the form a + b %i. See polartorect.

demo ("fft") displays a demonstration of the fft package.

Option variable: fortindent

Default value: 0

fortindent 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: fortran (expr)

Prints expr as a Fortran statement. The output line is indented with spaces. If the line is too long, fortran prints continuation lines. fortran prints the exponentiation operator ^ as **, and prints a complex number a + b %i in the form (a,b).

expr may be an equation. If so, fortran prints an assignment statement, assigning the right-hand side of the equation to the left-hand side. In particular, if the right-hand side of expr is the name of a matrix, then fortran prints an assignment statement for each element of the matrix.

If expr is not something recognized by fortran, the expression is printed in grind format without complaint. fortran does not know about lists, arrays, or functions.

fortindent controls the left margin of the printed lines. 0 is the normal margin (i.e., indented 6 spaces). Increasing fortindent causes expressions to be printed further to the right.

When fortspaces is true, fortran fills out each printed line with spaces to 80 columns.

fortran evaluates its arguments; quoting an argument defeats evaluation. fortran always returns done.

Examples:

 
(%i1) expr: (a + b)^12$
(%i2) fortran (expr);
      (b+a)**12                                                                 
(%o2)                         done
(%i3) fortran ('x=expr);
      x = (b+a)**12                                                             
(%o3)                         done
(%i4) fortran ('x=expand (expr));
      x = b**12+12*a*b**11+66*a**2*b**10+220*a**3*b**9+495*a**4*b**8+792        
     1   *a**5*b**7+924*a**6*b**6+792*a**7*b**5+495*a**8*b**4+220*a**9*b        
     2   **3+66*a**10*b**2+12*a**11*b+a**12                                     
(%o4)                         done
(%i5) fortran ('x=7+5*%i);
      x = (7,5)                                                                 
(%o5)                         done
(%i6) fortran ('x=[1,2,3,4]);
      x = [1,2,3,4]                                                             
(%o6)                         done
(%i7) f(x) := x^2$
(%i8) fortran (f);
      f                                                                         
(%o8)                         done
Option variable: fortspaces

Default value: false

When fortspaces is true, fortran fills out each printed line with spaces to 80 columns.

Function: horner (expr, x)
Function: horner (expr)

Returns a rearranged representation of expr as in Horner's rule, using x as the main variable if it is specified. x may be omitted in which case the main variable of the canonical rational expression form of expr is used.

horner sometimes improves stability if expr is to be numerically evaluated. It is also useful if Maxima is used to generate programs to be run in Fortran. See also stringout.

 
(%i1) expr: 1e-155*x^2 - 5.5*x + 5.2e155;
                           2
(%o1)            1.0E-155 x  - 5.5 x + 5.2E+155
(%i2) expr2: horner (%, x), keepfloat: true;
(%o2)            (1.0E-155 x - 5.5) x + 5.2E+155
(%i3) ev (expr, x=1e155);
Maxima encountered a Lisp error:

 floating point overflow

Automatically continuing.
To reenable the Lisp debugger set *debugger-hook* to nil.
(%i4) ev (expr2, x=1e155);
(%o4)                       7.0E+154
Function: find_root (f(x), x, a, b)
Function: find_root (f, a, b)

Finds the zero of function f as variable x varies over the range [a, b]. The function must have a different sign at each endpoint. If this condition is not met, the action of the function is governed by find_root_error. If find_root_error is true then an error occurs, otherwise the value of find_root_error is returned (thus for plotting find_root_error might be set to 0.0). Otherwise (given that Maxima can evaluate the first argument in the specified range, and that it is continuous) find_root is guaranteed to come up with the zero (or one of them if there is more than one zero). The accuracy of find_root is governed by find_root_abs and find_root_rel which must be non-negative floating point numbers. find_root will stop when the first arg evaluates to something less than or equal to find_root_abs or if successive approximants to the root differ by no more than find_root_rel * <one of the approximants>. The default values of find_root_abs and find_root_rel are 0.0 so find_root 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

 
find_root (sin(x) = x/2, x, %pi, 0.1);

is equivalent to

 
find_root (sin(x) = x/2, x, 0.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.

Examples:

 
(%i1) f(x) := sin(x) - x/2;
                                        x
(%o1)                  f(x) := sin(x) - -
                                        2
(%i2) find_root (sin(x) - x/2, x, 0.1, %pi);
(%o2)                   1.895494267033981
(%i3) find_root (sin(x) = x/2, x, 0.1, %pi);
(%o3)                   1.895494267033981
(%i4) find_root (f(x), x, 0.1, %pi);
(%o4)                   1.895494267033981
(%i5) find_root (f, 0.1, %pi);
(%o5)                   1.895494267033981
Option variable: find_root_abs

Default value: 0.0

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

Option variable: find_root_error

Default value: true

find_root_error governs the behavior of find_root. When find_root is called, it determines whether or not the function to be solved 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 find_root_error is true, then an error is signaled. If they are of like sign and find_root_error is not true, the value of find_root_error is returned. Thus for plotting, find_root_error might be set to 0.0.

Option variable: find_root_rel

Default value: 0.0

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


[ < ] [ > ]   [ << ] [ Up ] [ >> ]         [Top] [Contents] [Index] [ ? ]

23.4 Definitions for Fourier Series

Function: equalp (x, y)

Returns true if equal (x, y) otherwise false (doesn't give an error message like equal (x, y) would do in this case).

Function: remfun (f, expr)
Function: remfun (f, expr, x)

remfun (f, expr) replaces all occurrences of f (arg) by arg in expr.

remfun (f, expr, x) replaces all occurrences of f (arg) by arg in expr only if arg contains the variable x.

Function: funp (f, expr)
Function: funp (f, expr, x)

funp (f, expr) returns true if expr contains the function f.

funp (f, expr, x) returns true if expr contains the function f and the variable x is somewhere in the argument of one of the instances of f.

Function: absint (f, x, halfplane)
Function: absint (f, x)
Function: absint (f, x, a, b)

absint (f, x, halfplane) returns the indefinite integral of f with respect to x in the given halfplane (pos, neg, or both). f may contain expressions of the form abs (x), abs (sin (x)), abs (a) * exp (-abs (b) * abs (x)).

absint (f, x) is equivalent to absint (f, x, pos).

absint (f, x, a, b) returns the definite integral of f with respect to x from a to b. f may include absolute values.

Function: fourier (f, x, p)

Returns a list of the Fourier coefficients of f(x) defined on the interval [-%pi, %pi].

Function: foursimp (l)

Simplifies sin (n %pi) to 0 if sinnpiflag is true and cos (n %pi) to (-1)^n if cosnpiflag is true.

Option variable: sinnpiflag

Default value: true

See foursimp.

Option variable: cosnpiflag

Default value: true

See foursimp.

Function: fourexpand (l, x, p, limit)

Constructs and returns the Fourier series from the list of Fourier coefficients l up through limit terms (limit may be inf). x and p have same meaning as in fourier.

Function: fourcos (f, x, p)

Returns the Fourier cosine coefficients for f(x) defined on [0, %pi].

Function: foursin (f, x, p)

Returns the Fourier sine coefficients for f(x) defined on [0, %pi].

Function: totalfourier (f, x, p)

Returns fourexpand (foursimp (fourier (f, x, p)), x, p, 'inf).

Function: fourint (f, x)

Constructs and returns a list of the Fourier integral coefficients of f(x) defined on [minf, inf].

Function: fourintcos (f, x)

Returns the Fourier cosine integral coefficients for f(x) on [0, inf].

Function: fourintsin (f, x)

Returns the Fourier sine integral coefficients for f(x) on [0, inf].


[ << ] [ >> ]           [Top] [Contents] [Index] [ ? ]

This document was generated on March, 19 2006 using texi2html 1.76.