[ < ] | [ > ] |  | [ << ] | [ 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 `Functions and Variables for Fourier Series`.

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

## 23.3 Functions and Variables 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 (expr, x, a, b)
Function: find_root (f, a, b)
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.

  The accuracy of `find_root` is governed by `find_root_abs` and `find_root_rel`. `find_root` stops when the function in question evaluates to something less than or equal to `find_root_abs`, or if successive approximants x_0, x_1 differ by no more than `find_root_rel * 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. If this condition is not met, 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

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

[ < ] | [ > ] |  | [ << ] | [ 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).

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 `[-p, p]`.

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, p]`.

Function: foursin (f, x, p)

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

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 by Robert Dodier on August, 25 2007 using texi2html 1.76.