This is Info file gsl-ref.info, produced by Makeinfo version 1.67 from the input file gsl-ref.texi. INFO-DIR-SECTION Scientific software START-INFO-DIR-ENTRY * gsl-ref: (gsl-ref). GNU Scientific Library - Reference END-INFO-DIR-ENTRY This file documents the *GNU Scientific Library*. Copyright (C) 1996, 1997, 1998 The GSL Project. Permission is granted to make and distribute verbatim copies of this manual provided the copyright notice and this permission notice are preserved on all copies. Permission is granted to copy and distribute modified versions of this manual under the conditions for verbatim copying, provided that the entire resulting derived work is distributed under the terms of a permission notice identical to this one. Permission is granted to copy and distribute translations of this manual into another language, under the above conditions for modified versions, except that this permission notice may be stated in a translation approved by the Foundation. File: gsl-ref.info, Node: Example of using mixed-radix FFT routines for complex data, Prev: Mixed-radix FFT routines for complex data, Up: FFTs of complex data Example of using mixed-radix FFT routines for complex data ---------------------------------------------------------- Here is an example program which computes the FFT of a short pulse in a sample of length 630 (=2*3*3*5*7) using the mixed-radix algorithm. #include #include #include #include int main () { int i; const int n = 630 ; complex data[n]; gsl_fft_complex_wavetable wavetable; for (i = 0; i < n; i++) { data[i].real = 0.0 ; data[i].imag = 0.0 ; } data[0].real = 1.0 ; for (i = 1; i <= 10; i++) { data[i].real = data[n-i].real = 1.0 ; } for (i = 0; i < n; i++) { printf ("%d: %e %e\n", i, data[i].real, data[i].imag); } printf ("\n"); gsl_fft_complex_wavetable_alloc (n, &wavetable); gsl_fft_complex_init (n, &wavetable); for (i = 0; i < wavetable.nf; i++) { printf("# factor %d: %d\n", i, wavetable.factor[i]); } gsl_fft_complex_forward (data, n, &wavetable); for (i = 0; i < n; i++) { printf ("%d: %e %e\n", i, data[i].real, data[i].imag); } gsl_fft_complex_wavetable_free (&wavetable); } Note that we have assumed that the program is using the default `gsl' error handler (which calls `abort' for any errors). If you are not using a safe error handler you would need to check the return status of of all the `gsl' routines. File: gsl-ref.info, Node: FFTs of real data, Next: Test Signals, Prev: FFTs of complex data, Up: FFTs FFTs of real data ================= The functions for real data are similar to those for complex data. However, there is an important difference between forward and inverse transforms. The fourier transform of a real sequence is not real. It is a complex sequence with a special symmetry: z_k = z_{N-k}^* A sequence with this symmetry is called "conjugate-complex" or "half-complex". This different structure requires different storage-layouts for the forward transform (from real to half-complex) and inverse transform (from half-complex back to real). As a consequence the routines are divided into two sets: functions in `gsl_fft_real' which operate on real sequences and functions in `gsl_fft_halfcomplex' which operate on half-complex sequences. Functions in `gsl_fft_real' compute the frequency coefficients of a real sequence. The half-complex coefficients c of a real sequence x are given by fourier analysis, c_k = \sum_{j=0}^{N-1} x_k \exp(-2 \pi i j k /N) Functions in `gsl_fft_halfcomplex' compute inverse or backwards transforms. They reconstruct real sequences by fourier synthesis from their half-complex frequency coefficients, c, x_j = {1 \over N} \sum_{k=0}^{N-1} c_k \exp(2 \pi i j k /N) The symmetry of the half-complex sequence implies that only half of the complex numbers in the output need to be stored. The remaining half can be reconstructed using the half-complex symmetry condition. (This works for all lengths, even and odd. When the length is even the middle value, where k=N/2, is also real). Thus only N real numbers are required to store the half-complex sequence, and the transform of a real sequence can be stored in the same size array as the original data. The precise storage arrangements depend on the algorithm, and are different for radix-2 and mixed-radix routines. The radix-2 function operates in-place, which constrain the locations where each element can be stored. The restriction forces real and imaginary parts to be stored far apart. The mixed-radix algorithm does not have this restriction, and it stores the real and imaginary parts of a given term in neighboring locations. This is desirable for better locality of memory accesses. * Menu: * Radix-2 FFT routines for real data:: * Mixed-radix FFT routines for real data:: * Example of using mixed-radix FFT routines for real data:: File: gsl-ref.info, Node: Radix-2 FFT routines for real data, Next: Mixed-radix FFT routines for real data, Up: FFTs of real data Radix-2 FFT routines for real data ---------------------------------- This section describes radix-2 FFT algorithms for real data. They use the Cooley-Tukey algorithm to compute in-place FFTs for lengths which are a power of 2. The radix-2 FFT functions for real data are declared in the header files `gsl_fft_real.h' - Function: int gsl_fft_real_radix2 (double data[], const unsigned int n) This function computes an in-place radix-2 FFT of DATA, a real array of length N. The output is a half-complex sequence, which is stored in-place. The arrangement of the half-complex terms uses the following scheme: for k < N/2 the real part of the k-th term is stored in location k, and the corresponding imaginary part is stored in location N-k. Terms with k > N/2 can be reconstructed using the symmetry z_k = z^*_{N-k}. The terms for k=0 and k=N/2 are both purely real, and count as a special case. Their real parts are stored in locations 0 and N/2 respectively, while their imaginary parts which are zero are not stored. The following table shows the correspondence between the output DATA and the equivalent results obtained by considering the input data as a complex sequence with zero imaginary part, complex[0].real = data[0] complex[0].imag = 0 complex[1].real = data[1] complex[1].imag = data[N-1] ............... ................ complex[k].real = data[k] complex[k].imag = data[N-k] ............... ................ complex[N/2].real = data[N/2] complex[N/2].real = 0 ............... ................ complex[k'].real = data[k] k' = N - k complex[k'].imag = -data[N-k] ............... ................ complex[N-1].real = data[1] complex[N-1].imag = -data[N-1] The radix-2 FFT functions for halfcomplex data are declared in the header file `gsl_fft_halfcomplex.h'. - Function: int gsl_fft_halfcomplex_radix2_inverse (double data[], const unsigned int n) - Function: int gsl_fft_halfcomplex_radix2_backwards (double data[], const unsigned int n) These functions compute the inverse or backwards in-place radix-2 FFT of the half-complex sequence DATA, a real of length N stored according the output scheme used by `gsl_fft_real_radix2'. The result is a real array stored in natural order. File: gsl-ref.info, Node: Mixed-radix FFT routines for real data, Next: Example of using mixed-radix FFT routines for real data, Prev: Radix-2 FFT routines for real data, Up: FFTs of real data Mixed-radix FFT routines for real data -------------------------------------- This section describes mixed-radix FFT algorithms for real data. The mixed-radix functions work for FFTs of any length. They are a reimplementation of the real-FFT routines in the Fortran FFTPACK library by Paul Swarztrauber. The theory behind the algorithm is explained in the article `Fast Mixed-Radix Real Fourier Transforms' by Clive Temperton. The routines here use the same indexing scheme and basic algorithms as FFTPACK. The functions use the FFTPACK storage convention for half-complex sequences. In this convention the half-complex transform of a real sequence is stored in with frequencies in increasing order, starting at zero, with the real and imaginary parts of each frequency in neighboring locations. When a value is known to be real the imaginary part is not stored. The imaginary part of the zero-frequency component is never stored. It is known to be zero (since the zero frequency component is simply the sum of the input data (all real)). For a sequence of even length the imaginary part of the frequency n/2 is not stored either, since the symmetry z_k = z_{N-k}^* implies that this is purely real too. The storage scheme is best shown by some examples. The table below shows the output for an odd-length sequence, n=5. The two columns give the correspondence between the 5 values in the half-complex sequence returned by `gsl_fft_real', HALFCOMPLEX[] and the values COMPLEX[] that would be returned if the same real input sequence were passed to `gsl_fft_complex_backward' as a complex sequence (with imaginary parts set to `0'), complex[0].real = halfcomplex[0] complex[0].imag = 0 complex[1].real = halfcomplex[1] complex[1].imag = halfcomplex[2] complex[2].real = halfcomplex[3] complex[2].imag = halfcomplex[4] complex[3].real = halfcomplex[3] complex[3].imag = -halfcomplex[4] complex[4].real = halfcomplex[1] complex[4].imag = -halfcomplex[2] The upper elements of the COMPLEX array, `complex[3]' and `complex[4]' are filled in using the symmetry condition. The imaginary part of the zero-frequency term `complex[0].imag' is known to be zero by the symmetry. The next table shows the output for an even-length sequence, n=5 In the even case both the there are two values which are purely real, complex[0].real = halfcomplex[0] complex[0].imag = 0 complex[1].real = halfcomplex[1] complex[1].imag = halfcomplex[2] complex[2].real = halfcomplex[3] complex[2].imag = halfcomplex[4] complex[3].real = halfcomplex[5] complex[3].imag = 0 complex[4].real = halfcomplex[3] complex[4].imag = -halfcomplex[4] complex[5].real = halfcomplex[1] complex[5].imag = -halfcomplex[2] The upper elements of the COMPLEX array, `complex[4]' and `complex[5]' are be filled in using the symmetry condition. Both `complex[0].imag' and `complex[3].imag' are known to be zero. All these functions are declared in the header files `gsl_fft_real.h' and `gsl_fft_halfcomplex.h'. - Function: int gsl_fft_real_wavetable_alloc (unsigned int N, gsl_fft_real_wavetable * WAVETABLE); - Function: int gsl_fft_halfcomplex_wavetable_alloc (unsigned int N, gsl_fft_halfcomplex_wavetable * WAVETABLE); These functions allocate space for a scratch area of size N real elements and a trigonometric lookup table, of size n/2 complex elements. The relevant components of WAVETABLE are modified to point to the newly allocated memory. These functions return a value of `0' if no errors were detected, and `-1' in the case of error. The following `gsl_errno' conditions are defined for these functions: `GSL_EDOM' The length of the data N is not a positive integer (i.e. N is zero). `GSL_ENOMEM' The requested memory could not be allocated (`malloc' returned a null pointer). - Function: int gsl_fft_real_init (const unsigned int N, gsl_fft_real_wavetable * WAVETABLE); - Function: int gsl_fft_halfcomplex_init (const unsigned int N, gsl_fft_halfcomplex_wavetable * WAVETABLE); These functions initialize WAVETABLE. They first select a factorization of the length N into the implemented subtransforms, storing the details of the factorization in WAVETABLE. Using this factorization they then prepare a trigonometric lookup table in the memory previously allocated by `gsl_fft_real_wavetable_alloc' or `gsl_fft_halfcomplex_wavetable_alloc'. The wavetable is computed using direct calls to `sin' and `cos', for accuracy. It could be computed faster using recursion relations if necessary. If an application performs many FFTs of the same length then computing the wavetable is a one-off overhead which does not affect the final throughput. The wavetable structure can be used repeatedly for any transform of the same length. The table is not modified by calls to any of the other FFT functions. The same wavetable can be used for both forward and backward (or inverse) transforms of a given length. The functions return a value of `0' if no errors were detected, and `-1' in the case of error. The following `gsl_errno' conditions are defined for these functions: `GSL_EDOM' The length of the data N is not a positive integer (i.e. N is zero). `GSL_EFACTOR' The length N could not be factorized (this shouldn't happen). `GSL_EFAILED' A failure was detected in the wavetable generation. This could be caused by an inconsistency in a user-supplied WAVETABLE structure. - Function: int gsl_fft_real_wavetable_free (gsl_fft_real_wavetable * WAVETABLE); - Function: int gsl_fft_halfcomplex_wavetable_free (gsl_fft_halfcomplex_wavetable * WAVETABLE); These functions free the blocks of memory previously allocated by `gsl_fft_real_wavetable_alloc' or `gsl_fft_halfcomplex_wavetable_alloc' and pointed to by the components of WAVETABLE. The wavetable should be freed if no further FFTs of the same length will be needed. - Function: int gsl_fft_real (double DATA[], const unsigned int N, const gsl_fft_real_wavetable * WAVETABLE) - Function: int gsl_fft_halfcomplex (double DATA[], const unsigned int N, const gsl_fft_halfcomplex_wavetable * WAVETABLE) These functions compute the FFT of DATA, a real or half-complex array of length N, using a mixed radix decimation-in-frequency algorithm. For `gsl_fft_real' DATA is an array of time-ordered real data. For `gsl_fft_halfcomplex' DATA contains fourier coefficients in the half-complex ordering described above. The caller must supply a WAVETABLE containing the chosen factorization, trigonometric lookup tables and scratch area. The wavetable can be easily prepared using the functions `gsl_fft_real_alloc' and `gsl_fft_real_init' or `gsl_fft_halfcomplex_alloc' and `gsl_fft_halfcomplex_init'. There is no restriction on the length N. Efficient modules are provided for subtransforms of length 2, 3, 4 and 5. Any remaining factors are computed with a slow, O(n^2), general-n module. The functions return a value of `0' if no errors were detected, and `-1' in the case of error. The following `gsl_errno' conditions are defined for these functions: `GSL_EDOM' The length of the data N is not a positive integer (i.e. N is zero). `GSL_EINVAL' The length of the data N and the length used to compute the given WAVETABLE do not match. - Function: int gsl_fft_real_unpack (const double real_coefficient[], complex complex_coefficient[], const unsigned int n) This function converts a single real array, REAL_COEFFICIENT into an equivalent complex array, COMPLEX_COEFFICIENT, (with imaginary part set to zero), suitable for `gsl_fft_complex' routines. The algorithm for the conversion is simply, for (i = 0; i < n; i++) { complex_coefficient[i].real = real_coefficient[i]; complex_coefficient[i].imag = 0.0; } The function returns a value of `0' if no errors were detected, and `-1' in the case of error. There is only one `gsl_errno' condition defined for this function: `GSL_EDOM' The length of the data N is not a positive integer (i.e. N is zero). - Function: int gsl_fft_halfcomplex_unpack (const double halfcomplex_coefficient[], complex complex_coefficient[], const unsigned int n) This function converts HALFCOMPLEX_COEFFICIENT, an array of half-complex coefficients as returned by `gsl_fft_real', into an ordinary complex array, COMPLEX_COEFFICIENT. It fills in the complex array using the symmetry z_k = z_{N-k}^* to reconstruct the redundant elements. The algorithm for the conversion is, complex_coefficient[0].real = halfcomplex_coefficient[0]; complex_coefficient[0].imag = 0.0; for (i = 1; i < n - i; i++) { const double hc_real = halfcomplex_coefficient[2 * i - 1]; const double hc_imag = halfcomplex_coefficient[2 * i]; complex_coefficient[i].real = hc_real; complex_coefficient[i].imag = hc_imag; complex_coefficient[n - i].real = hc_real; complex_coefficient[n - i].imag = -hc_imag; } if (i == n - i) { complex_coefficient[i].real = halfcomplex_coefficient[n - 1]; complex_coefficient[i].imag = 0.0; } The function returns a value of `0' if no errors were detected, and `-1' in the case of error. There is only one `gsl_errno' condition defined for this function: `GSL_EDOM' The length of the data N is not a positive integer (i.e. N is zero). File: gsl-ref.info, Node: Example of using mixed-radix FFT routines for real data, Prev: Mixed-radix FFT routines for real data, Up: FFTs of real data Example of using mixed-radix FFT routines for real data ------------------------------------------------------- Here is an example program using `gsl_fft_real' and `gsl_fft_halfcomplex'. It generates a real signal in the shape of a square pulse. The pulse is fourier transformed to frequency space, and all but the lowest ten frequency components are removed from the array of fourier coefficients returned by `gsl_fft_real'. The remaining fourier coefficients are transformed back to the time-domain, to give a filtered version of the square pulse. Since fourier coefficients are stored using the half-complex symmetry both positive and negative frequencies are removed and the final filtered signal is also real. #include #include #include #include #include int main () { int i, n = 100; double data[n]; gsl_fft_real_wavetable real_wavetable; gsl_fft_halfcomplex_wavetable halfcomplex_wavetable; for (i = 0; i < n; i++) { data[i] = 0.0; } for (i = n / 3; i < 2 * n / 3; i++) { data[i] = 1.0; } for (i = 0; i < n; i++) { printf ("%d: %e\n", i, data[i]); } printf ("\n"); gsl_fft_real_wavetable_alloc (n, &real_wavetable); gsl_fft_real_init (n, &real_wavetable); gsl_fft_real (data, n, &real_wavetable); gsl_fft_real_wavetable_free (&real_wavetable); for (i = 11; i < n; i++) { data[i] = 0; } gsl_fft_halfcomplex_wavetable_alloc (n, &halfcomplex_wavetable); gsl_fft_halfcomplex_init (n, &halfcomplex_wavetable); gsl_fft_halfcomplex_inverse (data, n, &halfcomplex_wavetable); gsl_fft_halfcomplex_wavetable_free (&halfcomplex_wavetable); for (i = 0; i < n; i++) { printf ("%d: %e\n", i, data[i]); } } File: gsl-ref.info, Node: Test Signals, Next: FFT References and Further Reading, Prev: FFTs of real data, Up: FFTs Test Signals ============ This section describes functions for generating simple test signals, such as unit pulses, constant signals, complex exponentials (sine and cosine waves) and pairs of complex exponentials with different frequencies. The Fourier transforms for these signals can be calculated analytically. Each function returns both the test signal and its Fourier transform. There are also some routines for generating white noise. In this case the transform cannot be calculated analytically, and is computed using a brute-force discrete Fourier transform. This may be slow but it is easy to check and independent of all the other routines. These functions are declared in the header file `gsl_fft_signals.h'. - Function: int gsl_fft_signal_complex_pulse (const unsigned int t, const unsigned int n, const double z_real, const double z_imag, complex data[], complex fft[]) This function prepares a "delta-function" pulse of complex amplitude (Z_REAL,Z_IMAG) at time T in the complex array DATA of length N. The pulse only occupies one array element. All other elements of DATA are set to zero. The function also returns the fourier transform of the pulse in the complex array FFT, also of length N. Writing z as (z_{real} + i z_{imag}), the mathematical definition of the pulse is, z_k = z \delta_{kt}, where \delta_{kt}=1 if k=t, and 0 otherwise. The fourier transform is, x_j = \sum_{k=0}^{n-1} z \delta_{kt} \exp(-2\pi i j k/n) = z \exp(-2\pi i j t/n) - Function: int gsl_fft_signal_complex_constant (const unsigned int n, const double z_real, const double z_imag, complex data[], complex fft[]) This function prepares a constant complex signal of amplitude (Z_REAL,Z_IMAG) in the complex array DATA of length N. The function also returns the fourier transform of the constant in the complex array FFT, also of length N. The mathematical definition of the constant signal is, z_k = z and the fourier transform is, x_j = \sum_{k=0}^{n-1} z \exp(-2\pi i j k/n) = z n \delta_{j0} - Function: int gsl_fft_signal_complex_exp (const int f, const unsigned int n, const double z_real, const double z_imag, complex data[], complex fft[]) This function prepares a complex exponential signal of amplitude (Z_REAL,Z_IMAG) and frequency F in the complex array DATA of length N. The function also returns the fourier transform of the signal in the complex array FFT, also of length N. The mathematical definition of the complex exponential signal is, z_k = z \exp(2\pi i k f/n) and the fourier transform is, x_j = \sum_{k=0}^{n-1} z \exp(2\pi i k f/n) \exp(-2\pi i j k/n) = z n \delta_{j,f} - Function: int gsl_fft_signal_complex_exppair (const int f1, const int f2, const unsigned int n, const double z1_real, const double z1_imag, const double z2_real, const double z2_imag, complex data[], complex fft[]) This function prepares a pair of exponential signals of amplitude (Z1_REAL,Z1_IMAG) and (Z2_REAL,Z2_IMAG) with frequencies F1 and F2 in the complex array DATA of length N. The function also returns the fourier transform of the signal in the complex array FFT, also of length N. The mathematical definition of the pair of exponential signals is, z_k = z_1 \exp(2\pi i k f_1/n) + z_2 \exp(2\pi i k f_2/n) and the fourier transform is, x_j = z_1 n \delta_{j,f_1} + z_2 n \delta_{j,f_2} - Function: int gsl_fft_signal_complex_noise (const unsigned int n, complex data[], complex fft[]) This function prepares a complex white-noise signal in the complex array DATA of length N. The real and imaginary components of each array element are independent uniform distributed random values in the range [0,1] (computed using the system `rand' function). The function also returns the fourier transform of the signal in the complex array FFT, also of length N. The fourier transform is computed using the brute-force discrete fourier transform `gsl_dft_complex_forward'. Note that the direct computation of the discrete fourier transform has worse numerical stability than an FFT, so the array FFT is only useful in detecting large errors. - Function: int gsl_fft_signal_real_noise (const unsigned int n, complex data[], complex fft[]) This function generates real white noise instead of complex white noise. It is the same as `gsl_fft_signal_complex_noise' but with the imaginary parts of DATA explicitly set to zero. All these functions return a value of `0' if no errors were detected, and `-1' in the case of error. There is only one `gsl_errno' condition defined for these functions: `GSL_EDOM' The length of the data N is not a positive integer (i.e. N is zero). File: gsl-ref.info, Node: FFT References and Further Reading, Prev: Test Signals, Up: FFTs FFT References and Further Reading ================================== A good starting point for learning more about the FFT is the review article `Fast Fourier Transforms: A Tutorial Review and A State of the Art' by Duhamel and Vetterli, P. Duhamel and M. Vetterli. Fast fourier transforms: A tutorial review and a state of the art. `Signal Processing', 19:259-299, 1990. To find out about the algorithms used in the GSL routines you may want to consult the latex document `GSL FFT Algorithms' (it is included in GSL, as `doc/fftalgorithms.tex'). This has general information on FFTs and explicit derivations of the implementation for each routine. There are also references to the relevant literature. For convenience some of the more important references are reproduced below. There are several introductory books on the FFT with example programs, such as `The Fast Fourier Transform' by Brigham and `DFT/FFT and Convolution Algorithms' by Burrus and Parks, E. Oran Brigham. `The Fast Fourier Transform'. Prentice Hall, 1974. C. S. Burrus and T. W. Parks. `DFT/FFT and Convolution Algorithms'. Wiley, 1984. Both these introductory books cover the radix-2 FFT in some detail. The mixed-radix algorithm at the heart of the FFTPACK routines is reviewed in Clive Temperton's paper, Clive Temperton. Self-sorting mixed-radix fast fourier transforms. `Journal of Computational Physics', 52(1):1-23, 1983. The derivation of FFTs for real-valued data is explained in the following two articles, Henrik V. Sorenson, Douglas L. Jones, Michael T. Heideman, and C. Sidney Burrus. Real-valued fast fourier transform algorithms. `IEEE Transactions on Acoustics, Speech, and Signal Processing', ASSP-35(6):849-863, 1987. Clive Temperton. Fast mixed-radix real fourier transforms. `Journal of Computational Physics', 52:340-350, 1983. In 1979 the IEEE published a compendium of carefully-reviewed Fortran FFT programs in `Programs for Digital Signal Processing'. It is a useful reference for implementations of many different FFT algorithms, Digital Signal Processing Committee and IEEE Acoustics, Speech, and Signal Processing Committee, editors. `Programs for Digital Signal Processing'. IEEE Press, 1979. There is also an annotated bibliography of papers on the FFT and related topics by Burrus, C. S. Burrus. Notes on the FFT. The notes are available from `http://www-dsp.rice.edu/res/fft/fftnote.asc'. File: gsl-ref.info, Node: Root finding, Next: Simulated Annealing, Prev: FFTs, Up: Top Root finding ************ This chapter describes functions for finding a root of an arbitrary function which you provide. It discusses proper use and possible pitfalls of each function and gives an overview of the algorithms involved. The *GNU Scientific Library* provides two types of root finding functions: high-level functions which take few arguments and make most decisions for you, and low-level functions which take many arguments and can be precisely controlled. For most tasks, the simpler high-level functions will be sufficient, but the greater degree of control provided by the low-level functions is sometimes useful. For examples of GSL root finding functions in action, look at the source code for the testing program which verifies that the they work correctly (it is run when you do a `make check'). It is in the file `roots/test.c' and implements one way of properly handling errors. The header file `gsl_roots.h' contains prototypes for GSL root finding functions and some other related declarations. *FIXME: add a note about performance.* * Menu: * Root Finding Overview:: * High Level vs. Low Level:: * High Level Root Finding Functions:: * Low Level Root Finding Functions:: * Root Finder Error Handling:: File: gsl-ref.info, Node: Root Finding Overview, Next: High Level vs. Low Level, Up: Root finding Root Finding Overview ===================== *FIXME: Insert discussion of numerical root finding here.* GSL root finding functions operate on continous functions only. While it is not absolutely required that f have a root within the search region, numerical root finding functions should not be used haphazardly to check for the *existence* of roots. There are better ways to do this! Because it is so easy to create situations where numerical root finders go awry, it is a bad idea to throw a root finder at a function you do not know much about. Note that GSL root finding functions can only search for one root at a time. While their behavior in cases of multiple roots (several roots in the search area, not one root of degree greater than 1) is deterministic, it depends entirely on the algorithm and the function under search and is rather complex; it is therefore difficult to predict. *In most cases, no error will be reported if you try to find a root in an area where there is more than one.* File: gsl-ref.info, Node: High Level vs. Low Level, Next: High Level Root Finding Functions, Prev: Root Finding Overview, Up: Root finding High Level vs. Low Level ======================== *FIXME: fill in this section.* In general, you should use high level functions unless there is a definite reason to use low level functions. File: gsl-ref.info, Node: High Level Root Finding Functions, Next: Low Level Root Finding Functions, Prev: High Level vs. Low Level, Up: Root finding High Level Root Finding Functions ================================= GSL provides three high level root finding functions. They are simple to use, robust, and adequate for most tasks, but certain situations require more control than the high level root finding functions provide. * Menu: * Root Finder Exit Values:: How to tell if an error occured and how the root location is returned * Providing the Function to Search:: How to provide a function for the root finder to operate on * Automatic Control Decisions:: How high level root finders set parameters such as tolerance * High Level Functions:: Names and arguments of the high level root finding functions File: gsl-ref.info, Node: Root Finder Exit Values, Next: Providing the Function to Search, Up: High Level Root Finding Functions Root Finder Exit Values ----------------------- Since the return value of GSL root finding functions is reserved for the error status, you must provide storage for the location of the found root. - Function Argument: double * root A pointer to a place for the GSL root finder to store the location of the found root. This must be a valid pointer; GSL root finders will not allocate any memory for you. If a GSL root finder succeeds, it will return `0' and store the location of the found root in `*root'. If a GSL root finder fails, it will return `-1' and set `gsl_errno' to a diagnostic value. *Note Root Finder Error Handling::, for a discussion of possible error codes. Nothing useful will be stored in `*root' if the function failed. File: gsl-ref.info, Node: Providing the Function to Search, Next: Automatic Control Decisions, Prev: Root Finder Exit Values, Up: High Level Root Finding Functions Providing the Function to Search -------------------------------- You must provide a continous function of one variable for GSL root finder(s) to operate on, and, sometimes, its first derivative. Recall that when passing pointers to functions, you give the name of the function you are passing. For example: foo = i_take_a_function_pointer(my_function); - Function Argument: double (* f)(double) A pointer to the function whose root you are searching for. It is called by the root finding function many times during its search. It must be continous within the region of interest. Here is an example function which you could pass to a GSL root finder: double my_f (double x) { return sin (2 * x) + 2 * cos (x); } - Function Argument: double (* df)(double) A pointer to the first derivative of the function whose root you are searching for. If we were looking for a root of the function in the previous example, this is what we would use for `df': double my_df (double x) { return 2 * cos (2 * x) - 2 * sin (x); } - Function Argument: void (* fdf)(double *, double *, double, int, int) A pointer to a function which calculates both the value of the function under search and the value of its first derivative. Because many terms of a function and its derivative are the same, it is often faster to use this method as opposed to providing f(x) and f'(x) separately. However, it is more complicated. It stores f(x) in its first argument and f'(x) in its second. Here's an example where f(x) = 2\sin(2x)\cos(x): void my_fdf (double * y, double * yprime, double x, int y_wanted, int yprime_wanted) { double sin2x, cosx; sin2x = sin (2 * x); cosx = cos (x); if (y_wanted) *y = 2 * sin2x * cos (x); if (yprime_wanted) *yprime = 2 * sin2x * -sin (x) + 2 * cos (2 * x) * cosx); } File: gsl-ref.info, Node: Automatic Control Decisions, Next: High Level Functions, Prev: Providing the Function to Search, Up: High Level Root Finding Functions Automatic Control Decisions --------------------------- *FIXME: Coming soon...* File: gsl-ref.info, Node: High Level Functions, Prev: Automatic Control Decisions, Up: High Level Root Finding Functions High Level Functions -------------------- *FIXME: Coming soon...* File: gsl-ref.info, Node: Low Level Root Finding Functions, Next: Root Finder Error Handling, Prev: High Level Root Finding Functions, Up: Root finding Low Level Root Finding Functions ================================ GSL provides several low level root finding functions which are more complicated to use than the high level function but provide more control. The following nodes explain how to use low level root finding functions. Low level functions return errors and roots and are provided functions to search in the same manner as the high level functions; see *Note Root Finder Exit Values::, and *Note Providing the Function to Search::, respectively. * Menu: * Search Bounds and Guesses:: How to tell GSL where to search * Search Stopping Parameters:: How to tell GSL when to stop searching * Bisection:: GSL's implementation of bisection * False Position:: GSL's implementation of false position * Secant Method:: GSL's implementation of secant method * Newtons Method:: GSL's implementation of Newton's Method File: gsl-ref.info, Node: Search Bounds and Guesses, Next: Search Stopping Parameters, Up: Low Level Root Finding Functions Search Bounds and Guesses ------------------------- When using low level functions, you can specify and monitor the region being searched more precisely than you can when using high level functions. You provide either search bounds or one or two guesses; this node explains how search bounds and guesses work and how function arguments control them. Search bounds are the endpoints of the search interval which is iterated smaller and smaller until the length of the interval is smaller than the requested precision or one of the endpoints converges; a guess is an x value which is iterated around until the it is within the desired precision of a root. Two guesses behave similarly to one; there are just two x values wandering about instead of one. In low level functions, these arguments are defined as pointers to `double' rather than simply `double's for two reasons. First, if the root finding function fails, it is very useful to have the final values of your iterated variables available to help diagnose why it failed. Second, it makes it possible to preserve the state of the root finder, enabling it to be restarted in the same place if needed. A situation where this could be useful is if the function under search is very costly to evaluate. Pretend, for a moment, that you want to find a root of a function which takes several minutes to evaluate. Your program could chug away for a few iterations, then examine your mailbox to see if the number of flames about bogging down the machine is above a certain threshold. If it isn't, chug away for a few more iterations and repeat the process; if it is, depending on your temperament, your program could renice itself to a lower priority and mail back apologies before continuing, or it could mailbomb everyone who complained, spawn copies of itself endlessly, and dump garbage to `/dev/audio'. (The latter approach is not recommended.) Note that these arguments must be valid pointers; GSL root finders will not allocate any memory for you. - Low Level Function Argument: double * lower_bound - Low Level Function Argument: double * upper_bound The initial upper and lower bounds of the interval in which to search for a root. `lower_bound' must be less than `upper_bound'. These arguments are modified during execution of the root finding function; if you need to preserve their initial values, you must make copies of them. See the third paragraph of this node for the reasoning behind this behavior. - Low Level Function Argument: double * guess - Low Level Function Argument: double * guess2 One or two initial values for the guess(es) iterated by the root finding function. These arguments are modified during execution of the root finding function; if you need to preserve their initial values, you must make copies of them. See the third paragraph of this node for the reasoning behind this behavior. File: gsl-ref.info, Node: Search Stopping Parameters, Next: Bisection, Prev: Search Bounds and Guesses, Up: Low Level Root Finding Functions Search Stopping Parameters -------------------------- *FIXME: clean up this section and add documentation for max_deltay.* GSL root finding functions (and numerical root finding functions in general) stop when one of the following conditions is true: * A root has been found to within the user-specified precision. * A user-specified maximum number of iterations has executed. * An error has occured. Whenever you call a low level GSL root finding function, you must specify precisely absolute and/or relative tolerances and the maximum number of iterations. GSL decides that two values a and b. with relative tolerance R and absolute tolerance A, are close enough if the following relation is true: |a - b| \leq R \min(|a|,|b|) + A You can set either R or A to zero, but be aware that GSL will signal an error if the search moves into an area where both R and A are meaningless; assuming a and b are the endpoints of the region of interest, the following must be true or GSL will signal an error: R \min(|a|,|b|) + A \geq 10 \max(|a|,|b|) \times {\tt DBL\_EPSILON} (We introduce a buffer of 10 to protect against roundoff error.) For the sake of efficient resource use, do not ask for more precision than you need, especially if your function is costly to evaluate. - Low Level Function Argument: double abs_epsilon The maximum permissible absolute error in GSL root finder answers. The only static limit on `abs_epsilon' is that it must be positive; see above for other restrictions, however. - Low Level Function Argument: double rel_epsilon The maximum permissible relative error in GSL root finder answers. `rel_epsilon' must be greater than or equal to {\tt DBL\_EPSILON} \times 10 (note the buffer factor to protect against roundoff error). See above for additional non-static restrictions. - Function Argument: unsigned int max_iterations The maximum number of iterations a root finder is allowed to perform. This must be greater than or equal to 1, as performing a negative number of iterations is extremely difficult and not doing any iterations is rather useless. Do not set `max_iterations' too large. If there is a problem, you want to know about it as soon as possible; you don't want your program chugging away for many cycles in error. In addition, GSL root finding functions which extrapolate (Newton's Method, (*Note Newtons Method::, and Secant Method, *Note Secant Method::) accept an additional argument: - Function Argument: double max_step_size The maximum step size an extrapolating algorithm is allowed to take. This is to prevents the algorithm from landing on a place where the test function's derivative is very small and zooming off to infinity or into a different solution basin. *FIXME: talk about minimum value for max_step_size.* For example, if while solving \sin(x) = 0, x_n of Newton's Method (*note Newtons Method::.) landed on 1.570700000 (\pi/2 \approx 1.570796327), then x_{n+1} would be approximately -10000, which is definitely not what we wanted! We want the root finder to recognize this step as "too big" and flag an error. The alarm bell will ring if the following relation in true: |{{d} \over {dx}} f(x)| < |{{f(x)} \over {\tt max\_step\_size}}| Note that while Secant Method (*note Secant Method::.) does not deal with derivatives directly, when extrapolating it approximates them numerically. Do not set `max_step_size' too large; that will defeat its purpose. In the \sin(x) = 0 example, \pi would be a good value for `max_step_size'; any step larger than that would certainly be headed astray. A good understanding of the problem is especially important for `max_step_size'. File: gsl-ref.info, Node: Bisection, Next: False Position, Prev: Search Stopping Parameters, Up: Low Level Root Finding Functions Bisection --------- Bisection is a simple and robust method of finding roots of a function f; when its arguments are valid, it cannot fail. However, it is the slowest algorithm implemented by GSL, and it cannot find roots of even degree. Its convergence is linear. One begins the algorithm with an interval which is guaranteed by the Intermediate Value Theorem to contain a root: where a and b are the endpoints of the interval, f(a) must differ in sign from f(b). (If you're a bit fuzzy on the Intermediate Value Theorem, consult any elementary calculus textbook.) Each iteration, bisection chops its interval in half and discards the interval which does not contain a root. Once the interval is smaller than the requested epsilon, iteration stops and the root location is returned. (Often, root finding algorithms are better explained geometrically. The TeX version of this documentation contains an image illustrating bisection.) - Function: int gsl_root_bisection (double * ROOT, double (* F)(double), double * LOWER_BOUND, double * UPPER_BOUND, double REL_EPSILON, double ABS_EPSILON, unsigned int MAX_ITERATIONS, double MAX_DELTAY) Search for a zero of `f' using bisection. Returns `0' if successful, `-1' on error (*note Root Finder Error Handling::.). `f(*lower_bound)' and `f(*upper_bound)' must differ in sign. Arguments: `root' A place to store the root location once it is found. *Note Root Finder Exit Values::. `f' A user defined function to search for a root. *Note Providing the Function to Search::. `lower_bound, upper_bound' Lower and upper bounds of the interval to search. *Note Search Bounds and Guesses::. `rel_epsilon, abs_epsilon' Maximum permitted relative and absolute error. *Note Search Stopping Parameters::. `max_iterations' The maximum allowed number of iterations. *Note Search Stopping Parameters::. `max_deltay' The maximum allowed difference between `f(*lower_bound)' and `f(*upper_bound)'. *Note Search Stopping Parameters::. File: gsl-ref.info, Node: False Position, Next: Secant Method, Prev: Bisection, Up: Low Level Root Finding Functions False Position -------------- False position is a robust method of finding roots of a function f; if its arguments are valid, it cannot fail. However, it cannot find roots of even degree. Its convergence is linear, but it is usually faster than bisection. One begins the algorithm with an interval which is guaranteed by the Intermediate Value Theorem to contain a root: where a and b are the endpoints of the interval, f(a) must differ in sign from f(b). (If you're a bit fuzzy on the Intermediate Value Theorem, consult any elementary calculus textbook.) Each iteration, false position draws a line between f(a) and f(b); the x position where this line crosses the x axis is where the interval is split. The part of the interval which contains the root is taken to be the new interval, and the process is repeated until one of the following is true: |a - b| \leq \varepsilon a_n - a_{n-1} = 0 \;\;\;{\rm and}\;\;\; |b_n - b_{n-1}| \leq \varepsilon b_n - b_{n-1} = 0 \;\;\;{\rm and}\;\;\;|a_n - a_{n-1}| \leq \varepsilon (Often, root finding algorithms are better explained geometrically. The TeX version of this documentation contains an image illustrating false position.) - Function: int gsl_root_falsepos (double * ROOT, double (* F)(double), double * LOWER_BOUND, double * UPPER_BOUND, double REL_EPSILON, double ABS_EPSILON, unsigned int MAX_ITERATIONS, double MAX_DELTAY) Search for a zero of `f' using false position. Return `0' if successful, `-1' on error (*note Root Finder Error Handling::.. `f(*lower_bound)' and `f(*upper_bound)' must differ in sign. Arguments: `root' A place to store the root location once it is found. *Note Root Finder Exit Values::. `f' A user defined function to search for a root. *Note Providing the Function to Search::. `lower_bound, upper_bound' Lower and upper bounds of the interval to search. *Note Search Bounds and Guesses::. `rel_epsilon, abs_epsilon' Maximum permitted relative and absolute error. *Note Search Stopping Parameters::. `max_iterations' The maximum allowed number of iterations. *Note Search Stopping Parameters::. `max_deltay' The maximum allowed difference between `f(*lower_bound)' and `f(*upper_bound)'. *Note Search Stopping Parameters::. File: gsl-ref.info, Node: Secant Method, Next: Newtons Method, Prev: False Position, Up: Low Level Root Finding Functions Secant Method ------------- Secant Method is a somewhat fragile method of finding roots. On single roots, its convergence is of order (1 + \sqrt 5)/2 (approximately 1.62). On multiple roots, converges linearly. One begins the algorithm with two guesses for the value of the root, g_0 and g_1. The root may be either inside or outside the interval defined by g_0 and g_1. Each iteration, Secant Method draws a line through f(g_{n - 1}) and f(g_n). The x position where this line crosses the x axis becomes g_{n + 1}. g_{n - 1} is discarded, n is incremented, and the process is repeated until: |g_n - g_{n-1}| \leq \varepsilon Note that g_{n + 1} may be obtained by either interpolation or extrapolation and that Secant Method cannot fail during interpolation. Secant Method breaks in the same situations that Newton's Method does, though it is somewhat less sensitive. (*Note Newtons Method::.) (Often, root finding algorithms are better explained geometrically. The TeX version of this documentation contains an image illustrating Secant Method.) - Function: int gsl_root_secant (double * ROOT, double (* F)(double), double * GUESS, double * GUESS2, double EPSILON, unsigned int MAX_ITERATIONS, double MAX_STEP_SIZE) Search for a zero of `f' using Secant Method, with `*guess' and `*guess2' being the guesses. arguments. *Note Root Finder Error Handling::, for a discussion of possible error codes.