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

# 22. Numerical

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

## 22.1 Introduction to fast Fourier transform

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

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

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

Categories:  Package fft · 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`.

Categories:  Package fft · 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
```

Categories:  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 = realpart (y)
b = 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 : realpart (y) \$
(%i11) b : 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]
```

Categories:  Package fft

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

Categories:  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
```

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] [ ? ]

## 22.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.

Categories:  Fourier transform · Share packages · Package fourie

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

## 22.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).

Categories:  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.

Categories:  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.

Categories:  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.

Categories:  Package fourie · Integral calculus

Function: fourier (f, x, p)

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

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

Option variable: sinnpiflag

Default value: `true`

See `foursimp`.

Categories:  Package fourie

Option variable: cosnpiflag

Default value: `true`

See `foursimp`.

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

Categories:  Package fourie

Function: fourcos (f, x, p)

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

Categories:  Package fourie

Function: foursin (f, x, p)

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

Categories:  Package fourie

Function: totalfourier (f, x, p)

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

Categories:  Package fourie

Function: fourint (f, x)

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

Categories:  Package fourie

Function: fourintcos (f, x)

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

Categories:  Package fourie

Function: fourintsin (f, x)

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

Categories:  Package fourie

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

This document was generated by Oliver Kullmann on May, 18 2013 using texi2html 1.76.