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

# 23. Numerical

0

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

## 23.1 Introduction to fast Fourier transform

The `fft` package comprises functions for the numerical (not symbolic) computation of the fast Fourier transform.

·

@ref{Category: Fourier transform} · @ref{Category: Numerical methods} · @ref{Category: Share packages} · @ref{Category: Package fft}

0

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

## 23.2 Functions and Variables for fast Fourier transform

Function: polartorect (r, t)

  Translates complex values of the form `r %e^(%i t)` to the form `a + b %i`, where r is the magnitude and t is the phase. r and t are 1-dimensional arrays of the same size. The array size need not be a power of 2.

  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)

  `polartorect` is the inverse function of `recttopolar`.

  `load(fft)` loads this function. See also `fft`.

  ·

  @ref{Category: Package fft} · @ref{Category: Complex variables}

Function: recttopolar (a, b)

  Translates complex values of the form `a + b %i` to the form `r %e^(%i t)`, where a is the real part and b is the imaginary part. a and b are 1-dimensional arrays of the same size. The array size need not be a power of 2.

  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`.

  `recttopolar` is the inverse function of `polartorect`.

  `load(fft)` loads this function. See also `fft`.

  ·

  @ref{Category: Package fft} · @ref{Category: Complex variables}

Function: inverse_fft (y)

  Computes the inverse complex fast Fourier transform. y is a list or array (named or unnamed) which contains the data to transform. The number of elements must be a power of 2. The elements must be literal numbers (integers, rationals, floats, or bigfloats) or symbolic constants, or expressions `a + b*%i` where `a` and `b` are literal numbers or symbolic constants.

  `inverse_fft` returns a new object of the same type as y, which is not modified. Results are always computed as floats or expressions `a + b*%i` where `a` and `b` are floats.

  The inverse discrete Fourier transform is defined as follows. Let `x` be the output of the inverse transform. Then for `j` from 0 through `n - 1`,

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

  `load(fft)` loads this function.

  See also `fft` (forward transform), `recttopolar`, and `polartorect`.

  Examples:

  Real data.

      (%i1) load (fft) $
      (%i2) fpprintprec : 4 $
      (%i3) L : [1, 2, 3, 4, -1, -2, -3, -4] $
      (%i4) L1 : inverse_fft (L);
      (%o4) [0.0, 14.49 %i - .8284, 0.0, 2.485 %i + 4.828, 0.0,
                             4.828 - 2.485 %i, 0.0, - 14.49 %i - .8284]
      (%i5) L2 : fft (L1);
      (%o5) [1.0, 2.0 - 2.168L-19 %i, 3.0 - 7.525L-20 %i,
      4.0 - 4.256L-19 %i, - 1.0, 2.168L-19 %i - 2.0,
      7.525L-20 %i - 3.0, 4.256L-19 %i - 4.0]
      (%i6) lmax (abs (L2 - L));
      (%o6)                       3.545L-16

  Complex data.

      (%i1) load (fft) $
      (%i2) fpprintprec : 4 $
      (%i3) L : [1, 1 + %i, 1 - %i, -1, -1, 1 - %i, 1 + %i, 1] $
      (%i4) L1 : inverse_fft (L);
      (%o4) [4.0, 2.711L-19 %i + 4.0, 2.0 %i - 2.0,
      - 2.828 %i - 2.828, 0.0, 5.421L-20 %i + 4.0, - 2.0 %i - 2.0,
      2.828 %i + 2.828]
      (%i5) L2 : fft (L1);
      (%o5) [4.066E-20 %i + 1.0, 1.0 %i + 1.0, 1.0 - 1.0 %i,
      1.55L-19 %i - 1.0, - 4.066E-20 %i - 1.0, 1.0 - 1.0 %i,
      1.0 %i + 1.0, 1.0 - 7.368L-20 %i]
      (%i6) lmax (abs (L2 - L));
      (%o6)                       6.841L-17

  ·

  @ref{Category: Package fft}

Function: fft (x)

  Computes the complex fast Fourier transform. x is a list or array (named or unnamed) which contains the data to transform. The number of elements must be a power of 2. The elements must be literal numbers (integers, rationals, floats, or bigfloats) or symbolic constants, or expressions `a + b*%i` where `a` and `b` are literal numbers or symbolic constants.

  `fft` returns a new object of the same type as x, which is not modified. Results are always computed as floats or expressions `a + b*%i` where `a` and `b` are floats.

  The discrete Fourier transform is defined as follows. Let `y` be the output of the transform. Then for `k` from 0 through `n - 1`,

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

  When the data x are real, real coefficients `a` and `b` can be computed such that

      x[j] = sum (a[k] * cos (2*%pi*j*k / n) + b[k] * sin (2*%pi*j*k / n), k, 0, n/2)

  with

      a[0] = realpart (y[0])
      b[0] = 0

  and, for k from 1 through n/2 - 1,

      a[k] = realpart (y[k] + y[n - k])
      b[k] = imagpart (y[n - k] - y[k])

  and

      a[n/2] = realpart (y[n/2])
      b[n/2] = 0

  `load(fft)` loads this function.

  See also `inverse_fft` (inverse transform), `recttopolar`, and `polartorect`.

  Examples:

  Real data.

      (%i1) load (fft) $
      (%i2) fpprintprec : 4 $
      (%i3) L : [1, 2, 3, 4, -1, -2, -3, -4] $
      (%i4) L1 : fft (L);
      (%o4) [0.0, - 1.811 %i - .1036, 0.0, .6036 - .3107 %i, 0.0,
                               .3107 %i + .6036, 0.0, 1.811 %i - .1036]
      (%i5) L2 : inverse_fft (L1);
      (%o5) [1.0, 2.168L-19 %i + 2.0, 7.525L-20 %i + 3.0,
      4.256L-19 %i + 4.0, - 1.0, - 2.168L-19 %i - 2.0,
      - 7.525L-20 %i - 3.0, - 4.256L-19 %i - 4.0]
      (%i6) lmax (abs (L2 - L));
      (%o6)                       3.545L-16

  Complex data.

      (%i1) load (fft) $
      (%i2) fpprintprec : 4 $
      (%i3) L : [1, 1 + %i, 1 - %i, -1, -1, 1 - %i, 1 + %i, 1] $
      (%i4) L1 : fft (L);
      (%o4) [0.5, .3536 %i + .3536, - 0.25 %i - 0.25,
      0.5 - 6.776L-21 %i, 0.0, - .3536 %i - .3536, 0.25 %i - 0.25,
      0.5 - 3.388L-20 %i]
      (%i5) L2 : inverse_fft (L1);
      (%o5) [1.0 - 4.066E-20 %i, 1.0 %i + 1.0, 1.0 - 1.0 %i,
      - 1.008L-19 %i - 1.0, 4.066E-20 %i - 1.0, 1.0 - 1.0 %i,
      1.0 %i + 1.0, 1.947L-20 %i + 1.0]
      (%i6) lmax (abs (L2 - L));
      (%o6)                       6.83L-17

  Computation of sine and cosine coefficients.

      (%i1) load (fft) $
      (%i2) fpprintprec : 4 $
      (%i3) L : [1, 2, 3, 4, 5, 6, 7, 8] $
      (%i4) n : length (L) $
      (%i5) x : make_array (any, n) $
      (%i6) fillarray (x, L) $
      (%i7) y : fft (x) $
      (%i8) a : make_array (any, n/2 + 1) $
      (%i9) b : make_array (any, n/2 + 1) $
      (%i10) a[0] : realpart (y[0]) $
      (%i11) b[0] : 0 $
      (%i12) for k : 1 thru n/2 - 1 do
         (a[k] : realpart (y[k] + y[n - k]),
          b[k] : imagpart (y[n - k] - y[k]));
      (%o12)                        done
      (%i13) a[n/2] : y[n/2] $
      (%i14) b[n/2] : 0 $
      (%i15) listarray (a);
      (%o15)          [4.5, - 1.0, - 1.0, - 1.0, - 0.5]
      (%i16) listarray (b);
      (%o16)           [0, - 2.414, - 1.0, - .4142, 0]
      (%i17) f(j) := sum (a[k] * cos (2*%pi*j*k / n) + b[k] * sin (2*%pi*j*k / n), k, 0, n/2) $
      (%i18) makelist (float (f (j)), j, 0, n - 1);
      (%o18)      [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0]

  ·

  @ref{Category: Package fft}

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.

  ·

  @ref{Category: Translation and compilation}

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

  ·

  @ref{Category: Translation and compilation}

Option variable: fortspaces

  Default value: `false`

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

  ·

  @ref{Category: Translation and compilation}

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

  ·

  @ref{Category: Numerical methods}

Function: find_root (expr, x, a, b, [abserr, relerr])
Function: find_root (f, a, b, [abserr, relerr])
Function: bf_find_root (expr, x, a, b, [abserr, relerr])
Function: bf_find_root (f, a, b, [abserr, relerr])
Option variable: find_root_error
Option variable: find_root_abs
Option variable: find_root_rel

  Finds a root of the expression expr or the function f over the closed interval [a, b]. The expression expr may be an equation, in which case `find_root` seeks a root of `lhs(expr) - rhs(expr)`.

  Given that Maxima can evaluate expr or f over [a, b] and that expr or f is continuous, `find_root` is guaranteed to find the root, or one of the roots if there is more than one.

  `find_root` initially applies binary search. If the function in question appears to be smooth enough, `find_root` applies linear interpolation instead.

  `bf_find_root` is a bigfloat version of `find_root`. The function is computed using bigfloat arithmetic and a bigfloat result is returned. Otherwise, `bf_find_root` is identical to `find_root`, and the following description is equally applicable to `bf_find_root`.

  The accuracy of `find_root` is governed by `abserr` and `relerr`, which are optional keyword arguments to `find_root`. These keyword arguments take the form `key=val`. The keyword arguments are

  abserr

    Desired absolute error of function value at root. Default is `find_root_abs`.

  relerr

    Desired relative error of root. Default is `find_root_rel`.

  `find_root` stops when the function in question evaluates to something less than or equal to `abserr`, or if successive approximants x_0, x_1 differ by no more than `relerr * max(abs(x_0), abs(x_1))`. The default values of `find_root_abs` and `find_root_rel` are both zero.

  `find_root` expects the function in question to have a different sign at the endpoints of the search interval. When the function evaluates to a number at both endpoints and these numbers have the same sign, the behavior of `find_root` is governed by `find_root_error`. When `find_root_error` is `true`, `find_root` prints an error message. Otherwise `find_root` returns the value of `find_root_error`. The default value of `find_root_error` is `true`.

  If f evaluates to something other than a number at any step in the search algorithm, `find_root` returns a partially-evaluated `find_root` expression.

  The order of a and b is ignored; the region in which a root is sought is [min(a, b), max(a, b)].

  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
      (%i6) find_root (exp(x) = y, x, 0, 100);
                                  x
      (%o6)           find_root(%e  = y, x, 0.0, 100.0)
      (%i7) find_root (exp(x) = y, x, 0, 100), y = 10;
      (%o7)                   2.302585092994046
      (%i8) log (10.0);
      (%o8)                   2.302585092994046
      (%i9) fpprec:32;
      (%o9)                           32
      (%i10) bf_find_root (exp(x) = y, x, 0, 100), y = 10;
      (%o10)                  2.3025850929940456840179914546844b0
      (%i11) log(10b0);
      (%o11)                  2.3025850929940456840179914546844b0

  ·

  @ref{Category: Algebraic equations} · @ref{Category: Numerical methods}

Function: newton (expr, x, x_0, eps)

  Returns an approximate solution of `expr = 0` by Newton's method, considering expr to be a function of one variable, x. The search begins with `x = x_0` and proceeds until `abs(expr) < eps` (with expr evaluated at the current value of x).

  `newton` allows undefined variables to appear in expr, so long as the termination test `abs(expr) < eps` evaluates to `true` or `false`. Thus it is not necessary that expr evaluate to a number.

  `load(newton1)` loads this function.

  See also `realroots`, `allroots`, `find_root`, and `mnewton`.

  Examples:

      (%i1) load (newton1);
      (%o1) /usr/share/maxima/5.10.0cvs/share/numeric/newton1.mac
      (%i2) newton (cos (u), u, 1, 1/100);
      (%o2)                   1.570675277161251
      (%i3) ev (cos (u), u = %);
      (%o3)                 1.2104963335033528E-4
      (%i4) assume (a > 0);
      (%o4)                        [a > 0]
      (%i5) newton (x^2 - a^2, x, a/2, a^2/100);
      (%o5)                  1.00030487804878 a
      (%i6) ev (x^2 - a^2, x = %);
                                                 2
      (%o6)                6.098490481853958E-4 a

  ·

  @ref{Category: Algebraic equations} · @ref{Category: Numerical methods}

0

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

## 23.3 Introduction to Fourier series

The `fourie` package comprises functions for the symbolic computation of Fourier series. There are functions in the `fourie` package to calculate Fourier integral coefficients and some functions for manipulation of expressions.

·

@ref{Category: Fourier transform} · @ref{Category: Share packages} · @ref{Category: Package fourie}

0

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

## 23.4 Functions and Variables 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).

  ·

  @ref{Category: Package fourie}

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.

  ·

  @ref{Category: Package fourie}

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.

  ·

  @ref{Category: Package fourie}

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.

  ·

  @ref{Category: Package fourie} · @ref{Category: Integral calculus}

Function: fourier (f, x, p)

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

  ·

  @ref{Category: Package fourie}

Function: foursimp (l)

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

  ·

  @ref{Category: Package fourie} · @ref{Category: Trigonometric functions} · @ref{Category: Simplification functions}

Option variable: sinnpiflag

  Default value: `true`

  See `foursimp`.

  ·

  @ref{Category: Package fourie}

Option variable: cosnpiflag

  Default value: `true`

  See `foursimp`.

  ·

  @ref{Category: Package fourie}

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`.

  ·

  @ref{Category: Package fourie}

Function: fourcos (f, x, p)

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

  ·

  @ref{Category: Package fourie}

Function: foursin (f, x, p)

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

  ·

  @ref{Category: Package fourie}

Function: totalfourier (f, x, p)

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

  ·

  @ref{Category: Package fourie}

Function: fourint (f, x)

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

  ·

  @ref{Category: Package fourie}

Function: fourintcos (f, x)

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

  ·

  @ref{Category: Package fourie}

Function: fourintsin (f, x)

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

  ·

  @ref{Category: Package fourie}

0

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

This document was generated by Charlie & on January, 31 2011 using texi2html 1.76.