1 This is gsl-ref.info, produced by makeinfo version 4.8 from
4 INFO-DIR-SECTION Scientific software
6 * gsl-ref: (gsl-ref). GNU Scientific Library - Reference
9 Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2004,
10 2005, 2006, 2007 The GSL Team.
12 Permission is granted to copy, distribute and/or modify this document
13 under the terms of the GNU Free Documentation License, Version 1.2 or
14 any later version published by the Free Software Foundation; with the
15 Invariant Sections being "GNU General Public License" and "Free Software
16 Needs Free Documentation", the Front-Cover text being "A GNU Manual",
17 and with the Back-Cover Text being (a) (see below). A copy of the
18 license is included in the section entitled "GNU Free Documentation
21 (a) The Back-Cover Text is: "You have the freedom to copy and modify
25 File: gsl-ref.info, Node: Random Number Distribution Examples, Next: Random Number Distribution References and Further Reading, Prev: Shuffling and Sampling, Up: Random Number Distributions
30 The following program demonstrates the use of a random number generator
31 to produce variates from a distribution. It prints 10 samples from the
32 Poisson distribution with a mean of 3.
35 #include <gsl/gsl_rng.h>
36 #include <gsl/gsl_randist.h>
41 const gsl_rng_type * T;
47 /* create a generator chosen by the
48 environment variable GSL_RNG_TYPE */
53 r = gsl_rng_alloc (T);
55 /* print n random variates chosen from
56 the poisson distribution with mean
59 for (i = 0; i < n; i++)
61 unsigned int k = gsl_ran_poisson (r, mu);
70 If the library and header files are installed under `/usr/local' (the
71 default location) then the program can be compiled with these options,
73 $ gcc -Wall demo.c -lgsl -lgslcblas -lm
75 Here is the output of the program,
80 The variates depend on the seed used by the generator. The seed for the
81 default generator type `gsl_rng_default' can be changed with the
82 `GSL_RNG_SEED' environment variable to produce a different stream of
85 $ GSL_RNG_SEED=123 ./a.out
89 The following program generates a random walk in two dimensions.
92 #include <gsl/gsl_rng.h>
93 #include <gsl/gsl_randist.h>
99 double x = 0, y = 0, dx, dy;
101 const gsl_rng_type * T;
106 r = gsl_rng_alloc (T);
108 printf ("%g %g\n", x, y);
110 for (i = 0; i < 10; i++)
112 gsl_ran_dir_2d (r, &dx, &dy);
114 printf ("%g %g\n", x, y);
121 Here is the output from the program, three 10-step random walks from
124 The following program computes the upper and lower cumulative
125 distribution functions for the standard normal distribution at x=2.
128 #include <gsl/gsl_cdf.h>
136 P = gsl_cdf_ugaussian_P (x);
137 printf ("prob(x < %f) = %f\n", x, P);
139 Q = gsl_cdf_ugaussian_Q (x);
140 printf ("prob(x > %f) = %f\n", x, Q);
142 x = gsl_cdf_ugaussian_Pinv (P);
143 printf ("Pinv(%f) = %f\n", P, x);
145 x = gsl_cdf_ugaussian_Qinv (Q);
146 printf ("Qinv(%f) = %f\n", Q, x);
151 Here is the output of the program,
153 prob(x < 2.000000) = 0.977250
154 prob(x > 2.000000) = 0.022750
155 Pinv(0.977250) = 2.000000
156 Qinv(0.022750) = 2.000000
159 File: gsl-ref.info, Node: Random Number Distribution References and Further Reading, Prev: Random Number Distribution Examples, Up: Random Number Distributions
161 19.40 References and Further Reading
162 ====================================
164 For an encyclopaedic coverage of the subject readers are advised to
165 consult the book `Non-Uniform Random Variate Generation' by Luc
166 Devroye. It covers every imaginable distribution and provides hundreds
169 Luc Devroye, `Non-Uniform Random Variate Generation',
170 Springer-Verlag, ISBN 0-387-96305-7.
172 The subject of random variate generation is also reviewed by Knuth, who
173 describes algorithms for all the major distributions.
175 Donald E. Knuth, `The Art of Computer Programming: Seminumerical
176 Algorithms' (Vol 2, 3rd Ed, 1997), Addison-Wesley, ISBN 0201896842.
178 The Particle Data Group provides a short review of techniques for
179 generating distributions of random numbers in the "Monte Carlo" section
180 of its Annual Review of Particle Physics.
182 `Review of Particle Properties' R.M. Barnett et al., Physical
183 Review D54, 1 (1996) `http://pdg.lbl.gov/'.
185 The Review of Particle Physics is available online in postscript and pdf
188 An overview of methods used to compute cumulative distribution functions
189 can be found in `Statistical Computing' by W.J. Kennedy and J.E.
190 Gentle. Another general reference is `Elements of Statistical
191 Computing' by R.A. Thisted.
193 William E. Kennedy and James E. Gentle, `Statistical Computing'
194 (1980), Marcel Dekker, ISBN 0-8247-6898-1.
196 Ronald A. Thisted, `Elements of Statistical Computing' (1988),
197 Chapman & Hall, ISBN 0-412-01371-1.
199 The cumulative distribution functions for the Gaussian distribution are
200 based on the following papers,
202 `Rational Chebyshev Approximations Using Linear Equations', W.J.
203 Cody, W. Fraser, J.F. Hart. Numerische Mathematik 12, 242-251
206 `Rational Chebyshev Approximations for the Error Function', W.J.
207 Cody. Mathematics of Computation 23, n107, 631-637 (July 1969).
210 File: gsl-ref.info, Node: Statistics, Next: Histograms, Prev: Random Number Distributions, Up: Top
215 This chapter describes the statistical functions in the library. The
216 basic statistical functions include routines to compute the mean,
217 variance and standard deviation. More advanced functions allow you to
218 calculate absolute deviations, skewness, and kurtosis as well as the
219 median and arbitrary percentiles. The algorithms use recurrence
220 relations to compute average quantities in a stable way, without large
221 intermediate values that might overflow.
223 The functions are available in versions for datasets in the standard
224 floating-point and integer types. The versions for double precision
225 floating-point data have the prefix `gsl_stats' and are declared in the
226 header file `gsl_statistics_double.h'. The versions for integer data
227 have the prefix `gsl_stats_int' and are declared in the header file
228 `gsl_statistics_int.h'.
232 * Mean and standard deviation and variance::
233 * Absolute deviation::
234 * Higher moments (skewness and kurtosis)::
239 * Maximum and Minimum values::
240 * Median and Percentiles::
241 * Example statistical programs::
242 * Statistics References and Further Reading::
245 File: gsl-ref.info, Node: Mean and standard deviation and variance, Next: Absolute deviation, Up: Statistics
247 20.1 Mean, Standard Deviation and Variance
248 ==========================================
250 -- Function: double gsl_stats_mean (const double DATA[], size_t
252 This function returns the arithmetic mean of DATA, a dataset of
253 length N with stride STRIDE. The arithmetic mean, or "sample
254 mean", is denoted by \Hat\mu and defined as,
256 \Hat\mu = (1/N) \sum x_i
258 where x_i are the elements of the dataset DATA. For samples drawn
259 from a gaussian distribution the variance of \Hat\mu is \sigma^2 /
262 -- Function: double gsl_stats_variance (const double DATA[], size_t
264 This function returns the estimated, or "sample", variance of
265 DATA, a dataset of length N with stride STRIDE. The estimated
266 variance is denoted by \Hat\sigma^2 and is defined by,
268 \Hat\sigma^2 = (1/(N-1)) \sum (x_i - \Hat\mu)^2
270 where x_i are the elements of the dataset DATA. Note that the
271 normalization factor of 1/(N-1) results from the derivation of
272 \Hat\sigma^2 as an unbiased estimator of the population variance
273 \sigma^2. For samples drawn from a gaussian distribution the
274 variance of \Hat\sigma^2 itself is 2 \sigma^4 / N.
276 This function computes the mean via a call to `gsl_stats_mean'. If
277 you have already computed the mean then you can pass it directly to
278 `gsl_stats_variance_m'.
280 -- Function: double gsl_stats_variance_m (const double DATA[], size_t
281 STRIDE, size_t N, double MEAN)
282 This function returns the sample variance of DATA relative to the
283 given value of MEAN. The function is computed with \Hat\mu
284 replaced by the value of MEAN that you supply,
286 \Hat\sigma^2 = (1/(N-1)) \sum (x_i - mean)^2
288 -- Function: double gsl_stats_sd (const double DATA[], size_t STRIDE,
290 -- Function: double gsl_stats_sd_m (const double DATA[], size_t
291 STRIDE, size_t N, double MEAN)
292 The standard deviation is defined as the square root of the
293 variance. These functions return the square root of the
294 corresponding variance functions above.
296 -- Function: double gsl_stats_tss (const double DATA[], size_t STRIDE,
298 -- Function: double gsl_stats_tss_m (const double DATA[], size_t
299 STRIDE, size_t N, double MEAN)
300 These functions return the total sum of squares (TSS) of DATA about
301 the mean. For `gsl_stats_tss_m' the user-supplied value of MEAN
302 is used, and for `gsl_stats_tss' it is computed using
305 TSS = \sum (x_i - mean)^2
307 -- Function: double gsl_stats_variance_with_fixed_mean (const double
308 DATA[], size_t STRIDE, size_t N, double MEAN)
309 This function computes an unbiased estimate of the variance of
310 DATA when the population mean MEAN of the underlying distribution
311 is known _a priori_. In this case the estimator for the variance
312 uses the factor 1/N and the sample mean \Hat\mu is replaced by the
313 known population mean \mu,
315 \Hat\sigma^2 = (1/N) \sum (x_i - \mu)^2
317 -- Function: double gsl_stats_sd_with_fixed_mean (const double DATA[],
318 size_t STRIDE, size_t N, double MEAN)
319 This function calculates the standard deviation of DATA for a
320 fixed population mean MEAN. The result is the square root of the
321 corresponding variance function.
324 File: gsl-ref.info, Node: Absolute deviation, Next: Higher moments (skewness and kurtosis), Prev: Mean and standard deviation and variance, Up: Statistics
326 20.2 Absolute deviation
327 =======================
329 -- Function: double gsl_stats_absdev (const double DATA[], size_t
331 This function computes the absolute deviation from the mean of
332 DATA, a dataset of length N with stride STRIDE. The absolute
333 deviation from the mean is defined as,
335 absdev = (1/N) \sum |x_i - \Hat\mu|
337 where x_i are the elements of the dataset DATA. The absolute
338 deviation from the mean provides a more robust measure of the
339 width of a distribution than the variance. This function computes
340 the mean of DATA via a call to `gsl_stats_mean'.
342 -- Function: double gsl_stats_absdev_m (const double DATA[], size_t
343 STRIDE, size_t N, double MEAN)
344 This function computes the absolute deviation of the dataset DATA
345 relative to the given value of MEAN,
347 absdev = (1/N) \sum |x_i - mean|
349 This function is useful if you have already computed the mean of
350 DATA (and want to avoid recomputing it), or wish to calculate the
351 absolute deviation relative to another value (such as zero, or the
355 File: gsl-ref.info, Node: Higher moments (skewness and kurtosis), Next: Autocorrelation, Prev: Absolute deviation, Up: Statistics
357 20.3 Higher moments (skewness and kurtosis)
358 ===========================================
360 -- Function: double gsl_stats_skew (const double DATA[], size_t
362 This function computes the skewness of DATA, a dataset of length N
363 with stride STRIDE. The skewness is defined as,
365 skew = (1/N) \sum ((x_i - \Hat\mu)/\Hat\sigma)^3
367 where x_i are the elements of the dataset DATA. The skewness
368 measures the asymmetry of the tails of a distribution.
370 The function computes the mean and estimated standard deviation of
371 DATA via calls to `gsl_stats_mean' and `gsl_stats_sd'.
373 -- Function: double gsl_stats_skew_m_sd (const double DATA[], size_t
374 STRIDE, size_t N, double MEAN, double SD)
375 This function computes the skewness of the dataset DATA using the
376 given values of the mean MEAN and standard deviation SD,
378 skew = (1/N) \sum ((x_i - mean)/sd)^3
380 These functions are useful if you have already computed the mean
381 and standard deviation of DATA and want to avoid recomputing them.
383 -- Function: double gsl_stats_kurtosis (const double DATA[], size_t
385 This function computes the kurtosis of DATA, a dataset of length N
386 with stride STRIDE. The kurtosis is defined as,
388 kurtosis = ((1/N) \sum ((x_i - \Hat\mu)/\Hat\sigma)^4) - 3
390 The kurtosis measures how sharply peaked a distribution is,
391 relative to its width. The kurtosis is normalized to zero for a
392 gaussian distribution.
394 -- Function: double gsl_stats_kurtosis_m_sd (const double DATA[],
395 size_t STRIDE, size_t N, double MEAN, double SD)
396 This function computes the kurtosis of the dataset DATA using the
397 given values of the mean MEAN and standard deviation SD,
399 kurtosis = ((1/N) \sum ((x_i - mean)/sd)^4) - 3
401 This function is useful if you have already computed the mean and
402 standard deviation of DATA and want to avoid recomputing them.
405 File: gsl-ref.info, Node: Autocorrelation, Next: Covariance, Prev: Higher moments (skewness and kurtosis), Up: Statistics
410 -- Function: double gsl_stats_lag1_autocorrelation (const double
411 DATA[], const size_t STRIDE, const size_t N)
412 This function computes the lag-1 autocorrelation of the dataset
415 a_1 = {\sum_{i = 1}^{n} (x_{i} - \Hat\mu) (x_{i-1} - \Hat\mu)
417 \sum_{i = 1}^{n} (x_{i} - \Hat\mu) (x_{i} - \Hat\mu)}
419 -- Function: double gsl_stats_lag1_autocorrelation_m (const double
420 DATA[], const size_t STRIDE, const size_t N, const double
422 This function computes the lag-1 autocorrelation of the dataset
423 DATA using the given value of the mean MEAN.
427 File: gsl-ref.info, Node: Covariance, Next: Correlation, Prev: Autocorrelation, Up: Statistics
432 -- Function: double gsl_stats_covariance (const double DATA1[], const
433 size_t STRIDE1, const double DATA2[], const size_t STRIDE2,
435 This function computes the covariance of the datasets DATA1 and
436 DATA2 which must both be of the same length N.
438 covar = (1/(n - 1)) \sum_{i = 1}^{n} (x_i - \Hat x) (y_i - \Hat y)
440 -- Function: double gsl_stats_covariance_m (const double DATA1[],
441 const size_t STRIDE1, const double DATA2[], const size_t
442 STRIDE2, const size_t N, const double MEAN1, const double
444 This function computes the covariance of the datasets DATA1 and
445 DATA2 using the given values of the means, MEAN1 and MEAN2. This
446 is useful if you have already computed the means of DATA1 and
447 DATA2 and want to avoid recomputing them.
450 File: gsl-ref.info, Node: Correlation, Next: Weighted Samples, Prev: Covariance, Up: Statistics
455 -- Function: double gsl_stats_correlation (const double DATA1[], const
456 size_t STRIDE1, const double DATA2[], const size_t STRIDE2,
458 This function efficiently computes the Pearson correlation
459 coefficient between the datasets DATA1 and DATA2 which must both
460 be of the same length N.
461 r = cov(x, y) / (\Hat\sigma_x \Hat\sigma_y)
462 = {1/(n-1) \sum (x_i - \Hat x) (y_i - \Hat y)
464 \sqrt{1/(n-1) \sum (x_i - \Hat x)^2} \sqrt{1/(n-1) \sum (y_i - \Hat y)^2}
468 File: gsl-ref.info, Node: Weighted Samples, Next: Maximum and Minimum values, Prev: Correlation, Up: Statistics
470 20.7 Weighted Samples
471 =====================
473 The functions described in this section allow the computation of
474 statistics for weighted samples. The functions accept an array of
475 samples, x_i, with associated weights, w_i. Each sample x_i is
476 considered as having been drawn from a Gaussian distribution with
477 variance \sigma_i^2. The sample weight w_i is defined as the
478 reciprocal of this variance, w_i = 1/\sigma_i^2. Setting a weight to
479 zero corresponds to removing a sample from a dataset.
481 -- Function: double gsl_stats_wmean (const double W[], size_t WSTRIDE,
482 const double DATA[], size_t STRIDE, size_t N)
483 This function returns the weighted mean of the dataset DATA with
484 stride STRIDE and length N, using the set of weights W with stride
485 WSTRIDE and length N. The weighted mean is defined as,
487 \Hat\mu = (\sum w_i x_i) / (\sum w_i)
489 -- Function: double gsl_stats_wvariance (const double W[], size_t
490 WSTRIDE, const double DATA[], size_t STRIDE, size_t N)
491 This function returns the estimated variance of the dataset DATA
492 with stride STRIDE and length N, using the set of weights W with
493 stride WSTRIDE and length N. The estimated variance of a weighted
494 dataset is defined as,
496 \Hat\sigma^2 = ((\sum w_i)/((\sum w_i)^2 - \sum (w_i^2)))
497 \sum w_i (x_i - \Hat\mu)^2
499 Note that this expression reduces to an unweighted variance with
500 the familiar 1/(N-1) factor when there are N equal non-zero
503 -- Function: double gsl_stats_wvariance_m (const double W[], size_t
504 WSTRIDE, const double DATA[], size_t STRIDE, size_t N, double
506 This function returns the estimated variance of the weighted
507 dataset DATA using the given weighted mean WMEAN.
509 -- Function: double gsl_stats_wsd (const double W[], size_t WSTRIDE,
510 const double DATA[], size_t STRIDE, size_t N)
511 The standard deviation is defined as the square root of the
512 variance. This function returns the square root of the
513 corresponding variance function `gsl_stats_wvariance' above.
515 -- Function: double gsl_stats_wsd_m (const double W[], size_t WSTRIDE,
516 const double DATA[], size_t STRIDE, size_t N, double WMEAN)
517 This function returns the square root of the corresponding variance
518 function `gsl_stats_wvariance_m' above.
520 -- Function: double gsl_stats_wvariance_with_fixed_mean (const double
521 W[], size_t WSTRIDE, const double DATA[], size_t STRIDE,
522 size_t N, const double MEAN)
523 This function computes an unbiased estimate of the variance of
524 weighted dataset DATA when the population mean MEAN of the
525 underlying distribution is known _a priori_. In this case the
526 estimator for the variance replaces the sample mean \Hat\mu by the
527 known population mean \mu,
529 \Hat\sigma^2 = (\sum w_i (x_i - \mu)^2) / (\sum w_i)
531 -- Function: double gsl_stats_wsd_with_fixed_mean (const double W[],
532 size_t WSTRIDE, const double DATA[], size_t STRIDE, size_t N,
534 The standard deviation is defined as the square root of the
535 variance. This function returns the square root of the
536 corresponding variance function above.
538 -- Function: double gsl_stats_wtss (const double W[], const size_t
539 WSTRIDE, const double DATA[], size_t STRIDE, size_t N)
540 -- Function: double gsl_stats_wtss_m (const double W[], const size_t
541 WSTRIDE, const double DATA[], size_t STRIDE, size_t N, double
543 These functions return the weighted total sum of squares (TSS) of
544 DATA about the weighted mean. For `gsl_stats_wtss_m' the
545 user-supplied value of WMEAN is used, and for `gsl_stats_wtss' it
546 is computed using `gsl_stats_wmean'.
548 TSS = \sum w_i (x_i - wmean)^2
550 -- Function: double gsl_stats_wabsdev (const double W[], size_t
551 WSTRIDE, const double DATA[], size_t STRIDE, size_t N)
552 This function computes the weighted absolute deviation from the
553 weighted mean of DATA. The absolute deviation from the mean is
556 absdev = (\sum w_i |x_i - \Hat\mu|) / (\sum w_i)
558 -- Function: double gsl_stats_wabsdev_m (const double W[], size_t
559 WSTRIDE, const double DATA[], size_t STRIDE, size_t N, double
561 This function computes the absolute deviation of the weighted
562 dataset DATA about the given weighted mean WMEAN.
564 -- Function: double gsl_stats_wskew (const double W[], size_t WSTRIDE,
565 const double DATA[], size_t STRIDE, size_t N)
566 This function computes the weighted skewness of the dataset DATA.
568 skew = (\sum w_i ((x_i - xbar)/\sigma)^3) / (\sum w_i)
570 -- Function: double gsl_stats_wskew_m_sd (const double W[], size_t
571 WSTRIDE, const double DATA[], size_t STRIDE, size_t N, double
573 This function computes the weighted skewness of the dataset DATA
574 using the given values of the weighted mean and weighted standard
575 deviation, WMEAN and WSD.
577 -- Function: double gsl_stats_wkurtosis (const double W[], size_t
578 WSTRIDE, const double DATA[], size_t STRIDE, size_t N)
579 This function computes the weighted kurtosis of the dataset DATA.
581 kurtosis = ((\sum w_i ((x_i - xbar)/sigma)^4) / (\sum w_i)) - 3
583 -- Function: double gsl_stats_wkurtosis_m_sd (const double W[], size_t
584 WSTRIDE, const double DATA[], size_t STRIDE, size_t N, double
586 This function computes the weighted kurtosis of the dataset DATA
587 using the given values of the weighted mean and weighted standard
588 deviation, WMEAN and WSD.
591 File: gsl-ref.info, Node: Maximum and Minimum values, Next: Median and Percentiles, Prev: Weighted Samples, Up: Statistics
593 20.8 Maximum and Minimum values
594 ===============================
596 The following functions find the maximum and minimum values of a
597 dataset (or their indices). If the data contains `NaN's then a `NaN'
598 will be returned, since the maximum or minimum value is undefined. For
599 functions which return an index, the location of the first `NaN' in the
602 -- Function: double gsl_stats_max (const double DATA[], size_t STRIDE,
604 This function returns the maximum value in DATA, a dataset of
605 length N with stride STRIDE. The maximum value is defined as the
606 value of the element x_i which satisfies x_i >= x_j for all j.
608 If you want instead to find the element with the largest absolute
609 magnitude you will need to apply `fabs' or `abs' to your data
610 before calling this function.
612 -- Function: double gsl_stats_min (const double DATA[], size_t STRIDE,
614 This function returns the minimum value in DATA, a dataset of
615 length N with stride STRIDE. The minimum value is defined as the
616 value of the element x_i which satisfies x_i <= x_j for all j.
618 If you want instead to find the element with the smallest absolute
619 magnitude you will need to apply `fabs' or `abs' to your data
620 before calling this function.
622 -- Function: void gsl_stats_minmax (double * MIN, double * MAX, const
623 double DATA[], size_t STRIDE, size_t N)
624 This function finds both the minimum and maximum values MIN, MAX
625 in DATA in a single pass.
627 -- Function: size_t gsl_stats_max_index (const double DATA[], size_t
629 This function returns the index of the maximum value in DATA, a
630 dataset of length N with stride STRIDE. The maximum value is
631 defined as the value of the element x_i which satisfies x_i >= x_j
632 for all j. When there are several equal maximum elements then the
635 -- Function: size_t gsl_stats_min_index (const double DATA[], size_t
637 This function returns the index of the minimum value in DATA, a
638 dataset of length N with stride STRIDE. The minimum value is
639 defined as the value of the element x_i which satisfies x_i >= x_j
640 for all j. When there are several equal minimum elements then the
643 -- Function: void gsl_stats_minmax_index (size_t * MIN_INDEX, size_t *
644 MAX_INDEX, const double DATA[], size_t STRIDE, size_t N)
645 This function returns the indexes MIN_INDEX, MAX_INDEX of the
646 minimum and maximum values in DATA in a single pass.
649 File: gsl-ref.info, Node: Median and Percentiles, Next: Example statistical programs, Prev: Maximum and Minimum values, Up: Statistics
651 20.9 Median and Percentiles
652 ===========================
654 The median and percentile functions described in this section operate on
655 sorted data. For convenience we use "quantiles", measured on a scale
656 of 0 to 1, instead of percentiles (which use a scale of 0 to 100).
658 -- Function: double gsl_stats_median_from_sorted_data (const double
659 SORTED_DATA[], size_t STRIDE, size_t N)
660 This function returns the median value of SORTED_DATA, a dataset
661 of length N with stride STRIDE. The elements of the array must be
662 in ascending numerical order. There are no checks to see whether
663 the data are sorted, so the function `gsl_sort' should always be
666 When the dataset has an odd number of elements the median is the
667 value of element (n-1)/2. When the dataset has an even number of
668 elements the median is the mean of the two nearest middle values,
669 elements (n-1)/2 and n/2. Since the algorithm for computing the
670 median involves interpolation this function always returns a
671 floating-point number, even for integer data types.
673 -- Function: double gsl_stats_quantile_from_sorted_data (const double
674 SORTED_DATA[], size_t STRIDE, size_t N, double F)
675 This function returns a quantile value of SORTED_DATA, a
676 double-precision array of length N with stride STRIDE. The
677 elements of the array must be in ascending numerical order. The
678 quantile is determined by the F, a fraction between 0 and 1. For
679 example, to compute the value of the 75th percentile F should have
682 There are no checks to see whether the data are sorted, so the
683 function `gsl_sort' should always be used first.
685 The quantile is found by interpolation, using the formula
687 quantile = (1 - \delta) x_i + \delta x_{i+1}
689 where i is `floor'((n - 1)f) and \delta is (n-1)f - i.
691 Thus the minimum value of the array (`data[0*stride]') is given by
692 F equal to zero, the maximum value (`data[(n-1)*stride]') is given
693 by F equal to one and the median value is given by F equal to 0.5.
694 Since the algorithm for computing quantiles involves
695 interpolation this function always returns a floating-point
696 number, even for integer data types.
699 File: gsl-ref.info, Node: Example statistical programs, Next: Statistics References and Further Reading, Prev: Median and Percentiles, Up: Statistics
704 Here is a basic example of how to use the statistical functions:
707 #include <gsl/gsl_statistics.h>
712 double data[5] = {17.2, 18.1, 16.5, 18.3, 12.6};
713 double mean, variance, largest, smallest;
715 mean = gsl_stats_mean(data, 1, 5);
716 variance = gsl_stats_variance(data, 1, 5);
717 largest = gsl_stats_max(data, 1, 5);
718 smallest = gsl_stats_min(data, 1, 5);
720 printf ("The dataset is %g, %g, %g, %g, %g\n",
721 data[0], data[1], data[2], data[3], data[4]);
723 printf ("The sample mean is %g\n", mean);
724 printf ("The estimated variance is %g\n", variance);
725 printf ("The largest value is %g\n", largest);
726 printf ("The smallest value is %g\n", smallest);
730 The program should produce the following output,
732 The dataset is 17.2, 18.1, 16.5, 18.3, 12.6
733 The sample mean is 16.54
734 The estimated variance is 4.2984
735 The largest value is 18.3
736 The smallest value is 12.6
738 Here is an example using sorted data,
741 #include <gsl/gsl_sort.h>
742 #include <gsl/gsl_statistics.h>
747 double data[5] = {17.2, 18.1, 16.5, 18.3, 12.6};
748 double median, upperq, lowerq;
750 printf ("Original dataset: %g, %g, %g, %g, %g\n",
751 data[0], data[1], data[2], data[3], data[4]);
753 gsl_sort (data, 1, 5);
755 printf ("Sorted dataset: %g, %g, %g, %g, %g\n",
756 data[0], data[1], data[2], data[3], data[4]);
759 = gsl_stats_median_from_sorted_data (data,
763 = gsl_stats_quantile_from_sorted_data (data,
767 = gsl_stats_quantile_from_sorted_data (data,
771 printf ("The median is %g\n", median);
772 printf ("The upper quartile is %g\n", upperq);
773 printf ("The lower quartile is %g\n", lowerq);
777 This program should produce the following output,
779 Original dataset: 17.2, 18.1, 16.5, 18.3, 12.6
780 Sorted dataset: 12.6, 16.5, 17.2, 18.1, 18.3
782 The upper quartile is 18.1
783 The lower quartile is 16.5
786 File: gsl-ref.info, Node: Statistics References and Further Reading, Prev: Example statistical programs, Up: Statistics
788 20.11 References and Further Reading
789 ====================================
791 The standard reference for almost any topic in statistics is the
792 multi-volume `Advanced Theory of Statistics' by Kendall and Stuart.
794 Maurice Kendall, Alan Stuart, and J. Keith Ord. `The Advanced
795 Theory of Statistics' (multiple volumes) reprinted as `Kendall's
796 Advanced Theory of Statistics'. Wiley, ISBN 047023380X.
798 Many statistical concepts can be more easily understood by a Bayesian
799 approach. The following book by Gelman, Carlin, Stern and Rubin gives a
800 comprehensive coverage of the subject.
802 Andrew Gelman, John B. Carlin, Hal S. Stern, Donald B. Rubin.
803 `Bayesian Data Analysis'. Chapman & Hall, ISBN 0412039915.
805 For physicists the Particle Data Group provides useful reviews of
806 Probability and Statistics in the "Mathematical Tools" section of its
807 Annual Review of Particle Physics.
809 `Review of Particle Properties' R.M. Barnett et al., Physical
812 The Review of Particle Physics is available online at the website
813 `http://pdg.lbl.gov/'.
816 File: gsl-ref.info, Node: Histograms, Next: N-tuples, Prev: Statistics, Up: Top
821 This chapter describes functions for creating histograms. Histograms
822 provide a convenient way of summarizing the distribution of a set of
823 data. A histogram consists of a set of "bins" which count the number of
824 events falling into a given range of a continuous variable x. In GSL
825 the bins of a histogram contain floating-point numbers, so they can be
826 used to record both integer and non-integer distributions. The bins
827 can use arbitrary sets of ranges (uniformly spaced bins are the
828 default). Both one and two-dimensional histograms are supported.
830 Once a histogram has been created it can also be converted into a
831 probability distribution function. The library provides efficient
832 routines for selecting random samples from probability distributions.
833 This can be useful for generating simulations based on real data.
835 The functions are declared in the header files `gsl_histogram.h' and
840 * The histogram struct::
841 * Histogram allocation::
842 * Copying Histograms::
843 * Updating and accessing histogram elements::
844 * Searching histogram ranges::
845 * Histogram Statistics::
846 * Histogram Operations::
847 * Reading and writing histograms::
848 * Resampling from histograms::
849 * The histogram probability distribution struct::
850 * Example programs for histograms::
851 * Two dimensional histograms::
852 * The 2D histogram struct::
853 * 2D Histogram allocation::
854 * Copying 2D Histograms::
855 * Updating and accessing 2D histogram elements::
856 * Searching 2D histogram ranges::
857 * 2D Histogram Statistics::
858 * 2D Histogram Operations::
859 * Reading and writing 2D histograms::
860 * Resampling from 2D histograms::
861 * Example programs for 2D histograms::
864 File: gsl-ref.info, Node: The histogram struct, Next: Histogram allocation, Up: Histograms
866 21.1 The histogram struct
867 =========================
869 A histogram is defined by the following struct,
871 -- Data Type: gsl_histogram
873 This is the number of histogram bins
876 The ranges of the bins are stored in an array of N+1 elements
880 The counts for each bin are stored in an array of N elements
881 pointed to by BIN. The bins are floating-point numbers, so
882 you can increment them by non-integer values if necessary.
884 The range for BIN[i] is given by RANGE[i] to RANGE[i+1]. For n bins
885 there are n+1 entries in the array RANGE. Each bin is inclusive at the
886 lower end and exclusive at the upper end. Mathematically this means
887 that the bins are defined by the following inequality,
888 bin[i] corresponds to range[i] <= x < range[i+1]
890 Here is a diagram of the correspondence between ranges and bins on the
894 [ bin[0] )[ bin[1] )[ bin[2] )[ bin[3] )[ bin[4] )
895 ---|---------|---------|---------|---------|---------|--- x
896 r[0] r[1] r[2] r[3] r[4] r[5]
898 In this picture the values of the RANGE array are denoted by r. On the
899 left-hand side of each bin the square bracket `[' denotes an inclusive
900 lower bound (r <= x), and the round parentheses `)' on the right-hand
901 side denote an exclusive upper bound (x < r). Thus any samples which
902 fall on the upper end of the histogram are excluded. If you want to
903 include this value for the last bin you will need to add an extra bin
906 The `gsl_histogram' struct and its associated functions are defined
907 in the header file `gsl_histogram.h'.
910 File: gsl-ref.info, Node: Histogram allocation, Next: Copying Histograms, Prev: The histogram struct, Up: Histograms
912 21.2 Histogram allocation
913 =========================
915 The functions for allocating memory to a histogram follow the style of
916 `malloc' and `free'. In addition they also perform their own error
917 checking. If there is insufficient memory available to allocate a
918 histogram then the functions call the error handler (with an error
919 number of `GSL_ENOMEM') in addition to returning a null pointer. Thus
920 if you use the library error handler to abort your program then it
921 isn't necessary to check every histogram `alloc'.
923 -- Function: gsl_histogram * gsl_histogram_alloc (size_t N)
924 This function allocates memory for a histogram with N bins, and
925 returns a pointer to a newly created `gsl_histogram' struct. If
926 insufficient memory is available a null pointer is returned and the
927 error handler is invoked with an error code of `GSL_ENOMEM'. The
928 bins and ranges are not initialized, and should be prepared using
929 one of the range-setting functions below in order to make the
930 histogram ready for use.
932 -- Function: int gsl_histogram_set_ranges (gsl_histogram * H, const
933 double RANGE[], size_t SIZE)
934 This function sets the ranges of the existing histogram H using
935 the array RANGE of size SIZE. The values of the histogram bins
936 are reset to zero. The `range' array should contain the desired
937 bin limits. The ranges can be arbitrary, subject to the
938 restriction that they are monotonically increasing.
940 The following example shows how to create a histogram with
941 logarithmic bins with ranges [1,10), [10,100) and [100,1000).
943 gsl_histogram * h = gsl_histogram_alloc (3);
945 /* bin[0] covers the range 1 <= x < 10 */
946 /* bin[1] covers the range 10 <= x < 100 */
947 /* bin[2] covers the range 100 <= x < 1000 */
949 double range[4] = { 1.0, 10.0, 100.0, 1000.0 };
951 gsl_histogram_set_ranges (h, range, 4);
953 Note that the size of the RANGE array should be defined to be one
954 element bigger than the number of bins. The additional element is
955 required for the upper value of the final bin.
957 -- Function: int gsl_histogram_set_ranges_uniform (gsl_histogram * H,
958 double XMIN, double XMAX)
959 This function sets the ranges of the existing histogram H to cover
960 the range XMIN to XMAX uniformly. The values of the histogram
961 bins are reset to zero. The bin ranges are shown in the table
963 bin[0] corresponds to xmin <= x < xmin + d
964 bin[1] corresponds to xmin + d <= x < xmin + 2 d
966 bin[n-1] corresponds to xmin + (n-1)d <= x < xmax
968 where d is the bin spacing, d = (xmax-xmin)/n.
970 -- Function: void gsl_histogram_free (gsl_histogram * H)
971 This function frees the histogram H and all of the memory
975 File: gsl-ref.info, Node: Copying Histograms, Next: Updating and accessing histogram elements, Prev: Histogram allocation, Up: Histograms
977 21.3 Copying Histograms
978 =======================
980 -- Function: int gsl_histogram_memcpy (gsl_histogram * DEST, const
982 This function copies the histogram SRC into the pre-existing
983 histogram DEST, making DEST into an exact copy of SRC. The two
984 histograms must be of the same size.
986 -- Function: gsl_histogram * gsl_histogram_clone (const gsl_histogram
988 This function returns a pointer to a newly created histogram which
989 is an exact copy of the histogram SRC.
992 File: gsl-ref.info, Node: Updating and accessing histogram elements, Next: Searching histogram ranges, Prev: Copying Histograms, Up: Histograms
994 21.4 Updating and accessing histogram elements
995 ==============================================
997 There are two ways to access histogram bins, either by specifying an x
998 coordinate or by using the bin-index directly. The functions for
999 accessing the histogram through x coordinates use a binary search to
1000 identify the bin which covers the appropriate range.
1002 -- Function: int gsl_histogram_increment (gsl_histogram * H, double X)
1003 This function updates the histogram H by adding one (1.0) to the
1004 bin whose range contains the coordinate X.
1006 If X lies in the valid range of the histogram then the function
1007 returns zero to indicate success. If X is less than the lower
1008 limit of the histogram then the function returns `GSL_EDOM', and
1009 none of bins are modified. Similarly, if the value of X is greater
1010 than or equal to the upper limit of the histogram then the function
1011 returns `GSL_EDOM', and none of the bins are modified. The error
1012 handler is not called, however, since it is often necessary to
1013 compute histograms for a small range of a larger dataset, ignoring
1014 the values outside the range of interest.
1016 -- Function: int gsl_histogram_accumulate (gsl_histogram * H, double
1018 This function is similar to `gsl_histogram_increment' but increases
1019 the value of the appropriate bin in the histogram H by the
1020 floating-point number WEIGHT.
1022 -- Function: double gsl_histogram_get (const gsl_histogram * H, size_t
1024 This function returns the contents of the I-th bin of the histogram
1025 H. If I lies outside the valid range of indices for the histogram
1026 then the error handler is called with an error code of `GSL_EDOM'
1027 and the function returns 0.
1029 -- Function: int gsl_histogram_get_range (const gsl_histogram * H,
1030 size_t I, double * LOWER, double * UPPER)
1031 This function finds the upper and lower range limits of the I-th
1032 bin of the histogram H. If the index I is valid then the
1033 corresponding range limits are stored in LOWER and UPPER. The
1034 lower limit is inclusive (i.e. events with this coordinate are
1035 included in the bin) and the upper limit is exclusive (i.e. events
1036 with the coordinate of the upper limit are excluded and fall in the
1037 neighboring higher bin, if it exists). The function returns 0 to
1038 indicate success. If I lies outside the valid range of indices for
1039 the histogram then the error handler is called and the function
1040 returns an error code of `GSL_EDOM'.
1042 -- Function: double gsl_histogram_max (const gsl_histogram * H)
1043 -- Function: double gsl_histogram_min (const gsl_histogram * H)
1044 -- Function: size_t gsl_histogram_bins (const gsl_histogram * H)
1045 These functions return the maximum upper and minimum lower range
1046 limits and the number of bins of the histogram H. They provide a
1047 way of determining these values without accessing the
1048 `gsl_histogram' struct directly.
1050 -- Function: void gsl_histogram_reset (gsl_histogram * H)
1051 This function resets all the bins in the histogram H to zero.
1054 File: gsl-ref.info, Node: Searching histogram ranges, Next: Histogram Statistics, Prev: Updating and accessing histogram elements, Up: Histograms
1056 21.5 Searching histogram ranges
1057 ===============================
1059 The following functions are used by the access and update routines to
1060 locate the bin which corresponds to a given x coordinate.
1062 -- Function: int gsl_histogram_find (const gsl_histogram * H, double
1064 This function finds and sets the index I to the bin number which
1065 covers the coordinate X in the histogram H. The bin is located
1066 using a binary search. The search includes an optimization for
1067 histograms with uniform range, and will return the correct bin
1068 immediately in this case. If X is found in the range of the
1069 histogram then the function sets the index I and returns
1070 `GSL_SUCCESS'. If X lies outside the valid range of the histogram
1071 then the function returns `GSL_EDOM' and the error handler is
1075 File: gsl-ref.info, Node: Histogram Statistics, Next: Histogram Operations, Prev: Searching histogram ranges, Up: Histograms
1077 21.6 Histogram Statistics
1078 =========================
1080 -- Function: double gsl_histogram_max_val (const gsl_histogram * H)
1081 This function returns the maximum value contained in the histogram
1084 -- Function: size_t gsl_histogram_max_bin (const gsl_histogram * H)
1085 This function returns the index of the bin containing the maximum
1086 value. In the case where several bins contain the same maximum
1087 value the smallest index is returned.
1089 -- Function: double gsl_histogram_min_val (const gsl_histogram * H)
1090 This function returns the minimum value contained in the histogram
1093 -- Function: size_t gsl_histogram_min_bin (const gsl_histogram * H)
1094 This function returns the index of the bin containing the minimum
1095 value. In the case where several bins contain the same maximum
1096 value the smallest index is returned.
1098 -- Function: double gsl_histogram_mean (const gsl_histogram * H)
1099 This function returns the mean of the histogrammed variable, where
1100 the histogram is regarded as a probability distribution. Negative
1101 bin values are ignored for the purposes of this calculation. The
1102 accuracy of the result is limited by the bin width.
1104 -- Function: double gsl_histogram_sigma (const gsl_histogram * H)
1105 This function returns the standard deviation of the histogrammed
1106 variable, where the histogram is regarded as a probability
1107 distribution. Negative bin values are ignored for the purposes of
1108 this calculation. The accuracy of the result is limited by the bin
1111 -- Function: double gsl_histogram_sum (const gsl_histogram * H)
1112 This function returns the sum of all bin values. Negative bin
1113 values are included in the sum.
1116 File: gsl-ref.info, Node: Histogram Operations, Next: Reading and writing histograms, Prev: Histogram Statistics, Up: Histograms
1118 21.7 Histogram Operations
1119 =========================
1121 -- Function: int gsl_histogram_equal_bins_p (const gsl_histogram * H1,
1122 const gsl_histogram * H2)
1123 This function returns 1 if the all of the individual bin ranges of
1124 the two histograms are identical, and 0 otherwise.
1126 -- Function: int gsl_histogram_add (gsl_histogram * H1, const
1128 This function adds the contents of the bins in histogram H2 to the
1129 corresponding bins of histogram H1, i.e. h'_1(i) = h_1(i) +
1130 h_2(i). The two histograms must have identical bin ranges.
1132 -- Function: int gsl_histogram_sub (gsl_histogram * H1, const
1134 This function subtracts the contents of the bins in histogram H2
1135 from the corresponding bins of histogram H1, i.e. h'_1(i) = h_1(i)
1136 - h_2(i). The two histograms must have identical bin ranges.
1138 -- Function: int gsl_histogram_mul (gsl_histogram * H1, const
1140 This function multiplies the contents of the bins of histogram H1
1141 by the contents of the corresponding bins in histogram H2, i.e.
1142 h'_1(i) = h_1(i) * h_2(i). The two histograms must have identical
1145 -- Function: int gsl_histogram_div (gsl_histogram * H1, const
1147 This function divides the contents of the bins of histogram H1 by
1148 the contents of the corresponding bins in histogram H2, i.e.
1149 h'_1(i) = h_1(i) / h_2(i). The two histograms must have identical
1152 -- Function: int gsl_histogram_scale (gsl_histogram * H, double SCALE)
1153 This function multiplies the contents of the bins of histogram H
1154 by the constant SCALE, i.e. h'_1(i) = h_1(i) * scale.
1156 -- Function: int gsl_histogram_shift (gsl_histogram * H, double OFFSET)
1157 This function shifts the contents of the bins of histogram H by
1158 the constant OFFSET, i.e. h'_1(i) = h_1(i) + offset.
1161 File: gsl-ref.info, Node: Reading and writing histograms, Next: Resampling from histograms, Prev: Histogram Operations, Up: Histograms
1163 21.8 Reading and writing histograms
1164 ===================================
1166 The library provides functions for reading and writing histograms to a
1167 file as binary data or formatted text.
1169 -- Function: int gsl_histogram_fwrite (FILE * STREAM, const
1171 This function writes the ranges and bins of the histogram H to the
1172 stream STREAM in binary format. The return value is 0 for success
1173 and `GSL_EFAILED' if there was a problem writing to the file.
1174 Since the data is written in the native binary format it may not
1175 be portable between different architectures.
1177 -- Function: int gsl_histogram_fread (FILE * STREAM, gsl_histogram * H)
1178 This function reads into the histogram H from the open stream
1179 STREAM in binary format. The histogram H must be preallocated
1180 with the correct size since the function uses the number of bins
1181 in H to determine how many bytes to read. The return value is 0
1182 for success and `GSL_EFAILED' if there was a problem reading from
1183 the file. The data is assumed to have been written in the native
1184 binary format on the same architecture.
1186 -- Function: int gsl_histogram_fprintf (FILE * STREAM, const
1187 gsl_histogram * H, const char * RANGE_FORMAT, const char *
1189 This function writes the ranges and bins of the histogram H
1190 line-by-line to the stream STREAM using the format specifiers
1191 RANGE_FORMAT and BIN_FORMAT. These should be one of the `%g',
1192 `%e' or `%f' formats for floating point numbers. The function
1193 returns 0 for success and `GSL_EFAILED' if there was a problem
1194 writing to the file. The histogram output is formatted in three
1195 columns, and the columns are separated by spaces, like this,
1197 range[0] range[1] bin[0]
1198 range[1] range[2] bin[1]
1199 range[2] range[3] bin[2]
1201 range[n-1] range[n] bin[n-1]
1203 The values of the ranges are formatted using RANGE_FORMAT and the
1204 value of the bins are formatted using BIN_FORMAT. Each line
1205 contains the lower and upper limit of the range of the bins and the
1206 value of the bin itself. Since the upper limit of one bin is the
1207 lower limit of the next there is duplication of these values
1208 between lines but this allows the histogram to be manipulated with
1209 line-oriented tools.
1211 -- Function: int gsl_histogram_fscanf (FILE * STREAM, gsl_histogram *
1213 This function reads formatted data from the stream STREAM into the
1214 histogram H. The data is assumed to be in the three-column format
1215 used by `gsl_histogram_fprintf'. The histogram H must be
1216 preallocated with the correct length since the function uses the
1217 size of H to determine how many numbers to read. The function
1218 returns 0 for success and `GSL_EFAILED' if there was a problem
1219 reading from the file.
1222 File: gsl-ref.info, Node: Resampling from histograms, Next: The histogram probability distribution struct, Prev: Reading and writing histograms, Up: Histograms
1224 21.9 Resampling from histograms
1225 ===============================
1227 A histogram made by counting events can be regarded as a measurement of
1228 a probability distribution. Allowing for statistical error, the height
1229 of each bin represents the probability of an event where the value of x
1230 falls in the range of that bin. The probability distribution function
1231 has the one-dimensional form p(x)dx where,
1235 In this equation n_i is the number of events in the bin which contains
1236 x, w_i is the width of the bin and N is the total number of events.
1237 The distribution of events within each bin is assumed to be uniform.
1240 File: gsl-ref.info, Node: The histogram probability distribution struct, Next: Example programs for histograms, Prev: Resampling from histograms, Up: Histograms
1242 21.10 The histogram probability distribution struct
1243 ===================================================
1245 The probability distribution function for a histogram consists of a set
1246 of "bins" which measure the probability of an event falling into a
1247 given range of a continuous variable x. A probability distribution
1248 function is defined by the following struct, which actually stores the
1249 cumulative probability distribution function. This is the natural
1250 quantity for generating samples via the inverse transform method,
1251 because there is a one-to-one mapping between the cumulative
1252 probability distribution and the range [0,1]. It can be shown that by
1253 taking a uniform random number in this range and finding its
1254 corresponding coordinate in the cumulative probability distribution we
1255 obtain samples with the desired probability distribution.
1257 -- Data Type: gsl_histogram_pdf
1259 This is the number of bins used to approximate the probability
1260 distribution function.
1263 The ranges of the bins are stored in an array of N+1 elements
1264 pointed to by RANGE.
1267 The cumulative probability for the bins is stored in an array
1268 of N elements pointed to by SUM.
1270 The following functions allow you to create a `gsl_histogram_pdf'
1271 struct which represents this probability distribution and generate
1272 random samples from it.
1274 -- Function: gsl_histogram_pdf * gsl_histogram_pdf_alloc (size_t N)
1275 This function allocates memory for a probability distribution with
1276 N bins and returns a pointer to a newly initialized
1277 `gsl_histogram_pdf' struct. If insufficient memory is available a
1278 null pointer is returned and the error handler is invoked with an
1279 error code of `GSL_ENOMEM'.
1281 -- Function: int gsl_histogram_pdf_init (gsl_histogram_pdf * P, const
1283 This function initializes the probability distribution P with the
1284 contents of the histogram H. If any of the bins of H are negative
1285 then the error handler is invoked with an error code of `GSL_EDOM'
1286 because a probability distribution cannot contain negative values.
1288 -- Function: void gsl_histogram_pdf_free (gsl_histogram_pdf * P)
1289 This function frees the probability distribution function P and
1290 all of the memory associated with it.
1292 -- Function: double gsl_histogram_pdf_sample (const gsl_histogram_pdf
1294 This function uses R, a uniform random number between zero and
1295 one, to compute a single random sample from the probability
1296 distribution P. The algorithm used to compute the sample s is
1297 given by the following formula,
1299 s = range[i] + delta * (range[i+1] - range[i])
1301 where i is the index which satisfies sum[i] <= r < sum[i+1] and
1302 delta is (r - sum[i])/(sum[i+1] - sum[i]).
1305 File: gsl-ref.info, Node: Example programs for histograms, Next: Two dimensional histograms, Prev: The histogram probability distribution struct, Up: Histograms
1307 21.11 Example programs for histograms
1308 =====================================
1310 The following program shows how to make a simple histogram of a column
1311 of numerical data supplied on `stdin'. The program takes three
1312 arguments, specifying the upper and lower bounds of the histogram and
1313 the number of bins. It then reads numbers from `stdin', one line at a
1314 time, and adds them to the histogram. When there is no more data to
1315 read it prints out the accumulated histogram using
1316 `gsl_histogram_fprintf'.
1320 #include <gsl/gsl_histogram.h>
1323 main (int argc, char **argv)
1330 printf ("Usage: gsl-histogram xmin xmax n\n"
1331 "Computes a histogram of the data "
1332 "on stdin using n bins from xmin "
1343 gsl_histogram * h = gsl_histogram_alloc (n);
1344 gsl_histogram_set_ranges_uniform (h, a, b);
1346 while (fscanf (stdin, "%lg", &x) == 1)
1348 gsl_histogram_increment (h, x);
1350 gsl_histogram_fprintf (stdout, h, "%g", "%g");
1351 gsl_histogram_free (h);
1356 Here is an example of the program in use. We generate 10000 random
1357 samples from a Cauchy distribution with a width of 30 and histogram
1358 them over the range -100 to 100, using 200 bins.
1360 $ gsl-randist 0 10000 cauchy 30
1361 | gsl-histogram -100 100 200 > histogram.dat
1363 A plot of the resulting histogram shows the familiar shape of the
1364 Cauchy distribution and the fluctuations caused by the finite sample
1367 $ awk '{print $1, $3 ; print $2, $3}' histogram.dat
1371 File: gsl-ref.info, Node: Two dimensional histograms, Next: The 2D histogram struct, Prev: Example programs for histograms, Up: Histograms
1373 21.12 Two dimensional histograms
1374 ================================
1376 A two dimensional histogram consists of a set of "bins" which count the
1377 number of events falling in a given area of the (x,y) plane. The
1378 simplest way to use a two dimensional histogram is to record
1379 two-dimensional position information, n(x,y). Another possibility is
1380 to form a "joint distribution" by recording related variables. For
1381 example a detector might record both the position of an event (x) and
1382 the amount of energy it deposited E. These could be histogrammed as
1383 the joint distribution n(x,E).
1386 File: gsl-ref.info, Node: The 2D histogram struct, Next: 2D Histogram allocation, Prev: Two dimensional histograms, Up: Histograms
1388 21.13 The 2D histogram struct
1389 =============================
1391 Two dimensional histograms are defined by the following struct,
1393 -- Data Type: gsl_histogram2d
1395 This is the number of histogram bins in the x and y
1399 The ranges of the bins in the x-direction are stored in an
1400 array of NX + 1 elements pointed to by XRANGE.
1403 The ranges of the bins in the y-direction are stored in an
1404 array of NY + 1 elements pointed to by YRANGE.
1407 The counts for each bin are stored in an array pointed to by
1408 BIN. The bins are floating-point numbers, so you can
1409 increment them by non-integer values if necessary. The array
1410 BIN stores the two dimensional array of bins in a single
1411 block of memory according to the mapping `bin(i,j)' = `bin[i
1414 The range for `bin(i,j)' is given by `xrange[i]' to `xrange[i+1]' in
1415 the x-direction and `yrange[j]' to `yrange[j+1]' in the y-direction.
1416 Each bin is inclusive at the lower end and exclusive at the upper end.
1417 Mathematically this means that the bins are defined by the following
1419 bin(i,j) corresponds to xrange[i] <= x < xrange[i+1]
1420 and yrange[j] <= y < yrange[j+1]
1422 Note that any samples which fall on the upper sides of the histogram are
1423 excluded. If you want to include these values for the side bins you
1424 will need to add an extra row or column to your histogram.
1426 The `gsl_histogram2d' struct and its associated functions are
1427 defined in the header file `gsl_histogram2d.h'.
1430 File: gsl-ref.info, Node: 2D Histogram allocation, Next: Copying 2D Histograms, Prev: The 2D histogram struct, Up: Histograms
1432 21.14 2D Histogram allocation
1433 =============================
1435 The functions for allocating memory to a 2D histogram follow the style
1436 of `malloc' and `free'. In addition they also perform their own error
1437 checking. If there is insufficient memory available to allocate a
1438 histogram then the functions call the error handler (with an error
1439 number of `GSL_ENOMEM') in addition to returning a null pointer. Thus
1440 if you use the library error handler to abort your program then it
1441 isn't necessary to check every 2D histogram `alloc'.
1443 -- Function: gsl_histogram2d * gsl_histogram2d_alloc (size_t NX,
1445 This function allocates memory for a two-dimensional histogram with
1446 NX bins in the x direction and NY bins in the y direction. The
1447 function returns a pointer to a newly created `gsl_histogram2d'
1448 struct. If insufficient memory is available a null pointer is
1449 returned and the error handler is invoked with an error code of
1450 `GSL_ENOMEM'. The bins and ranges must be initialized with one of
1451 the functions below before the histogram is ready for use.
1453 -- Function: int gsl_histogram2d_set_ranges (gsl_histogram2d * H,
1454 const double XRANGE[], size_t XSIZE, const double YRANGE[],
1456 This function sets the ranges of the existing histogram H using
1457 the arrays XRANGE and YRANGE of size XSIZE and YSIZE respectively.
1458 The values of the histogram bins are reset to zero.
1460 -- Function: int gsl_histogram2d_set_ranges_uniform (gsl_histogram2d *
1461 H, double XMIN, double XMAX, double YMIN, double YMAX)
1462 This function sets the ranges of the existing histogram H to cover
1463 the ranges XMIN to XMAX and YMIN to YMAX uniformly. The values of
1464 the histogram bins are reset to zero.
1466 -- Function: void gsl_histogram2d_free (gsl_histogram2d * H)
1467 This function frees the 2D histogram H and all of the memory
1471 File: gsl-ref.info, Node: Copying 2D Histograms, Next: Updating and accessing 2D histogram elements, Prev: 2D Histogram allocation, Up: Histograms
1473 21.15 Copying 2D Histograms
1474 ===========================
1476 -- Function: int gsl_histogram2d_memcpy (gsl_histogram2d * DEST, const
1477 gsl_histogram2d * SRC)
1478 This function copies the histogram SRC into the pre-existing
1479 histogram DEST, making DEST into an exact copy of SRC. The two
1480 histograms must be of the same size.
1482 -- Function: gsl_histogram2d * gsl_histogram2d_clone (const
1483 gsl_histogram2d * SRC)
1484 This function returns a pointer to a newly created histogram which
1485 is an exact copy of the histogram SRC.
1488 File: gsl-ref.info, Node: Updating and accessing 2D histogram elements, Next: Searching 2D histogram ranges, Prev: Copying 2D Histograms, Up: Histograms
1490 21.16 Updating and accessing 2D histogram elements
1491 ==================================================
1493 You can access the bins of a two-dimensional histogram either by
1494 specifying a pair of (x,y) coordinates or by using the bin indices
1495 (i,j) directly. The functions for accessing the histogram through
1496 (x,y) coordinates use binary searches in the x and y directions to
1497 identify the bin which covers the appropriate range.
1499 -- Function: int gsl_histogram2d_increment (gsl_histogram2d * H,
1501 This function updates the histogram H by adding one (1.0) to the
1502 bin whose x and y ranges contain the coordinates (X,Y).
1504 If the point (x,y) lies inside the valid ranges of the histogram
1505 then the function returns zero to indicate success. If (x,y) lies
1506 outside the limits of the histogram then the function returns
1507 `GSL_EDOM', and none of the bins are modified. The error handler
1508 is not called, since it is often necessary to compute histograms
1509 for a small range of a larger dataset, ignoring any coordinates
1510 outside the range of interest.
1512 -- Function: int gsl_histogram2d_accumulate (gsl_histogram2d * H,
1513 double X, double Y, double WEIGHT)
1514 This function is similar to `gsl_histogram2d_increment' but
1515 increases the value of the appropriate bin in the histogram H by
1516 the floating-point number WEIGHT.
1518 -- Function: double gsl_histogram2d_get (const gsl_histogram2d * H,
1520 This function returns the contents of the (I,J)-th bin of the
1521 histogram H. If (I,J) lies outside the valid range of indices for
1522 the histogram then the error handler is called with an error code
1523 of `GSL_EDOM' and the function returns 0.
1525 -- Function: int gsl_histogram2d_get_xrange (const gsl_histogram2d *
1526 H, size_t I, double * XLOWER, double * XUPPER)
1527 -- Function: int gsl_histogram2d_get_yrange (const gsl_histogram2d *
1528 H, size_t J, double * YLOWER, double * YUPPER)
1529 These functions find the upper and lower range limits of the I-th
1530 and J-th bins in the x and y directions of the histogram H. The
1531 range limits are stored in XLOWER and XUPPER or YLOWER and YUPPER.
1532 The lower limits are inclusive (i.e. events with these
1533 coordinates are included in the bin) and the upper limits are
1534 exclusive (i.e. events with the value of the upper limit are not
1535 included and fall in the neighboring higher bin, if it exists).
1536 The functions return 0 to indicate success. If I or J lies
1537 outside the valid range of indices for the histogram then the
1538 error handler is called with an error code of `GSL_EDOM'.
1540 -- Function: double gsl_histogram2d_xmax (const gsl_histogram2d * H)
1541 -- Function: double gsl_histogram2d_xmin (const gsl_histogram2d * H)
1542 -- Function: size_t gsl_histogram2d_nx (const gsl_histogram2d * H)
1543 -- Function: double gsl_histogram2d_ymax (const gsl_histogram2d * H)
1544 -- Function: double gsl_histogram2d_ymin (const gsl_histogram2d * H)
1545 -- Function: size_t gsl_histogram2d_ny (const gsl_histogram2d * H)
1546 These functions return the maximum upper and minimum lower range
1547 limits and the number of bins for the x and y directions of the
1548 histogram H. They provide a way of determining these values
1549 without accessing the `gsl_histogram2d' struct directly.
1551 -- Function: void gsl_histogram2d_reset (gsl_histogram2d * H)
1552 This function resets all the bins of the histogram H to zero.
1555 File: gsl-ref.info, Node: Searching 2D histogram ranges, Next: 2D Histogram Statistics, Prev: Updating and accessing 2D histogram elements, Up: Histograms
1557 21.17 Searching 2D histogram ranges
1558 ===================================
1560 The following functions are used by the access and update routines to
1561 locate the bin which corresponds to a given (x,y) coordinate.
1563 -- Function: int gsl_histogram2d_find (const gsl_histogram2d * H,
1564 double X, double Y, size_t * I, size_t * J)
1565 This function finds and sets the indices I and J to the to the bin
1566 which covers the coordinates (X,Y). The bin is located using a
1567 binary search. The search includes an optimization for histograms
1568 with uniform ranges, and will return the correct bin immediately
1569 in this case. If (x,y) is found then the function sets the indices
1570 (I,J) and returns `GSL_SUCCESS'. If (x,y) lies outside the valid
1571 range of the histogram then the function returns `GSL_EDOM' and
1572 the error handler is invoked.
1575 File: gsl-ref.info, Node: 2D Histogram Statistics, Next: 2D Histogram Operations, Prev: Searching 2D histogram ranges, Up: Histograms
1577 21.18 2D Histogram Statistics
1578 =============================
1580 -- Function: double gsl_histogram2d_max_val (const gsl_histogram2d * H)
1581 This function returns the maximum value contained in the histogram
1584 -- Function: void gsl_histogram2d_max_bin (const gsl_histogram2d * H,
1585 size_t * I, size_t * J)
1586 This function finds the indices of the bin containing the maximum
1587 value in the histogram H and stores the result in (I,J). In the
1588 case where several bins contain the same maximum value the first
1589 bin found is returned.
1591 -- Function: double gsl_histogram2d_min_val (const gsl_histogram2d * H)
1592 This function returns the minimum value contained in the histogram
1595 -- Function: void gsl_histogram2d_min_bin (const gsl_histogram2d * H,
1596 size_t * I, size_t * J)
1597 This function finds the indices of the bin containing the minimum
1598 value in the histogram H and stores the result in (I,J). In the
1599 case where several bins contain the same maximum value the first
1600 bin found is returned.
1602 -- Function: double gsl_histogram2d_xmean (const gsl_histogram2d * H)
1603 This function returns the mean of the histogrammed x variable,
1604 where the histogram is regarded as a probability distribution.
1605 Negative bin values are ignored for the purposes of this
1608 -- Function: double gsl_histogram2d_ymean (const gsl_histogram2d * H)
1609 This function returns the mean of the histogrammed y variable,
1610 where the histogram is regarded as a probability distribution.
1611 Negative bin values are ignored for the purposes of this
1614 -- Function: double gsl_histogram2d_xsigma (const gsl_histogram2d * H)
1615 This function returns the standard deviation of the histogrammed x
1616 variable, where the histogram is regarded as a probability
1617 distribution. Negative bin values are ignored for the purposes of
1620 -- Function: double gsl_histogram2d_ysigma (const gsl_histogram2d * H)
1621 This function returns the standard deviation of the histogrammed y
1622 variable, where the histogram is regarded as a probability
1623 distribution. Negative bin values are ignored for the purposes of
1626 -- Function: double gsl_histogram2d_cov (const gsl_histogram2d * H)
1627 This function returns the covariance of the histogrammed x and y
1628 variables, where the histogram is regarded as a probability
1629 distribution. Negative bin values are ignored for the purposes of
1632 -- Function: double gsl_histogram2d_sum (const gsl_histogram2d * H)
1633 This function returns the sum of all bin values. Negative bin
1634 values are included in the sum.
1637 File: gsl-ref.info, Node: 2D Histogram Operations, Next: Reading and writing 2D histograms, Prev: 2D Histogram Statistics, Up: Histograms
1639 21.19 2D Histogram Operations
1640 =============================
1642 -- Function: int gsl_histogram2d_equal_bins_p (const gsl_histogram2d *
1643 H1, const gsl_histogram2d * H2)
1644 This function returns 1 if all the individual bin ranges of the two
1645 histograms are identical, and 0 otherwise.
1647 -- Function: int gsl_histogram2d_add (gsl_histogram2d * H1, const
1648 gsl_histogram2d * H2)
1649 This function adds the contents of the bins in histogram H2 to the
1650 corresponding bins of histogram H1, i.e. h'_1(i,j) = h_1(i,j) +
1651 h_2(i,j). The two histograms must have identical bin ranges.
1653 -- Function: int gsl_histogram2d_sub (gsl_histogram2d * H1, const
1654 gsl_histogram2d * H2)
1655 This function subtracts the contents of the bins in histogram H2
1656 from the corresponding bins of histogram H1, i.e. h'_1(i,j) =
1657 h_1(i,j) - h_2(i,j). The two histograms must have identical bin
1660 -- Function: int gsl_histogram2d_mul (gsl_histogram2d * H1, const
1661 gsl_histogram2d * H2)
1662 This function multiplies the contents of the bins of histogram H1
1663 by the contents of the corresponding bins in histogram H2, i.e.
1664 h'_1(i,j) = h_1(i,j) * h_2(i,j). The two histograms must have
1665 identical bin ranges.
1667 -- Function: int gsl_histogram2d_div (gsl_histogram2d * H1, const
1668 gsl_histogram2d * H2)
1669 This function divides the contents of the bins of histogram H1 by
1670 the contents of the corresponding bins in histogram H2, i.e.
1671 h'_1(i,j) = h_1(i,j) / h_2(i,j). The two histograms must have
1672 identical bin ranges.
1674 -- Function: int gsl_histogram2d_scale (gsl_histogram2d * H, double
1676 This function multiplies the contents of the bins of histogram H
1677 by the constant SCALE, i.e. h'_1(i,j) = h_1(i,j) scale.
1679 -- Function: int gsl_histogram2d_shift (gsl_histogram2d * H, double
1681 This function shifts the contents of the bins of histogram H by
1682 the constant OFFSET, i.e. h'_1(i,j) = h_1(i,j) + offset.
1685 File: gsl-ref.info, Node: Reading and writing 2D histograms, Next: Resampling from 2D histograms, Prev: 2D Histogram Operations, Up: Histograms
1687 21.20 Reading and writing 2D histograms
1688 =======================================
1690 The library provides functions for reading and writing two dimensional
1691 histograms to a file as binary data or formatted text.
1693 -- Function: int gsl_histogram2d_fwrite (FILE * STREAM, const
1694 gsl_histogram2d * H)
1695 This function writes the ranges and bins of the histogram H to the
1696 stream STREAM in binary format. The return value is 0 for success
1697 and `GSL_EFAILED' if there was a problem writing to the file.
1698 Since the data is written in the native binary format it may not
1699 be portable between different architectures.
1701 -- Function: int gsl_histogram2d_fread (FILE * STREAM, gsl_histogram2d
1703 This function reads into the histogram H from the stream STREAM in
1704 binary format. The histogram H must be preallocated with the
1705 correct size since the function uses the number of x and y bins in
1706 H to determine how many bytes to read. The return value is 0 for
1707 success and `GSL_EFAILED' if there was a problem reading from the
1708 file. The data is assumed to have been written in the native
1709 binary format on the same architecture.
1711 -- Function: int gsl_histogram2d_fprintf (FILE * STREAM, const
1712 gsl_histogram2d * H, const char * RANGE_FORMAT, const char *
1714 This function writes the ranges and bins of the histogram H
1715 line-by-line to the stream STREAM using the format specifiers
1716 RANGE_FORMAT and BIN_FORMAT. These should be one of the `%g',
1717 `%e' or `%f' formats for floating point numbers. The function
1718 returns 0 for success and `GSL_EFAILED' if there was a problem
1719 writing to the file. The histogram output is formatted in five
1720 columns, and the columns are separated by spaces, like this,
1722 xrange[0] xrange[1] yrange[0] yrange[1] bin(0,0)
1723 xrange[0] xrange[1] yrange[1] yrange[2] bin(0,1)
1724 xrange[0] xrange[1] yrange[2] yrange[3] bin(0,2)
1726 xrange[0] xrange[1] yrange[ny-1] yrange[ny] bin(0,ny-1)
1728 xrange[1] xrange[2] yrange[0] yrange[1] bin(1,0)
1729 xrange[1] xrange[2] yrange[1] yrange[2] bin(1,1)
1730 xrange[1] xrange[2] yrange[1] yrange[2] bin(1,2)
1732 xrange[1] xrange[2] yrange[ny-1] yrange[ny] bin(1,ny-1)
1736 xrange[nx-1] xrange[nx] yrange[0] yrange[1] bin(nx-1,0)
1737 xrange[nx-1] xrange[nx] yrange[1] yrange[2] bin(nx-1,1)
1738 xrange[nx-1] xrange[nx] yrange[1] yrange[2] bin(nx-1,2)
1740 xrange[nx-1] xrange[nx] yrange[ny-1] yrange[ny] bin(nx-1,ny-1)
1742 Each line contains the lower and upper limits of the bin and the
1743 contents of the bin. Since the upper limits of the each bin are
1744 the lower limits of the neighboring bins there is duplication of
1745 these values but this allows the histogram to be manipulated with
1746 line-oriented tools.
1748 -- Function: int gsl_histogram2d_fscanf (FILE * STREAM,
1749 gsl_histogram2d * H)
1750 This function reads formatted data from the stream STREAM into the
1751 histogram H. The data is assumed to be in the five-column format
1752 used by `gsl_histogram2d_fprintf'. The histogram H must be
1753 preallocated with the correct lengths since the function uses the
1754 sizes of H to determine how many numbers to read. The function
1755 returns 0 for success and `GSL_EFAILED' if there was a problem
1756 reading from the file.
1759 File: gsl-ref.info, Node: Resampling from 2D histograms, Next: Example programs for 2D histograms, Prev: Reading and writing 2D histograms, Up: Histograms
1761 21.21 Resampling from 2D histograms
1762 ===================================
1764 As in the one-dimensional case, a two-dimensional histogram made by
1765 counting events can be regarded as a measurement of a probability
1766 distribution. Allowing for statistical error, the height of each bin
1767 represents the probability of an event where (x,y) falls in the range
1768 of that bin. For a two-dimensional histogram the probability
1769 distribution takes the form p(x,y) dx dy where,
1771 p(x,y) = n_{ij}/ (N A_{ij})
1773 In this equation n_{ij} is the number of events in the bin which
1774 contains (x,y), A_{ij} is the area of the bin and N is the total number
1775 of events. The distribution of events within each bin is assumed to be
1778 -- Data Type: gsl_histogram2d_pdf
1780 This is the number of histogram bins used to approximate the
1781 probability distribution function in the x and y directions.
1784 The ranges of the bins in the x-direction are stored in an
1785 array of NX + 1 elements pointed to by XRANGE.
1788 The ranges of the bins in the y-direction are stored in an
1789 array of NY + 1 pointed to by YRANGE.
1792 The cumulative probability for the bins is stored in an array
1793 of NX*NY elements pointed to by SUM.
1795 The following functions allow you to create a `gsl_histogram2d_pdf'
1796 struct which represents a two dimensional probability distribution and
1797 generate random samples from it.
1799 -- Function: gsl_histogram2d_pdf * gsl_histogram2d_pdf_alloc (size_t
1801 This function allocates memory for a two-dimensional probability
1802 distribution of size NX-by-NY and returns a pointer to a newly
1803 initialized `gsl_histogram2d_pdf' struct. If insufficient memory
1804 is available a null pointer is returned and the error handler is
1805 invoked with an error code of `GSL_ENOMEM'.
1807 -- Function: int gsl_histogram2d_pdf_init (gsl_histogram2d_pdf * P,
1808 const gsl_histogram2d * H)
1809 This function initializes the two-dimensional probability
1810 distribution calculated P from the histogram H. If any of the
1811 bins of H are negative then the error handler is invoked with an
1812 error code of `GSL_EDOM' because a probability distribution cannot
1813 contain negative values.
1815 -- Function: void gsl_histogram2d_pdf_free (gsl_histogram2d_pdf * P)
1816 This function frees the two-dimensional probability distribution
1817 function P and all of the memory associated with it.
1819 -- Function: int gsl_histogram2d_pdf_sample (const gsl_histogram2d_pdf
1820 * P, double R1, double R2, double * X, double * Y)
1821 This function uses two uniform random numbers between zero and one,
1822 R1 and R2, to compute a single random sample from the
1823 two-dimensional probability distribution P.
1826 File: gsl-ref.info, Node: Example programs for 2D histograms, Prev: Resampling from 2D histograms, Up: Histograms
1828 21.22 Example programs for 2D histograms
1829 ========================================
1831 This program demonstrates two features of two-dimensional histograms.
1832 First a 10-by-10 two-dimensional histogram is created with x and y
1833 running from 0 to 1. Then a few sample points are added to the
1834 histogram, at (0.3,0.3) with a height of 1, at (0.8,0.1) with a height
1835 of 5 and at (0.7,0.9) with a height of 0.5. This histogram with three
1836 events is used to generate a random sample of 1000 simulated events,
1837 which are printed out.
1840 #include <gsl/gsl_rng.h>
1841 #include <gsl/gsl_histogram2d.h>
1846 const gsl_rng_type * T;
1849 gsl_histogram2d * h = gsl_histogram2d_alloc (10, 10);
1851 gsl_histogram2d_set_ranges_uniform (h,
1855 gsl_histogram2d_accumulate (h, 0.3, 0.3, 1);
1856 gsl_histogram2d_accumulate (h, 0.8, 0.1, 5);
1857 gsl_histogram2d_accumulate (h, 0.7, 0.9, 0.5);
1859 gsl_rng_env_setup ();
1861 T = gsl_rng_default;
1862 r = gsl_rng_alloc (T);
1866 gsl_histogram2d_pdf * p
1867 = gsl_histogram2d_pdf_alloc (h->nx, h->ny);
1869 gsl_histogram2d_pdf_init (p, h);
1871 for (i = 0; i < 1000; i++) {
1873 double u = gsl_rng_uniform (r);
1874 double v = gsl_rng_uniform (r);
1876 gsl_histogram2d_pdf_sample (p, u, v, &x, &y);
1878 printf ("%g %g\n", x, y);
1881 gsl_histogram2d_pdf_free (p);
1884 gsl_histogram2d_free (h);
1892 File: gsl-ref.info, Node: N-tuples, Next: Monte Carlo Integration, Prev: Histograms, Up: Top
1897 This chapter describes functions for creating and manipulating
1898 "ntuples", sets of values associated with events. The ntuples are
1899 stored in files. Their values can be extracted in any combination and
1900 "booked" in a histogram using a selection function.
1902 The values to be stored are held in a user-defined data structure,
1903 and an ntuple is created associating this data structure with a file.
1904 The values are then written to the file (normally inside a loop) using
1905 the ntuple functions described below.
1907 A histogram can be created from ntuple data by providing a selection
1908 function and a value function. The selection function specifies whether
1909 an event should be included in the subset to be analyzed or not. The
1910 value function computes the entry to be added to the histogram for each
1913 All the ntuple functions are defined in the header file
1918 * The ntuple struct::
1919 * Creating ntuples::
1920 * Opening an existing ntuple file::
1922 * Reading ntuples ::
1923 * Closing an ntuple file::
1924 * Histogramming ntuple values::
1925 * Example ntuple programs::
1926 * Ntuple References and Further Reading::
1929 File: gsl-ref.info, Node: The ntuple struct, Next: Creating ntuples, Up: N-tuples
1931 22.1 The ntuple struct
1932 ======================
1934 Ntuples are manipulated using the `gsl_ntuple' struct. This struct
1935 contains information on the file where the ntuple data is stored, a
1936 pointer to the current ntuple data row and the size of the user-defined
1946 File: gsl-ref.info, Node: Creating ntuples, Next: Opening an existing ntuple file, Prev: The ntuple struct, Up: N-tuples
1948 22.2 Creating ntuples
1949 =====================
1951 -- Function: gsl_ntuple * gsl_ntuple_create (char * FILENAME, void *
1952 NTUPLE_DATA, size_t SIZE)
1953 This function creates a new write-only ntuple file FILENAME for
1954 ntuples of size SIZE and returns a pointer to the newly created
1955 ntuple struct. Any existing file with the same name is truncated
1956 to zero length and overwritten. A pointer to memory for the
1957 current ntuple row NTUPLE_DATA must be supplied--this is used to
1958 copy ntuples in and out of the file.
1961 File: gsl-ref.info, Node: Opening an existing ntuple file, Next: Writing ntuples, Prev: Creating ntuples, Up: N-tuples
1963 22.3 Opening an existing ntuple file
1964 ====================================
1966 -- Function: gsl_ntuple * gsl_ntuple_open (char * FILENAME, void *
1967 NTUPLE_DATA, size_t SIZE)
1968 This function opens an existing ntuple file FILENAME for reading
1969 and returns a pointer to a corresponding ntuple struct. The
1970 ntuples in the file must have size SIZE. A pointer to memory for
1971 the current ntuple row NTUPLE_DATA must be supplied--this is used
1972 to copy ntuples in and out of the file.
1975 File: gsl-ref.info, Node: Writing ntuples, Next: Reading ntuples, Prev: Opening an existing ntuple file, Up: N-tuples
1977 22.4 Writing ntuples
1978 ====================
1980 -- Function: int gsl_ntuple_write (gsl_ntuple * NTUPLE)
1981 This function writes the current ntuple NTUPLE->NTUPLE_DATA of
1982 size NTUPLE->SIZE to the corresponding file.
1984 -- Function: int gsl_ntuple_bookdata (gsl_ntuple * NTUPLE)
1985 This function is a synonym for `gsl_ntuple_write'.
1988 File: gsl-ref.info, Node: Reading ntuples, Next: Closing an ntuple file, Prev: Writing ntuples, Up: N-tuples
1990 22.5 Reading ntuples
1991 ====================
1993 -- Function: int gsl_ntuple_read (gsl_ntuple * NTUPLE)
1994 This function reads the current row of the ntuple file for NTUPLE
1995 and stores the values in NTUPLE->DATA.
1998 File: gsl-ref.info, Node: Closing an ntuple file, Next: Histogramming ntuple values, Prev: Reading ntuples, Up: N-tuples
2000 22.6 Closing an ntuple file
2001 ===========================
2003 -- Function: int gsl_ntuple_close (gsl_ntuple * NTUPLE)
2004 This function closes the ntuple file NTUPLE and frees its
2005 associated allocated memory.
2008 File: gsl-ref.info, Node: Histogramming ntuple values, Next: Example ntuple programs, Prev: Closing an ntuple file, Up: N-tuples
2010 22.7 Histogramming ntuple values
2011 ================================
2013 Once an ntuple has been created its contents can be histogrammed in
2014 various ways using the function `gsl_ntuple_project'. Two user-defined
2015 functions must be provided, a function to select events and a function
2016 to compute scalar values. The selection function and the value function
2017 both accept the ntuple row as a first argument and other parameters as
2020 The "selection function" determines which ntuple rows are selected
2021 for histogramming. It is defined by the following struct,
2023 int (* function) (void * ntuple_data, void * params);
2025 } gsl_ntuple_select_fn;
2027 The struct component FUNCTION should return a non-zero value for each
2028 ntuple row that is to be included in the histogram.
2030 The "value function" computes scalar values for those ntuple rows
2031 selected by the selection function,
2033 double (* function) (void * ntuple_data, void * params);
2035 } gsl_ntuple_value_fn;
2037 In this case the struct component FUNCTION should return the value to
2038 be added to the histogram for the ntuple row.
2040 -- Function: int gsl_ntuple_project (gsl_histogram * H, gsl_ntuple *
2041 NTUPLE, gsl_ntuple_value_fn * VALUE_FUNC,
2042 gsl_ntuple_select_fn * SELECT_FUNC)
2043 This function updates the histogram H from the ntuple NTUPLE using
2044 the functions VALUE_FUNC and SELECT_FUNC. For each ntuple row
2045 where the selection function SELECT_FUNC is non-zero the
2046 corresponding value of that row is computed using the function
2047 VALUE_FUNC and added to the histogram. Those ntuple rows where
2048 SELECT_FUNC returns zero are ignored. New entries are added to
2049 the histogram, so subsequent calls can be used to accumulate
2050 further data in the same histogram.
2053 File: gsl-ref.info, Node: Example ntuple programs, Next: Ntuple References and Further Reading, Prev: Histogramming ntuple values, Up: N-tuples
2058 The following example programs demonstrate the use of ntuples in
2059 managing a large dataset. The first program creates a set of 10,000
2060 simulated "events", each with 3 associated values (x,y,z). These are
2061 generated from a gaussian distribution with unit variance, for
2062 demonstration purposes, and written to the ntuple file `test.dat'.
2064 #include <gsl/gsl_ntuple.h>
2065 #include <gsl/gsl_rng.h>
2066 #include <gsl/gsl_randist.h>
2078 const gsl_rng_type * T;
2081 struct data ntuple_row;
2085 = gsl_ntuple_create ("test.dat", &ntuple_row,
2086 sizeof (ntuple_row));
2088 gsl_rng_env_setup ();
2090 T = gsl_rng_default;
2091 r = gsl_rng_alloc (T);
2093 for (i = 0; i < 10000; i++)
2095 ntuple_row.x = gsl_ran_ugaussian (r);
2096 ntuple_row.y = gsl_ran_ugaussian (r);
2097 ntuple_row.z = gsl_ran_ugaussian (r);
2099 gsl_ntuple_write (ntuple);
2102 gsl_ntuple_close (ntuple);
2108 The next program analyses the ntuple data in the file `test.dat'. The
2109 analysis procedure is to compute the squared-magnitude of each event,
2110 E^2=x^2+y^2+z^2, and select only those which exceed a lower limit of
2111 1.5. The selected events are then histogrammed using their E^2 values.
2114 #include <gsl/gsl_ntuple.h>
2115 #include <gsl/gsl_histogram.h>
2124 int sel_func (void *ntuple_data, void *params);
2125 double val_func (void *ntuple_data, void *params);
2130 struct data ntuple_row;
2133 = gsl_ntuple_open ("test.dat", &ntuple_row,
2134 sizeof (ntuple_row));
2137 gsl_ntuple_select_fn S;
2138 gsl_ntuple_value_fn V;
2140 gsl_histogram *h = gsl_histogram_alloc (100);
2141 gsl_histogram_set_ranges_uniform(h, 0.0, 10.0);
2143 S.function = &sel_func;
2146 V.function = &val_func;
2149 gsl_ntuple_project (h, ntuple, &V, &S);
2150 gsl_histogram_fprintf (stdout, h, "%f", "%f");
2151 gsl_histogram_free (h);
2152 gsl_ntuple_close (ntuple);
2158 sel_func (void *ntuple_data, void *params)
2160 struct data * data = (struct data *) ntuple_data;
2161 double x, y, z, E2, scale;
2162 scale = *(double *) params;
2168 E2 = x * x + y * y + z * z;
2174 val_func (void *ntuple_data, void *params)
2176 struct data * data = (struct data *) ntuple_data;
2183 return x * x + y * y + z * z;
2186 The following plot shows the distribution of the selected events.
2187 Note the cut-off at the lower bound.
2190 File: gsl-ref.info, Node: Ntuple References and Further Reading, Prev: Example ntuple programs, Up: N-tuples
2192 22.9 References and Further Reading
2193 ===================================
2195 Further information on the use of ntuples can be found in the
2196 documentation for the CERN packages PAW and HBOOK (available online).
2199 File: gsl-ref.info, Node: Monte Carlo Integration, Next: Simulated Annealing, Prev: N-tuples, Up: Top
2201 23 Monte Carlo Integration
2202 **************************
2204 This chapter describes routines for multidimensional Monte Carlo
2205 integration. These include the traditional Monte Carlo method and
2206 adaptive algorithms such as VEGAS and MISER which use importance
2207 sampling and stratified sampling techniques. Each algorithm computes an
2208 estimate of a multidimensional definite integral of the form,
2210 I = \int_xl^xu dx \int_yl^yu dy ... f(x, y, ...)
2212 over a hypercubic region ((x_l,x_u), (y_l,y_u), ...) using a fixed
2213 number of function calls. The routines also provide a statistical
2214 estimate of the error on the result. This error estimate should be
2215 taken as a guide rather than as a strict error bound--random sampling
2216 of the region may not uncover all the important features of the
2217 function, resulting in an underestimate of the error.
2219 The functions are defined in separate header files for each routine,
2220 `gsl_monte_plain.h', `gsl_monte_miser.h' and `gsl_monte_vegas.h'.
2224 * Monte Carlo Interface::
2225 * PLAIN Monte Carlo::
2228 * Monte Carlo Examples::
2229 * Monte Carlo Integration References and Further Reading::
2232 File: gsl-ref.info, Node: Monte Carlo Interface, Next: PLAIN Monte Carlo, Up: Monte Carlo Integration
2237 All of the Monte Carlo integration routines use the same general form of
2238 interface. There is an allocator to allocate memory for control
2239 variables and workspace, a routine to initialize those control
2240 variables, the integrator itself, and a function to free the space when
2243 Each integration function requires a random number generator to be
2244 supplied, and returns an estimate of the integral and its standard
2245 deviation. The accuracy of the result is determined by the number of
2246 function calls specified by the user. If a known level of accuracy is
2247 required this can be achieved by calling the integrator several times
2248 and averaging the individual results until the desired accuracy is
2251 Random sample points used within the Monte Carlo routines are always
2252 chosen strictly within the integration region, so that endpoint
2253 singularities are automatically avoided.
2255 The function to be integrated has its own datatype, defined in the
2256 header file `gsl_monte.h'.
2258 -- Data Type: gsl_monte_function
2259 This data type defines a general function with parameters for Monte
2262 `double (* f) (double * X, size_t DIM, void * PARAMS)'
2263 this function should return the value f(x,params) for the
2264 argument X and parameters PARAMS, where X is an array of size
2265 DIM giving the coordinates of the point where the function is
2269 the number of dimensions for X.
2272 a pointer to the parameters of the function.
2274 Here is an example for a quadratic function in two dimensions,
2276 f(x,y) = a x^2 + b x y + c y^2
2278 with a = 3, b = 2, c = 1. The following code defines a
2279 `gsl_monte_function' `F' which you could pass to an integrator:
2281 struct my_f_params { double a; double b; double c; };
2284 my_f (double x[], size_t dim, void * p) {
2285 struct my_f_params * fp = (struct my_f_params *)p;
2289 fprintf (stderr, "error: dim != 2");
2293 return fp->a * x[0] * x[0]
2294 + fp->b * x[0] * x[1]
2295 + fp->c * x[1] * x[1];
2298 gsl_monte_function F;
2299 struct my_f_params params = { 3.0, 2.0, 1.0 };
2305 The function f(x) can be evaluated using the following macro,
2307 #define GSL_MONTE_FN_EVAL(F,x)
2308 (*((F)->f))(x,(F)->dim,(F)->params)
2311 File: gsl-ref.info, Node: PLAIN Monte Carlo, Next: MISER, Prev: Monte Carlo Interface, Up: Monte Carlo Integration
2313 23.2 PLAIN Monte Carlo
2314 ======================
2316 The plain Monte Carlo algorithm samples points randomly from the
2317 integration region to estimate the integral and its error. Using this
2318 algorithm the estimate of the integral E(f; N) for N randomly
2319 distributed points x_i is given by,
2321 E(f; N) = = V <f> = (V / N) \sum_i^N f(x_i)
2323 where V is the volume of the integration region. The error on this
2324 estimate \sigma(E;N) is calculated from the estimated variance of the
2327 \sigma^2 (E; N) = (V / N) \sum_i^N (f(x_i) - <f>)^2.
2329 For large N this variance decreases asymptotically as \Var(f)/N, where
2330 \Var(f) is the true variance of the function over the integration
2331 region. The error estimate itself should decrease as
2332 \sigma(f)/\sqrt{N}. The familiar law of errors decreasing as
2333 1/\sqrt{N} applies--to reduce the error by a factor of 10 requires a
2334 100-fold increase in the number of sample points.
2336 The functions described in this section are declared in the header
2337 file `gsl_monte_plain.h'.
2339 -- Function: gsl_monte_plain_state * gsl_monte_plain_alloc (size_t DIM)
2340 This function allocates and initializes a workspace for Monte Carlo
2341 integration in DIM dimensions.
2343 -- Function: int gsl_monte_plain_init (gsl_monte_plain_state* S)
2344 This function initializes a previously allocated integration state.
2345 This allows an existing workspace to be reused for different
2348 -- Function: int gsl_monte_plain_integrate (gsl_monte_function * F,
2349 double * XL, double * XU, size_t DIM, size_t CALLS, gsl_rng *
2350 R, gsl_monte_plain_state * S, double * RESULT, double *
2352 This routines uses the plain Monte Carlo algorithm to integrate the
2353 function F over the DIM-dimensional hypercubic region defined by
2354 the lower and upper limits in the arrays XL and XU, each of size
2355 DIM. The integration uses a fixed number of function calls CALLS,
2356 and obtains random sampling points using the random number
2357 generator R. A previously allocated workspace S must be supplied.
2358 The result of the integration is returned in RESULT, with an
2359 estimated absolute error ABSERR.
2361 -- Function: void gsl_monte_plain_free (gsl_monte_plain_state * S)
2362 This function frees the memory associated with the integrator state
2366 File: gsl-ref.info, Node: MISER, Next: VEGAS, Prev: PLAIN Monte Carlo, Up: Monte Carlo Integration
2371 The MISER algorithm of Press and Farrar is based on recursive
2372 stratified sampling. This technique aims to reduce the overall
2373 integration error by concentrating integration points in the regions of
2376 The idea of stratified sampling begins with the observation that for
2377 two disjoint regions a and b with Monte Carlo estimates of the integral
2378 E_a(f) and E_b(f) and variances \sigma_a^2(f) and \sigma_b^2(f), the
2379 variance \Var(f) of the combined estimate E(f) = (1/2) (E_a(f) + E_b(f))
2382 \Var(f) = (\sigma_a^2(f) / 4 N_a) + (\sigma_b^2(f) / 4 N_b).
2384 It can be shown that this variance is minimized by distributing the
2387 N_a / (N_a + N_b) = \sigma_a / (\sigma_a + \sigma_b).
2389 Hence the smallest error estimate is obtained by allocating sample
2390 points in proportion to the standard deviation of the function in each
2393 The MISER algorithm proceeds by bisecting the integration region
2394 along one coordinate axis to give two sub-regions at each step. The
2395 direction is chosen by examining all d possible bisections and
2396 selecting the one which will minimize the combined variance of the two
2397 sub-regions. The variance in the sub-regions is estimated by sampling
2398 with a fraction of the total number of points available to the current
2399 step. The same procedure is then repeated recursively for each of the
2400 two half-spaces from the best bisection. The remaining sample points are
2401 allocated to the sub-regions using the formula for N_a and N_b. This
2402 recursive allocation of integration points continues down to a
2403 user-specified depth where each sub-region is integrated using a plain
2404 Monte Carlo estimate. These individual values and their error
2405 estimates are then combined upwards to give an overall result and an
2406 estimate of its error.
2408 The functions described in this section are declared in the header
2409 file `gsl_monte_miser.h'.
2411 -- Function: gsl_monte_miser_state * gsl_monte_miser_alloc (size_t DIM)
2412 This function allocates and initializes a workspace for Monte Carlo
2413 integration in DIM dimensions. The workspace is used to maintain
2414 the state of the integration.
2416 -- Function: int gsl_monte_miser_init (gsl_monte_miser_state* S)
2417 This function initializes a previously allocated integration state.
2418 This allows an existing workspace to be reused for different
2421 -- Function: int gsl_monte_miser_integrate (gsl_monte_function * F,
2422 double * XL, double * XU, size_t DIM, size_t CALLS, gsl_rng *
2423 R, gsl_monte_miser_state * S, double * RESULT, double *
2425 This routines uses the MISER Monte Carlo algorithm to integrate the
2426 function F over the DIM-dimensional hypercubic region defined by
2427 the lower and upper limits in the arrays XL and XU, each of size
2428 DIM. The integration uses a fixed number of function calls CALLS,
2429 and obtains random sampling points using the random number
2430 generator R. A previously allocated workspace S must be supplied.
2431 The result of the integration is returned in RESULT, with an
2432 estimated absolute error ABSERR.
2434 -- Function: void gsl_monte_miser_free (gsl_monte_miser_state * S)
2435 This function frees the memory associated with the integrator state
2438 The MISER algorithm has several configurable parameters. The
2439 following variables can be accessed through the `gsl_monte_miser_state'
2442 -- Variable: double estimate_frac
2443 This parameter specifies the fraction of the currently available
2444 number of function calls which are allocated to estimating the
2445 variance at each recursive step. The default value is 0.1.
2447 -- Variable: size_t min_calls
2448 This parameter specifies the minimum number of function calls
2449 required for each estimate of the variance. If the number of
2450 function calls allocated to the estimate using ESTIMATE_FRAC falls
2451 below MIN_CALLS then MIN_CALLS are used instead. This ensures
2452 that each estimate maintains a reasonable level of accuracy. The
2453 default value of MIN_CALLS is `16 * dim'.
2455 -- Variable: size_t min_calls_per_bisection
2456 This parameter specifies the minimum number of function calls
2457 required to proceed with a bisection step. When a recursive step
2458 has fewer calls available than MIN_CALLS_PER_BISECTION it performs
2459 a plain Monte Carlo estimate of the current sub-region and
2460 terminates its branch of the recursion. The default value of this
2461 parameter is `32 * min_calls'.
2463 -- Variable: double alpha
2464 This parameter controls how the estimated variances for the two
2465 sub-regions of a bisection are combined when allocating points.
2466 With recursive sampling the overall variance should scale better
2467 than 1/N, since the values from the sub-regions will be obtained
2468 using a procedure which explicitly minimizes their variance. To
2469 accommodate this behavior the MISER algorithm allows the total
2470 variance to depend on a scaling parameter \alpha,
2472 \Var(f) = {\sigma_a \over N_a^\alpha} + {\sigma_b \over N_b^\alpha}.
2474 The authors of the original paper describing MISER recommend the
2475 value \alpha = 2 as a good choice, obtained from numerical
2476 experiments, and this is used as the default value in this
2479 -- Variable: double dither
2480 This parameter introduces a random fractional variation of size
2481 DITHER into each bisection, which can be used to break the
2482 symmetry of integrands which are concentrated near the exact
2483 center of the hypercubic integration region. The default value of
2484 dither is zero, so no variation is introduced. If needed, a
2485 typical value of DITHER is 0.1.
2488 File: gsl-ref.info, Node: VEGAS, Next: Monte Carlo Examples, Prev: MISER, Up: Monte Carlo Integration
2493 The VEGAS algorithm of Lepage is based on importance sampling. It
2494 samples points from the probability distribution described by the
2495 function |f|, so that the points are concentrated in the regions that
2496 make the largest contribution to the integral.
2498 In general, if the Monte Carlo integral of f is sampled with points
2499 distributed according to a probability distribution described by the
2500 function g, we obtain an estimate E_g(f; N),
2502 E_g(f; N) = E(f/g; N)
2504 with a corresponding variance,
2506 \Var_g(f; N) = \Var(f/g; N).
2508 If the probability distribution is chosen as g = |f|/I(|f|) then it can
2509 be shown that the variance V_g(f; N) vanishes, and the error in the
2510 estimate will be zero. In practice it is not possible to sample from
2511 the exact distribution g for an arbitrary function, so importance
2512 sampling algorithms aim to produce efficient approximations to the
2513 desired distribution.
2515 The VEGAS algorithm approximates the exact distribution by making a
2516 number of passes over the integration region while histogramming the
2517 function f. Each histogram is used to define a sampling distribution
2518 for the next pass. Asymptotically this procedure converges to the
2519 desired distribution. In order to avoid the number of histogram bins
2520 growing like K^d the probability distribution is approximated by a
2521 separable function: g(x_1, x_2, ...) = g_1(x_1) g_2(x_2) ... so that
2522 the number of bins required is only Kd. This is equivalent to locating
2523 the peaks of the function from the projections of the integrand onto
2524 the coordinate axes. The efficiency of VEGAS depends on the validity
2525 of this assumption. It is most efficient when the peaks of the
2526 integrand are well-localized. If an integrand can be rewritten in a
2527 form which is approximately separable this will increase the efficiency
2528 of integration with VEGAS.
2530 VEGAS incorporates a number of additional features, and combines both
2531 stratified sampling and importance sampling. The integration region is
2532 divided into a number of "boxes", with each box getting a fixed number
2533 of points (the goal is 2). Each box can then have a fractional number
2534 of bins, but if the ratio of bins-per-box is less than two, Vegas
2535 switches to a kind variance reduction (rather than importance sampling).
2537 -- Function: gsl_monte_vegas_state * gsl_monte_vegas_alloc (size_t DIM)
2538 This function allocates and initializes a workspace for Monte Carlo
2539 integration in DIM dimensions. The workspace is used to maintain
2540 the state of the integration.
2542 -- Function: int gsl_monte_vegas_init (gsl_monte_vegas_state* S)
2543 This function initializes a previously allocated integration state.
2544 This allows an existing workspace to be reused for different
2547 -- Function: int gsl_monte_vegas_integrate (gsl_monte_function * F,
2548 double * XL, double * XU, size_t DIM, size_t CALLS, gsl_rng *
2549 R, gsl_monte_vegas_state * S, double * RESULT, double *
2551 This routines uses the VEGAS Monte Carlo algorithm to integrate the
2552 function F over the DIM-dimensional hypercubic region defined by
2553 the lower and upper limits in the arrays XL and XU, each of size
2554 DIM. The integration uses a fixed number of function calls CALLS,
2555 and obtains random sampling points using the random number
2556 generator R. A previously allocated workspace S must be supplied.
2557 The result of the integration is returned in RESULT, with an
2558 estimated absolute error ABSERR. The result and its error
2559 estimate are based on a weighted average of independent samples.
2560 The chi-squared per degree of freedom for the weighted average is
2561 returned via the state struct component, S->CHISQ, and must be
2562 consistent with 1 for the weighted average to be reliable.
2564 -- Function: void gsl_monte_vegas_free (gsl_monte_vegas_state * S)
2565 This function frees the memory associated with the integrator state
2568 The VEGAS algorithm computes a number of independent estimates of the
2569 integral internally, according to the `iterations' parameter described
2570 below, and returns their weighted average. Random sampling of the
2571 integrand can occasionally produce an estimate where the error is zero,
2572 particularly if the function is constant in some regions. An estimate
2573 with zero error causes the weighted average to break down and must be
2574 handled separately. In the original Fortran implementations of VEGAS
2575 the error estimate is made non-zero by substituting a small value
2576 (typically `1e-30'). The implementation in GSL differs from this and
2577 avoids the use of an arbitrary constant--it either assigns the value a
2578 weight which is the average weight of the preceding estimates or
2579 discards it according to the following procedure,
2581 current estimate has zero error, weighted average has finite error
2582 The current estimate is assigned a weight which is the average
2583 weight of the preceding estimates.
2585 current estimate has finite error, previous estimates had zero error
2586 The previous estimates are discarded and the weighted averaging
2587 procedure begins with the current estimate.
2589 current estimate has zero error, previous estimates had zero error
2590 The estimates are averaged using the arithmetic mean, but no error
2593 The VEGAS algorithm is highly configurable. The following variables
2594 can be accessed through the `gsl_monte_vegas_state' struct,
2596 -- Variable: double result
2597 -- Variable: double sigma
2598 These parameters contain the raw value of the integral RESULT and
2599 its error SIGMA from the last iteration of the algorithm.
2601 -- Variable: double chisq
2602 This parameter gives the chi-squared per degree of freedom for the
2603 weighted estimate of the integral. The value of CHISQ should be
2604 close to 1. A value of CHISQ which differs significantly from 1
2605 indicates that the values from different iterations are
2606 inconsistent. In this case the weighted error will be
2607 under-estimated, and further iterations of the algorithm are
2608 needed to obtain reliable results.
2610 -- Variable: double alpha
2611 The parameter `alpha' controls the stiffness of the rebinning
2612 algorithm. It is typically set between one and two. A value of
2613 zero prevents rebinning of the grid. The default value is 1.5.
2615 -- Variable: size_t iterations
2616 The number of iterations to perform for each call to the routine.
2617 The default value is 5 iterations.
2619 -- Variable: int stage
2620 Setting this determines the "stage" of the calculation. Normally,
2621 `stage = 0' which begins with a new uniform grid and empty weighted
2622 average. Calling vegas with `stage = 1' retains the grid from the
2623 previous run but discards the weighted average, so that one can
2624 "tune" the grid using a relatively small number of points and then
2625 do a large run with `stage = 1' on the optimized grid. Setting
2626 `stage = 2' keeps the grid and the weighted average from the
2627 previous run, but may increase (or decrease) the number of
2628 histogram bins in the grid depending on the number of calls
2629 available. Choosing `stage = 3' enters at the main loop, so that
2630 nothing is changed, and is equivalent to performing additional
2631 iterations in a previous call.
2633 -- Variable: int mode
2634 The possible choices are `GSL_VEGAS_MODE_IMPORTANCE',
2635 `GSL_VEGAS_MODE_STRATIFIED', `GSL_VEGAS_MODE_IMPORTANCE_ONLY'.
2636 This determines whether VEGAS will use importance sampling or
2637 stratified sampling, or whether it can pick on its own. In low
2638 dimensions VEGAS uses strict stratified sampling (more precisely,
2639 stratified sampling is chosen if there are fewer than 2 bins per
2642 -- Variable: int verbose
2643 -- Variable: FILE * ostream
2644 These parameters set the level of information printed by VEGAS. All
2645 information is written to the stream OSTREAM. The default setting
2646 of VERBOSE is `-1', which turns off all output. A VERBOSE value
2647 of `0' prints summary information about the weighted average and
2648 final result, while a value of `1' also displays the grid
2649 coordinates. A value of `2' prints information from the rebinning
2650 procedure for each iteration.
2653 File: gsl-ref.info, Node: Monte Carlo Examples, Next: Monte Carlo Integration References and Further Reading, Prev: VEGAS, Up: Monte Carlo Integration
2658 The example program below uses the Monte Carlo routines to estimate the
2659 value of the following 3-dimensional integral from the theory of random
2662 I = \int_{-pi}^{+pi} {dk_x/(2 pi)}
2663 \int_{-pi}^{+pi} {dk_y/(2 pi)}
2664 \int_{-pi}^{+pi} {dk_z/(2 pi)}
2665 1 / (1 - cos(k_x)cos(k_y)cos(k_z)).
2667 The analytic value of this integral can be shown to be I =
2668 \Gamma(1/4)^4/(4 \pi^3) = 1.393203929685676859.... The integral gives
2669 the mean time spent at the origin by a random walk on a body-centered
2670 cubic lattice in three dimensions.
2672 For simplicity we will compute the integral over the region (0,0,0)
2673 to (\pi,\pi,\pi) and multiply by 8 to obtain the full result. The
2674 integral is slowly varying in the middle of the region but has
2675 integrable singularities at the corners (0,0,0), (0,\pi,\pi),
2676 (\pi,0,\pi) and (\pi,\pi,0). The Monte Carlo routines only select
2677 points which are strictly within the integration region and so no
2678 special measures are needed to avoid these singularities.
2681 #include <gsl/gsl_math.h>
2682 #include <gsl/gsl_monte.h>
2683 #include <gsl/gsl_monte_plain.h>
2684 #include <gsl/gsl_monte_miser.h>
2685 #include <gsl/gsl_monte_vegas.h>
2687 /* Computation of the integral,
2689 I = int (dx dy dz)/(2pi)^3 1/(1-cos(x)cos(y)cos(z))
2691 over (-pi,-pi,-pi) to (+pi, +pi, +pi). The exact answer
2692 is Gamma(1/4)^4/(4 pi^3). This example is taken from
2693 C.Itzykson, J.M.Drouffe, "Statistical Field Theory -
2694 Volume 1", Section 1.1, p21, which cites the original
2695 paper M.L.Glasser, I.J.Zucker, Proc.Natl.Acad.Sci.USA 74
2698 /* For simplicity we compute the integral over the region
2699 (0,0,0) -> (pi,pi,pi) and multiply by 8 */
2701 double exact = 1.3932039296856768591842462603255;
2704 g (double *k, size_t dim, void *params)
2706 double A = 1.0 / (M_PI * M_PI * M_PI);
2707 return A / (1.0 - cos (k[0]) * cos (k[1]) * cos (k[2]));
2711 display_results (char *title, double result, double error)
2713 printf ("%s ==================\n", title);
2714 printf ("result = % .6f\n", result);
2715 printf ("sigma = % .6f\n", error);
2716 printf ("exact = % .6f\n", exact);
2717 printf ("error = % .6f = %.1g sigma\n", result - exact,
2718 fabs (result - exact) / error);
2726 double xl[3] = { 0, 0, 0 };
2727 double xu[3] = { M_PI, M_PI, M_PI };
2729 const gsl_rng_type *T;
2732 gsl_monte_function G = { &g, 3, 0 };
2734 size_t calls = 500000;
2736 gsl_rng_env_setup ();
2738 T = gsl_rng_default;
2739 r = gsl_rng_alloc (T);
2742 gsl_monte_plain_state *s = gsl_monte_plain_alloc (3);
2743 gsl_monte_plain_integrate (&G, xl, xu, 3, calls, r, s,
2745 gsl_monte_plain_free (s);
2747 display_results ("plain", res, err);
2751 gsl_monte_miser_state *s = gsl_monte_miser_alloc (3);
2752 gsl_monte_miser_integrate (&G, xl, xu, 3, calls, r, s,
2754 gsl_monte_miser_free (s);
2756 display_results ("miser", res, err);
2760 gsl_monte_vegas_state *s = gsl_monte_vegas_alloc (3);
2762 gsl_monte_vegas_integrate (&G, xl, xu, 3, 10000, r, s,
2764 display_results ("vegas warm-up", res, err);
2766 printf ("converging...\n");
2770 gsl_monte_vegas_integrate (&G, xl, xu, 3, calls/5, r, s,
2772 printf ("result = % .6f sigma = % .6f "
2773 "chisq/dof = %.1f\n", res, err, s->chisq);
2775 while (fabs (s->chisq - 1.0) > 0.5);
2777 display_results ("vegas final", res, err);
2779 gsl_monte_vegas_free (s);
2787 With 500,000 function calls the plain Monte Carlo algorithm achieves a
2788 fractional error of 0.6%. The estimated error `sigma' is consistent
2789 with the actual error, and the computed result differs from the true
2790 result by about one standard deviation,
2792 plain ==================
2796 error = -0.007337 = 0.9 sigma
2798 The MISER algorithm reduces the error by a factor of two, and also
2799 correctly estimates the error,
2801 miser ==================
2805 error = -0.002548 = 0.7 sigma
2807 In the case of the VEGAS algorithm the program uses an initial warm-up
2808 run of 10,000 function calls to prepare, or "warm up", the grid. This
2809 is followed by a main run with five iterations of 100,000 function
2810 calls. The chi-squared per degree of freedom for the five iterations are
2811 checked for consistency with 1, and the run is repeated if the results
2812 have not converged. In this case the estimates are consistent on the
2815 vegas warm-up ==================
2819 error = -0.006278 = 2 sigma
2821 result = 1.392957 sigma = 0.000452 chisq/dof = 1.1
2822 vegas final ==================
2826 error = -0.000247 = 0.5 sigma
2828 If the value of `chisq' had differed significantly from 1 it would
2829 indicate inconsistent results, with a correspondingly underestimated
2830 error. The final estimate from VEGAS (using a similar number of
2831 function calls) is significantly more accurate than the other two
2835 File: gsl-ref.info, Node: Monte Carlo Integration References and Further Reading, Prev: Monte Carlo Examples, Up: Monte Carlo Integration
2837 23.6 References and Further Reading
2838 ===================================
2840 The MISER algorithm is described in the following article by Press and
2843 W.H. Press, G.R. Farrar, `Recursive Stratified Sampling for
2844 Multidimensional Monte Carlo Integration', Computers in Physics,
2845 v4 (1990), pp190-195.
2847 The VEGAS algorithm is described in the following papers,
2849 G.P. Lepage, `A New Algorithm for Adaptive Multidimensional
2850 Integration', Journal of Computational Physics 27, 192-203, (1978)
2852 G.P. Lepage, `VEGAS: An Adaptive Multi-dimensional Integration
2853 Program', Cornell preprint CLNS 80-447, March 1980
2856 File: gsl-ref.info, Node: Simulated Annealing, Next: Ordinary Differential Equations, Prev: Monte Carlo Integration, Up: Top
2858 24 Simulated Annealing
2859 **********************
2861 Stochastic search techniques are used when the structure of a space is
2862 not well understood or is not smooth, so that techniques like Newton's
2863 method (which requires calculating Jacobian derivative matrices) cannot
2864 be used. In particular, these techniques are frequently used to solve
2865 combinatorial optimization problems, such as the traveling salesman
2868 The goal is to find a point in the space at which a real valued
2869 "energy function" (or "cost function") is minimized. Simulated
2870 annealing is a minimization technique which has given good results in
2871 avoiding local minima; it is based on the idea of taking a random walk
2872 through the space at successively lower temperatures, where the
2873 probability of taking a step is given by a Boltzmann distribution.
2875 The functions described in this chapter are declared in the header
2880 * Simulated Annealing algorithm::
2881 * Simulated Annealing functions::
2882 * Examples with Simulated Annealing::
2883 * Simulated Annealing References and Further Reading::
2886 File: gsl-ref.info, Node: Simulated Annealing algorithm, Next: Simulated Annealing functions, Up: Simulated Annealing
2888 24.1 Simulated Annealing algorithm
2889 ==================================
2891 The simulated annealing algorithm takes random walks through the problem
2892 space, looking for points with low energies; in these random walks, the
2893 probability of taking a step is determined by the Boltzmann
2896 p = e^{-(E_{i+1} - E_i)/(kT)}
2898 if E_{i+1} > E_i, and p = 1 when E_{i+1} <= E_i.
2900 In other words, a step will occur if the new energy is lower. If
2901 the new energy is higher, the transition can still occur, and its
2902 likelihood is proportional to the temperature T and inversely
2903 proportional to the energy difference E_{i+1} - E_i.
2905 The temperature T is initially set to a high value, and a random
2906 walk is carried out at that temperature. Then the temperature is
2907 lowered very slightly according to a "cooling schedule", for example: T
2908 -> T/mu_T where \mu_T is slightly greater than 1.
2910 The slight probability of taking a step that gives higher energy is
2911 what allows simulated annealing to frequently get out of local minima.
2914 File: gsl-ref.info, Node: Simulated Annealing functions, Next: Examples with Simulated Annealing, Prev: Simulated Annealing algorithm, Up: Simulated Annealing
2916 24.2 Simulated Annealing functions
2917 ==================================
2919 -- Function: void gsl_siman_solve (const gsl_rng * R, void * X0_P,
2920 gsl_siman_Efunc_t EF, gsl_siman_step_t TAKE_STEP,
2921 gsl_siman_metric_t DISTANCE, gsl_siman_print_t
2922 PRINT_POSITION, gsl_siman_copy_t COPYFUNC,
2923 gsl_siman_copy_construct_t COPY_CONSTRUCTOR,
2924 gsl_siman_destroy_t DESTRUCTOR, size_t ELEMENT_SIZE,
2925 gsl_siman_params_t PARAMS)
2926 This function performs a simulated annealing search through a given
2927 space. The space is specified by providing the functions EF and
2928 DISTANCE. The simulated annealing steps are generated using the
2929 random number generator R and the function TAKE_STEP.
2931 The starting configuration of the system should be given by X0_P.
2932 The routine offers two modes for updating configurations, a
2933 fixed-size mode and a variable-size mode. In the fixed-size mode
2934 the configuration is stored as a single block of memory of size
2935 ELEMENT_SIZE. Copies of this configuration are created, copied
2936 and destroyed internally using the standard library functions
2937 `malloc', `memcpy' and `free'. The function pointers COPYFUNC,
2938 COPY_CONSTRUCTOR and DESTRUCTOR should be null pointers in
2939 fixed-size mode. In the variable-size mode the functions
2940 COPYFUNC, COPY_CONSTRUCTOR and DESTRUCTOR are used to create, copy
2941 and destroy configurations internally. The variable ELEMENT_SIZE
2942 should be zero in the variable-size mode.
2944 The PARAMS structure (described below) controls the run by
2945 providing the temperature schedule and other tunable parameters to
2948 On exit the best result achieved during the search is placed in
2949 `*X0_P'. If the annealing process has been successful this should
2950 be a good approximation to the optimal point in the space.
2952 If the function pointer PRINT_POSITION is not null, a debugging
2953 log will be printed to `stdout' with the following columns:
2955 #-iter #-evals temperature position energy best_energy
2957 and the output of the function PRINT_POSITION itself. If
2958 PRINT_POSITION is null then no information is printed.
2960 The simulated annealing routines require several user-specified
2961 functions to define the configuration space and energy function. The
2962 prototypes for these functions are given below.
2964 -- Data Type: gsl_siman_Efunc_t
2965 This function type should return the energy of a configuration XP.
2967 double (*gsl_siman_Efunc_t) (void *xp)
2969 -- Data Type: gsl_siman_step_t
2970 This function type should modify the configuration XP using a
2971 random step taken from the generator R, up to a maximum distance of
2974 void (*gsl_siman_step_t) (const gsl_rng *r, void *xp,
2977 -- Data Type: gsl_siman_metric_t
2978 This function type should return the distance between two
2979 configurations XP and YP.
2981 double (*gsl_siman_metric_t) (void *xp, void *yp)
2983 -- Data Type: gsl_siman_print_t
2984 This function type should print the contents of the configuration
2987 void (*gsl_siman_print_t) (void *xp)
2989 -- Data Type: gsl_siman_copy_t
2990 This function type should copy the configuration SOURCE into DEST.
2992 void (*gsl_siman_copy_t) (void *source, void *dest)
2994 -- Data Type: gsl_siman_copy_construct_t
2995 This function type should create a new copy of the configuration
2998 void * (*gsl_siman_copy_construct_t) (void *xp)
3000 -- Data Type: gsl_siman_destroy_t
3001 This function type should destroy the configuration XP, freeing its
3004 void (*gsl_siman_destroy_t) (void *xp)
3006 -- Data Type: gsl_siman_params_t
3007 These are the parameters that control a run of `gsl_siman_solve'.
3008 This structure contains all the information needed to control the
3009 search, beyond the energy function, the step function and the
3013 The number of points to try for each step.
3016 The number of iterations at each temperature.
3019 The maximum step size in the random walk.
3021 `double k, t_initial, mu_t, t_min'
3022 The parameters of the Boltzmann distribution and cooling
3026 File: gsl-ref.info, Node: Examples with Simulated Annealing, Next: Simulated Annealing References and Further Reading, Prev: Simulated Annealing functions, Up: Simulated Annealing
3031 The simulated annealing package is clumsy, and it has to be because it
3032 is written in C, for C callers, and tries to be polymorphic at the same
3033 time. But here we provide some examples which can be pasted into your
3034 application with little change and should make things easier.
3039 * Traveling Salesman Problem::
3042 File: gsl-ref.info, Node: Trivial example, Next: Traveling Salesman Problem, Up: Examples with Simulated Annealing
3044 24.3.1 Trivial example
3045 ----------------------
3047 The first example, in one dimensional cartesian space, sets up an energy
3048 function which is a damped sine wave; this has many local minima, but
3049 only one global minimum, somewhere between 1.0 and 1.5. The initial
3050 guess given is 15.5, which is several local minima away from the global
3056 #include <gsl/gsl_siman.h>
3058 /* set up parameters for this simulated annealing run */
3060 /* how many points do we try before stepping */
3063 /* how many iterations for each T? */
3064 #define ITERS_FIXED_T 1000
3066 /* max step size in random walk */
3067 #define STEP_SIZE 1.0
3069 /* Boltzmann constant */
3072 /* initial temperature */
3073 #define T_INITIAL 0.008
3075 /* damping factor for temperature */
3077 #define T_MIN 2.0e-6
3079 gsl_siman_params_t params
3080 = {N_TRIES, ITERS_FIXED_T, STEP_SIZE,
3081 K, T_INITIAL, MU_T, T_MIN};
3083 /* now some functions to test in one dimension */
3086 double x = * ((double *) xp);
3088 return exp(-pow((x-1.0),2.0))*sin(8*x);
3091 double M1(void *xp, void *yp)
3093 double x = *((double *) xp);
3094 double y = *((double *) yp);
3099 void S1(const gsl_rng * r, void *xp, double step_size)
3101 double old_x = *((double *) xp);
3104 double u = gsl_rng_uniform(r);
3105 new_x = u * 2 * step_size - step_size + old_x;
3107 memcpy(xp, &new_x, sizeof(new_x));
3112 printf ("%12g", *((double *) xp));
3116 main(int argc, char *argv[])
3118 const gsl_rng_type * T;
3121 double x_initial = 15.5;
3123 gsl_rng_env_setup();
3125 T = gsl_rng_default;
3126 r = gsl_rng_alloc(T);
3128 gsl_siman_solve(r, &x_initial, E1, S1, M1, P1,
3130 sizeof(double), params);
3136 Here are a couple of plots that are generated by running
3137 `siman_test' in the following way:
3139 $ ./siman_test | awk '!/^#/ {print $1, $4}'
3140 | graph -y 1.34 1.4 -W0 -X generation -Y position
3141 | plot -Tps > siman-test.eps
3142 $ ./siman_test | awk '!/^#/ {print $1, $5}'
3143 | graph -y -0.88 -0.83 -W0 -X generation -Y energy
3144 | plot -Tps > siman-energy.eps
3147 File: gsl-ref.info, Node: Traveling Salesman Problem, Prev: Trivial example, Up: Examples with Simulated Annealing
3149 24.3.2 Traveling Salesman Problem
3150 ---------------------------------
3152 The TSP ("Traveling Salesman Problem") is the classic combinatorial
3153 optimization problem. I have provided a very simple version of it,
3154 based on the coordinates of twelve cities in the southwestern United
3155 States. This should maybe be called the "Flying Salesman Problem",
3156 since I am using the great-circle distance between cities, rather than
3157 the driving distance. Also: I assume the earth is a sphere, so I don't
3158 use geoid distances.
3160 The `gsl_siman_solve' routine finds a route which is 3490.62
3161 Kilometers long; this is confirmed by an exhaustive search of all
3162 possible routes with the same initial city.
3164 The full code can be found in `siman/siman_tsp.c', but I include
3165 here some plots generated in the following way:
3167 $ ./siman_tsp > tsp.output
3168 $ grep -v "^#" tsp.output
3169 | awk '{print $1, $NF}'
3170 | graph -y 3300 6500 -W0 -X generation -Y distance
3171 -L "TSP - 12 southwest cities"
3172 | plot -Tps > 12-cities.eps
3173 $ grep initial_city_coord tsp.output
3174 | awk '{print $2, $3}'
3175 | graph -X "longitude (- means west)" -Y "latitude"
3176 -L "TSP - initial-order" -f 0.03 -S 1 0.1
3177 | plot -Tps > initial-route.eps
3178 $ grep final_city_coord tsp.output
3179 | awk '{print $2, $3}'
3180 | graph -X "longitude (- means west)" -Y "latitude"
3181 -L "TSP - final-order" -f 0.03 -S 1 0.1
3182 | plot -Tps > final-route.eps
3184 This is the output showing the initial order of the cities; longitude is
3185 negative, since it is west and I want the plot to look like a map.
3187 # initial coordinates of cities (longitude and latitude)
3188 ###initial_city_coord: -105.95 35.68 Santa Fe
3189 ###initial_city_coord: -112.07 33.54 Phoenix
3190 ###initial_city_coord: -106.62 35.12 Albuquerque
3191 ###initial_city_coord: -103.2 34.41 Clovis
3192 ###initial_city_coord: -107.87 37.29 Durango
3193 ###initial_city_coord: -96.77 32.79 Dallas
3194 ###initial_city_coord: -105.92 35.77 Tesuque
3195 ###initial_city_coord: -107.84 35.15 Grants
3196 ###initial_city_coord: -106.28 35.89 Los Alamos
3197 ###initial_city_coord: -106.76 32.34 Las Cruces
3198 ###initial_city_coord: -108.58 37.35 Cortez
3199 ###initial_city_coord: -108.74 35.52 Gallup
3200 ###initial_city_coord: -105.95 35.68 Santa Fe
3202 The optimal route turns out to be:
3204 # final coordinates of cities (longitude and latitude)
3205 ###final_city_coord: -105.95 35.68 Santa Fe
3206 ###final_city_coord: -103.2 34.41 Clovis
3207 ###final_city_coord: -96.77 32.79 Dallas
3208 ###final_city_coord: -106.76 32.34 Las Cruces
3209 ###final_city_coord: -112.07 33.54 Phoenix
3210 ###final_city_coord: -108.74 35.52 Gallup
3211 ###final_city_coord: -108.58 37.35 Cortez
3212 ###final_city_coord: -107.87 37.29 Durango
3213 ###final_city_coord: -107.84 35.15 Grants
3214 ###final_city_coord: -106.62 35.12 Albuquerque
3215 ###final_city_coord: -106.28 35.89 Los Alamos
3216 ###final_city_coord: -105.92 35.77 Tesuque
3217 ###final_city_coord: -105.95 35.68 Santa Fe
3219 Here's a plot of the cost function (energy) versus generation (point in
3220 the calculation at which a new temperature is set) for this problem:
3223 File: gsl-ref.info, Node: Simulated Annealing References and Further Reading, Prev: Examples with Simulated Annealing, Up: Simulated Annealing
3225 24.4 References and Further Reading
3226 ===================================
3228 Further information is available in the following book,
3230 `Modern Heuristic Techniques for Combinatorial Problems', Colin R.
3231 Reeves (ed.), McGraw-Hill, 1995 (ISBN 0-07-709239-2).
3234 File: gsl-ref.info, Node: Ordinary Differential Equations, Next: Interpolation, Prev: Simulated Annealing, Up: Top
3236 25 Ordinary Differential Equations
3237 **********************************
3239 This chapter describes functions for solving ordinary differential
3240 equation (ODE) initial value problems. The library provides a variety
3241 of low-level methods, such as Runge-Kutta and Bulirsch-Stoer routines,
3242 and higher-level components for adaptive step-size control. The
3243 components can be combined by the user to achieve the desired solution,
3244 with full access to any intermediate steps.
3246 These functions are declared in the header file `gsl_odeiv.h'.
3250 * Defining the ODE System::
3251 * Stepping Functions::
3252 * Adaptive Step-size Control::
3254 * ODE Example programs::
3255 * ODE References and Further Reading::
3258 File: gsl-ref.info, Node: Defining the ODE System, Next: Stepping Functions, Up: Ordinary Differential Equations
3260 25.1 Defining the ODE System
3261 ============================
3263 The routines solve the general n-dimensional first-order system,
3265 dy_i(t)/dt = f_i(t, y_1(t), ..., y_n(t))
3267 for i = 1, \dots, n. The stepping functions rely on the vector of
3268 derivatives f_i and the Jacobian matrix, J_{ij} = df_i(t,y(t)) / dy_j.
3269 A system of equations is defined using the `gsl_odeiv_system' datatype.
3271 -- Data Type: gsl_odeiv_system
3272 This data type defines a general ODE system with arbitrary
3275 `int (* function) (double t, const double y[], double dydt[], void * params)'
3276 This function should store the vector elements
3277 f_i(t,y,params) in the array DYDT, for arguments (T,Y) and
3278 parameters PARAMS. The function should return `GSL_SUCCESS'
3279 if the calculation was completed successfully. Any other
3280 return value indicates an error.
3282 `int (* jacobian) (double t, const double y[], double * dfdy, double dfdt[], void * params);'
3283 This function should store the vector of derivative elements
3284 df_i(t,y,params)/dt in the array DFDT and the Jacobian matrix
3285 J_{ij} in the array DFDY, regarded as a row-ordered matrix
3286 `J(i,j) = dfdy[i * dimension + j]' where `dimension' is the
3287 dimension of the system. The function should return
3288 `GSL_SUCCESS' if the calculation was completed successfully.
3289 Any other return value indicates an error.
3291 Some of the simpler solver algorithms do not make use of the
3292 Jacobian matrix, so it is not always strictly necessary to
3293 provide it (the `jacobian' element of the struct can be
3294 replaced by a null pointer for those algorithms). However,
3295 it is useful to provide the Jacobian to allow the solver
3296 algorithms to be interchanged--the best algorithms make use
3300 This is the dimension of the system of equations.
3303 This is a pointer to the arbitrary parameters of the system.
3306 File: gsl-ref.info, Node: Stepping Functions, Next: Adaptive Step-size Control, Prev: Defining the ODE System, Up: Ordinary Differential Equations
3308 25.2 Stepping Functions
3309 =======================
3311 The lowest level components are the "stepping functions" which advance
3312 a solution from time t to t+h for a fixed step-size h and estimate the
3313 resulting local error.
3315 -- Function: gsl_odeiv_step * gsl_odeiv_step_alloc (const
3316 gsl_odeiv_step_type * T, size_t DIM)
3317 This function returns a pointer to a newly allocated instance of a
3318 stepping function of type T for a system of DIM dimensions.
3320 -- Function: int gsl_odeiv_step_reset (gsl_odeiv_step * S)
3321 This function resets the stepping function S. It should be used
3322 whenever the next use of S will not be a continuation of a
3325 -- Function: void gsl_odeiv_step_free (gsl_odeiv_step * S)
3326 This function frees all the memory associated with the stepping
3329 -- Function: const char * gsl_odeiv_step_name (const gsl_odeiv_step *
3331 This function returns a pointer to the name of the stepping
3332 function. For example,
3334 printf ("step method is '%s'\n",
3335 gsl_odeiv_step_name (s));
3337 would print something like `step method is 'rkf45''.
3339 -- Function: unsigned int gsl_odeiv_step_order (const gsl_odeiv_step *
3341 This function returns the order of the stepping function on the
3342 previous step. This order can vary if the stepping function
3345 -- Function: int gsl_odeiv_step_apply (gsl_odeiv_step * S, double T,
3346 double H, double Y[], double YERR[], const double DYDT_IN[],
3347 double DYDT_OUT[], const gsl_odeiv_system * DYDT)
3348 This function applies the stepping function S to the system of
3349 equations defined by DYDT, using the step size H to advance the
3350 system from time T and state Y to time T+H. The new state of the
3351 system is stored in Y on output, with an estimate of the absolute
3352 error in each component stored in YERR. If the argument DYDT_IN
3353 is not null it should point an array containing the derivatives
3354 for the system at time T on input. This is optional as the
3355 derivatives will be computed internally if they are not provided,
3356 but allows the reuse of existing derivative information. On
3357 output the new derivatives of the system at time T+H will be
3358 stored in DYDT_OUT if it is not null.
3360 If the user-supplied functions defined in the system DYDT return a
3361 status other than `GSL_SUCCESS' the step will be aborted. In this
3362 case, the elements of Y will be restored to their pre-step values
3363 and the error code from the user-supplied function will be
3364 returned. The step-size H will be set to the step-size which
3365 caused the error. If the function is called again with a smaller
3366 step-size, e.g. H/10, it should be possible to get closer to any
3367 singularity. To distinguish between error codes from the
3368 user-supplied functions and those from `gsl_odeiv_step_apply'
3369 itself, any user-defined return values should be distinct from the
3370 standard GSL error codes.
3372 The following algorithms are available,
3374 -- Step Type: gsl_odeiv_step_rk2
3375 Embedded Runge-Kutta (2, 3) method.
3377 -- Step Type: gsl_odeiv_step_rk4
3378 4th order (classical) Runge-Kutta. The error estimate is obtained
3379 by halving the step-size. For more efficient estimate of the
3380 error, use the Runge-Kutta-Fehlberg method described below.
3382 -- Step Type: gsl_odeiv_step_rkf45
3383 Embedded Runge-Kutta-Fehlberg (4, 5) method. This method is a good
3384 general-purpose integrator.
3386 -- Step Type: gsl_odeiv_step_rkck
3387 Embedded Runge-Kutta Cash-Karp (4, 5) method.
3389 -- Step Type: gsl_odeiv_step_rk8pd
3390 Embedded Runge-Kutta Prince-Dormand (8,9) method.
3392 -- Step Type: gsl_odeiv_step_rk2imp
3393 Implicit 2nd order Runge-Kutta at Gaussian points.
3395 -- Step Type: gsl_odeiv_step_rk4imp
3396 Implicit 4th order Runge-Kutta at Gaussian points.
3398 -- Step Type: gsl_odeiv_step_bsimp
3399 Implicit Bulirsch-Stoer method of Bader and Deuflhard. This
3400 algorithm requires the Jacobian.
3402 -- Step Type: gsl_odeiv_step_gear1
3403 M=1 implicit Gear method.
3405 -- Step Type: gsl_odeiv_step_gear2
3406 M=2 implicit Gear method.
3409 File: gsl-ref.info, Node: Adaptive Step-size Control, Next: Evolution, Prev: Stepping Functions, Up: Ordinary Differential Equations
3411 25.3 Adaptive Step-size Control
3412 ===============================
3414 The control function examines the proposed change to the solution
3415 produced by a stepping function and attempts to determine the optimal
3416 step-size for a user-specified level of error.
3418 -- Function: gsl_odeiv_control * gsl_odeiv_control_standard_new
3419 (double EPS_ABS, double EPS_REL, double A_Y, double A_DYDT)
3420 The standard control object is a four parameter heuristic based on
3421 absolute and relative errors EPS_ABS and EPS_REL, and scaling
3422 factors A_Y and A_DYDT for the system state y(t) and derivatives
3425 The step-size adjustment procedure for this method begins by
3426 computing the desired error level D_i for each component,
3428 D_i = eps_abs + eps_rel * (a_y |y_i| + a_dydt h |y'_i|)
3430 and comparing it with the observed error E_i = |yerr_i|. If the
3431 observed error E exceeds the desired error level D by more than
3432 10% for any component then the method reduces the step-size by an
3435 h_new = h_old * S * (E/D)^(-1/q)
3437 where q is the consistency order of the method (e.g. q=4 for 4(5)
3438 embedded RK), and S is a safety factor of 0.9. The ratio E/D is
3439 taken to be the maximum of the ratios E_i/D_i.
3441 If the observed error E is less than 50% of the desired error
3442 level D for the maximum ratio E_i/D_i then the algorithm takes the
3443 opportunity to increase the step-size to bring the error in line
3444 with the desired level,
3446 h_new = h_old * S * (E/D)^(-1/(q+1))
3448 This encompasses all the standard error scaling methods. To avoid
3449 uncontrolled changes in the stepsize, the overall scaling factor is
3450 limited to the range 1/5 to 5.
3452 -- Function: gsl_odeiv_control * gsl_odeiv_control_y_new (double
3453 EPS_ABS, double EPS_REL)
3454 This function creates a new control object which will keep the
3455 local error on each step within an absolute error of EPS_ABS and
3456 relative error of EPS_REL with respect to the solution y_i(t).
3457 This is equivalent to the standard control object with A_Y=1 and
3460 -- Function: gsl_odeiv_control * gsl_odeiv_control_yp_new (double
3461 EPS_ABS, double EPS_REL)
3462 This function creates a new control object which will keep the
3463 local error on each step within an absolute error of EPS_ABS and
3464 relative error of EPS_REL with respect to the derivatives of the
3465 solution y'_i(t). This is equivalent to the standard control
3466 object with A_Y=0 and A_DYDT=1.
3468 -- Function: gsl_odeiv_control * gsl_odeiv_control_scaled_new (double
3469 EPS_ABS, double EPS_REL, double A_Y, double A_DYDT, const
3470 double SCALE_ABS[], size_t DIM)
3471 This function creates a new control object which uses the same
3472 algorithm as `gsl_odeiv_control_standard_new' but with an absolute
3473 error which is scaled for each component by the array SCALE_ABS.
3474 The formula for D_i for this control object is,
3476 D_i = eps_abs * s_i + eps_rel * (a_y |y_i| + a_dydt h |y'_i|)
3478 where s_i is the i-th component of the array SCALE_ABS. The same
3479 error control heuristic is used by the Matlab ODE suite.
3481 -- Function: gsl_odeiv_control * gsl_odeiv_control_alloc (const
3482 gsl_odeiv_control_type * T)
3483 This function returns a pointer to a newly allocated instance of a
3484 control function of type T. This function is only needed for
3485 defining new types of control functions. For most purposes the
3486 standard control functions described above should be sufficient.
3488 -- Function: int gsl_odeiv_control_init (gsl_odeiv_control * C, double
3489 EPS_ABS, double EPS_REL, double A_Y, double A_DYDT)
3490 This function initializes the control function C with the
3491 parameters EPS_ABS (absolute error), EPS_REL (relative error), A_Y
3492 (scaling factor for y) and A_DYDT (scaling factor for derivatives).
3494 -- Function: void gsl_odeiv_control_free (gsl_odeiv_control * C)
3495 This function frees all the memory associated with the control
3498 -- Function: int gsl_odeiv_control_hadjust (gsl_odeiv_control * C,
3499 gsl_odeiv_step * S, const double Y[], const double YERR[],
3500 const double DYDT[], double * H)
3501 This function adjusts the step-size H using the control function
3502 C, and the current values of Y, YERR and DYDT. The stepping
3503 function STEP is also needed to determine the order of the method.
3504 If the error in the y-values YERR is found to be too large then
3505 the step-size H is reduced and the function returns
3506 `GSL_ODEIV_HADJ_DEC'. If the error is sufficiently small then H
3507 may be increased and `GSL_ODEIV_HADJ_INC' is returned. The
3508 function returns `GSL_ODEIV_HADJ_NIL' if the step-size is
3509 unchanged. The goal of the function is to estimate the largest
3510 step-size which satisfies the user-specified accuracy requirements
3511 for the current point.
3513 -- Function: const char * gsl_odeiv_control_name (const
3514 gsl_odeiv_control * C)
3515 This function returns a pointer to the name of the control
3516 function. For example,
3518 printf ("control method is '%s'\n",
3519 gsl_odeiv_control_name (c));
3521 would print something like `control method is 'standard''
3524 File: gsl-ref.info, Node: Evolution, Next: ODE Example programs, Prev: Adaptive Step-size Control, Up: Ordinary Differential Equations
3529 The highest level of the system is the evolution function which combines
3530 the results of a stepping function and control function to reliably
3531 advance the solution forward over an interval (t_0, t_1). If the
3532 control function signals that the step-size should be decreased the
3533 evolution function backs out of the current step and tries the proposed
3534 smaller step-size. This process is continued until an acceptable
3537 -- Function: gsl_odeiv_evolve * gsl_odeiv_evolve_alloc (size_t DIM)
3538 This function returns a pointer to a newly allocated instance of an
3539 evolution function for a system of DIM dimensions.
3541 -- Function: int gsl_odeiv_evolve_apply (gsl_odeiv_evolve * E,
3542 gsl_odeiv_control * CON, gsl_odeiv_step * STEP, const
3543 gsl_odeiv_system * DYDT, double * T, double T1, double * H,
3545 This function advances the system (E, DYDT) from time T and
3546 position Y using the stepping function STEP. The new time and
3547 position are stored in T and Y on output. The initial step-size
3548 is taken as H, but this will be modified using the control
3549 function C to achieve the appropriate error bound if necessary.
3550 The routine may make several calls to STEP in order to determine
3551 the optimum step-size. An estimate of the local error for the step
3552 can be obtained from the components of the array `E->yerr[]'. If
3553 the step-size has been changed the value of H will be modified on
3554 output. The maximum time T1 is guaranteed not to be exceeded by
3555 the time-step. On the final time-step the value of T will be set
3558 If the user-supplied functions defined in the system DYDT return a
3559 status other than `GSL_SUCCESS' the step will be aborted. In this
3560 case, T and Y will be restored to their pre-step values and the
3561 error code from the user-supplied function will be returned. To
3562 distinguish between error codes from the user-supplied functions
3563 and those from `gsl_odeiv_evolve_apply' itself, any user-defined
3564 return values should be distinct from the standard GSL error codes.
3566 -- Function: int gsl_odeiv_evolve_reset (gsl_odeiv_evolve * E)
3567 This function resets the evolution function E. It should be used
3568 whenever the next use of E will not be a continuation of a
3571 -- Function: void gsl_odeiv_evolve_free (gsl_odeiv_evolve * E)
3572 This function frees all the memory associated with the evolution
3575 Where a system has discontinuous changes in the derivatives at known
3576 times it is advisable to evolve the system between each discontinuity
3577 in sequence. For example, if a step-change in an external driving
3578 force occurs at times t_a, t_b, t_c, \dots then evolving over the
3579 ranges (t_0,t_a), (t_a,t_b), ..., (t_c,t_1) is more efficient than
3580 using the single range (t_0,t_1).
3582 Evolving the system directly through a discontinuity with a strict
3583 tolerance may result in extremely small steps being taken at the edge
3584 of the discontinuity (e.g. down to the limit of machine precision). In
3585 this case it may be necessary to impose a minimum step size `hmin'
3586 suitable for the problem:
3590 gsl_odeiv_evolve_apply (e, c, s, &sys, &t, t1, &h, y);
3591 if (h < hmin) { h = hmin; } ;
3593 The value of H returned by `gsl_odeiv_evolve_apply' is always a
3594 suggested value and can be modified whenever needed.
3597 File: gsl-ref.info, Node: ODE Example programs, Next: ODE References and Further Reading, Prev: Evolution, Up: Ordinary Differential Equations
3602 The following program solves the second-order nonlinear Van der Pol
3603 oscillator equation,
3605 x''(t) + \mu x'(t) (x(t)^2 - 1) + x(t) = 0
3607 This can be converted into a first order system suitable for use with
3608 the routines described in this chapter by introducing a separate
3609 variable for the velocity, y = x'(t),
3612 y' = -x + \mu y (1-x^2)
3614 The program begins by defining functions for these derivatives and
3618 #include <gsl/gsl_errno.h>
3619 #include <gsl/gsl_matrix.h>
3620 #include <gsl/gsl_odeiv.h>
3623 func (double t, const double y[], double f[],
3626 double mu = *(double *)params;
3628 f[1] = -y[0] - mu*y[1]*(y[0]*y[0] - 1);
3633 jac (double t, const double y[], double *dfdy,
3634 double dfdt[], void *params)
3636 double mu = *(double *)params;
3637 gsl_matrix_view dfdy_mat
3638 = gsl_matrix_view_array (dfdy, 2, 2);
3639 gsl_matrix * m = &dfdy_mat.matrix;
3640 gsl_matrix_set (m, 0, 0, 0.0);
3641 gsl_matrix_set (m, 0, 1, 1.0);
3642 gsl_matrix_set (m, 1, 0, -2.0*mu*y[0]*y[1] - 1.0);
3643 gsl_matrix_set (m, 1, 1, -mu*(y[0]*y[0] - 1.0));
3652 const gsl_odeiv_step_type * T
3653 = gsl_odeiv_step_rk8pd;
3656 = gsl_odeiv_step_alloc (T, 2);
3657 gsl_odeiv_control * c
3658 = gsl_odeiv_control_y_new (1e-6, 0.0);
3659 gsl_odeiv_evolve * e
3660 = gsl_odeiv_evolve_alloc (2);
3663 gsl_odeiv_system sys = {func, jac, 2, &mu};
3665 double t = 0.0, t1 = 100.0;
3667 double y[2] = { 1.0, 0.0 };
3671 int status = gsl_odeiv_evolve_apply (e, c, s,
3676 if (status != GSL_SUCCESS)
3679 printf ("%.5e %.5e %.5e\n", t, y[0], y[1]);
3682 gsl_odeiv_evolve_free (e);
3683 gsl_odeiv_control_free (c);
3684 gsl_odeiv_step_free (s);
3688 For functions with multiple parameters, the appropriate information can
3689 be passed in through the PARAMS argument using a pointer to a struct.
3691 The main loop of the program evolves the solution from (y, y') = (1,
3692 0) at t=0 to t=100. The step-size h is automatically adjusted by the
3693 controller to maintain an absolute accuracy of 10^{-6} in the function
3696 To obtain the values at user-specified positions, rather than those
3697 chosen automatically by the control function, the main loop can be
3698 modified to advance the solution from one chosen point to the next.
3699 For example, the following main loop prints the solution at the points
3700 t_i = 0, 1, 2, \dots, 100,
3702 for (i = 1; i <= 100; i++)
3704 double ti = i * t1 / 100.0;
3708 gsl_odeiv_evolve_apply (e, c, s,
3714 printf ("%.5e %.5e %.5e\n", t, y[0], y[1]);
3717 Note that arbitrary values of t_i can be used for each stage of the
3718 integration. The equally spaced points in this example are just used
3721 It is also possible to work with a non-adaptive integrator, using
3722 only the stepping function itself. The following program uses the `rk4'
3723 fourth-order Runge-Kutta stepping function with a fixed stepsize of
3729 const gsl_odeiv_step_type * T
3730 = gsl_odeiv_step_rk4;
3733 = gsl_odeiv_step_alloc (T, 2);
3736 gsl_odeiv_system sys = {func, jac, 2, &mu};
3738 double t = 0.0, t1 = 100.0;
3740 double y[2] = { 1.0, 0.0 }, y_err[2];
3741 double dydt_in[2], dydt_out[2];
3743 /* initialise dydt_in from system parameters */
3744 GSL_ODEIV_FN_EVAL(&sys, t, y, dydt_in);
3748 int status = gsl_odeiv_step_apply (s, t, h,
3754 if (status != GSL_SUCCESS)
3757 dydt_in[0] = dydt_out[0];
3758 dydt_in[1] = dydt_out[1];
3762 printf ("%.5e %.5e %.5e\n", t, y[0], y[1]);
3765 gsl_odeiv_step_free (s);
3769 The derivatives must be initialized for the starting point t=0 before
3770 the first step is taken. Subsequent steps use the output derivatives
3771 DYDT_OUT as inputs to the next step by copying their values into
3775 File: gsl-ref.info, Node: ODE References and Further Reading, Prev: ODE Example programs, Up: Ordinary Differential Equations
3777 25.6 References and Further Reading
3778 ===================================
3780 Many of the basic Runge-Kutta formulas can be found in the Handbook of
3781 Mathematical Functions,
3783 Abramowitz & Stegun (eds.), `Handbook of Mathematical Functions',
3786 The implicit Bulirsch-Stoer algorithm `bsimp' is described in the
3789 G. Bader and P. Deuflhard, "A Semi-Implicit Mid-Point Rule for
3790 Stiff Systems of Ordinary Differential Equations.", Numer. Math.
3794 File: gsl-ref.info, Node: Interpolation, Next: Numerical Differentiation, Prev: Ordinary Differential Equations, Up: Top
3799 This chapter describes functions for performing interpolation. The
3800 library provides a variety of interpolation methods, including Cubic
3801 splines and Akima splines. The interpolation types are interchangeable,
3802 allowing different methods to be used without recompiling.
3803 Interpolations can be defined for both normal and periodic boundary
3804 conditions. Additional functions are available for computing
3805 derivatives and integrals of interpolating functions.
3807 The functions described in this section are declared in the header
3808 files `gsl_interp.h' and `gsl_spline.h'.
3812 * Introduction to Interpolation::
3813 * Interpolation Functions::
3814 * Interpolation Types::
3815 * Index Look-up and Acceleration::
3816 * Evaluation of Interpolating Functions::
3817 * Higher-level Interface::
3818 * Interpolation Example programs::
3819 * Interpolation References and Further Reading::
3822 File: gsl-ref.info, Node: Introduction to Interpolation, Next: Interpolation Functions, Up: Interpolation
3827 Given a set of data points (x_1, y_1) \dots (x_n, y_n) the routines
3828 described in this section compute a continuous interpolating function
3829 y(x) such that y(x_i) = y_i. The interpolation is piecewise smooth,
3830 and its behavior at the end-points is determined by the type of
3834 File: gsl-ref.info, Node: Interpolation Functions, Next: Interpolation Types, Prev: Introduction to Interpolation, Up: Interpolation
3836 26.2 Interpolation Functions
3837 ============================
3839 The interpolation function for a given dataset is stored in a
3840 `gsl_interp' object. These are created by the following functions.
3842 -- Function: gsl_interp * gsl_interp_alloc (const gsl_interp_type * T,
3844 This function returns a pointer to a newly allocated interpolation
3845 object of type T for SIZE data-points.
3847 -- Function: int gsl_interp_init (gsl_interp * INTERP, const double
3848 XA[], const double YA[], size_t SIZE)
3849 This function initializes the interpolation object INTERP for the
3850 data (XA,YA) where XA and YA are arrays of size SIZE. The
3851 interpolation object (`gsl_interp') does not save the data arrays
3852 XA and YA and only stores the static state computed from the data.
3853 The XA data array is always assumed to be strictly ordered, with
3854 increasing x values; the behavior for other arrangements is not
3857 -- Function: void gsl_interp_free (gsl_interp * INTERP)
3858 This function frees the interpolation object INTERP.
3861 File: gsl-ref.info, Node: Interpolation Types, Next: Index Look-up and Acceleration, Prev: Interpolation Functions, Up: Interpolation
3863 26.3 Interpolation Types
3864 ========================
3866 The interpolation library provides five interpolation types:
3868 -- Interpolation Type: gsl_interp_linear
3869 Linear interpolation. This interpolation method does not require
3870 any additional memory.
3872 -- Interpolation Type: gsl_interp_polynomial
3873 Polynomial interpolation. This method should only be used for
3874 interpolating small numbers of points because polynomial
3875 interpolation introduces large oscillations, even for well-behaved
3876 datasets. The number of terms in the interpolating polynomial is
3877 equal to the number of points.
3879 -- Interpolation Type: gsl_interp_cspline
3880 Cubic spline with natural boundary conditions. The resulting
3881 curve is piecewise cubic on each interval, with matching first and
3882 second derivatives at the supplied data-points. The second
3883 derivative is chosen to be zero at the first point and last point.
3885 -- Interpolation Type: gsl_interp_cspline_periodic
3886 Cubic spline with periodic boundary conditions. The resulting
3887 curve is piecewise cubic on each interval, with matching first and
3888 second derivatives at the supplied data-points. The derivatives
3889 at the first and last points are also matched. Note that the last
3890 point in the data must have the same y-value as the first point,
3891 otherwise the resulting periodic interpolation will have a
3892 discontinuity at the boundary.
3895 -- Interpolation Type: gsl_interp_akima
3896 Non-rounded Akima spline with natural boundary conditions. This
3897 method uses the non-rounded corner algorithm of Wodicka.
3899 -- Interpolation Type: gsl_interp_akima_periodic
3900 Non-rounded Akima spline with periodic boundary conditions. This
3901 method uses the non-rounded corner algorithm of Wodicka.
3903 The following related functions are available:
3905 -- Function: const char * gsl_interp_name (const gsl_interp * INTERP)
3906 This function returns the name of the interpolation type used by
3907 INTERP. For example,
3909 printf ("interp uses '%s' interpolation.\n",
3910 gsl_interp_name (interp));
3912 would print something like,
3914 interp uses 'cspline' interpolation.
3916 -- Function: unsigned int gsl_interp_min_size (const gsl_interp *
3918 This function returns the minimum number of points required by the
3919 interpolation type of INTERP. For example, Akima spline
3920 interpolation requires a minimum of 5 points.
3923 File: gsl-ref.info, Node: Index Look-up and Acceleration, Next: Evaluation of Interpolating Functions, Prev: Interpolation Types, Up: Interpolation
3925 26.4 Index Look-up and Acceleration
3926 ===================================
3928 The state of searches can be stored in a `gsl_interp_accel' object,
3929 which is a kind of iterator for interpolation lookups. It caches the
3930 previous value of an index lookup. When the subsequent interpolation
3931 point falls in the same interval its index value can be returned
3934 -- Function: size_t gsl_interp_bsearch (const double X_ARRAY[], double
3935 X, size_t INDEX_LO, size_t INDEX_HI)
3936 This function returns the index i of the array X_ARRAY such that
3937 `x_array[i] <= x < x_array[i+1]'. The index is searched for in
3938 the range [INDEX_LO,INDEX_HI]. An inline version of this function
3939 is used when `HAVE_INLINE' is defined.
3941 -- Function: gsl_interp_accel * gsl_interp_accel_alloc (void)
3942 This function returns a pointer to an accelerator object, which is
3943 a kind of iterator for interpolation lookups. It tracks the state
3944 of lookups, thus allowing for application of various acceleration
3947 -- Function: size_t gsl_interp_accel_find (gsl_interp_accel * A, const
3948 double X_ARRAY[], size_t SIZE, double X)
3949 This function performs a lookup action on the data array X_ARRAY
3950 of size SIZE, using the given accelerator A. This is how lookups
3951 are performed during evaluation of an interpolation. The function
3952 returns an index i such that `x_array[i] <= x < x_array[i+1]'. An
3953 inline version of this function is used when `HAVE_INLINE' is
3956 -- Function: void gsl_interp_accel_free (gsl_interp_accel* ACC)
3957 This function frees the accelerator object ACC.
3960 File: gsl-ref.info, Node: Evaluation of Interpolating Functions, Next: Higher-level Interface, Prev: Index Look-up and Acceleration, Up: Interpolation
3962 26.5 Evaluation of Interpolating Functions
3963 ==========================================
3965 -- Function: double gsl_interp_eval (const gsl_interp * INTERP, const
3966 double XA[], const double YA[], double X, gsl_interp_accel *
3968 -- Function: int gsl_interp_eval_e (const gsl_interp * INTERP, const
3969 double XA[], const double YA[], double X, gsl_interp_accel *
3971 These functions return the interpolated value of Y for a given
3972 point X, using the interpolation object INTERP, data arrays XA and
3973 YA and the accelerator ACC.
3975 -- Function: double gsl_interp_eval_deriv (const gsl_interp * INTERP,
3976 const double XA[], const double YA[], double X,
3977 gsl_interp_accel * ACC)
3978 -- Function: int gsl_interp_eval_deriv_e (const gsl_interp * INTERP,
3979 const double XA[], const double YA[], double X,
3980 gsl_interp_accel * ACC, double * D)
3981 These functions return the derivative D of an interpolated
3982 function for a given point X, using the interpolation object
3983 INTERP, data arrays XA and YA and the accelerator ACC.
3985 -- Function: double gsl_interp_eval_deriv2 (const gsl_interp * INTERP,
3986 const double XA[], const double YA[], double X,
3987 gsl_interp_accel * ACC)
3988 -- Function: int gsl_interp_eval_deriv2_e (const gsl_interp * INTERP,
3989 const double XA[], const double YA[], double X,
3990 gsl_interp_accel * ACC, double * D2)
3991 These functions return the second derivative D2 of an interpolated
3992 function for a given point X, using the interpolation object
3993 INTERP, data arrays XA and YA and the accelerator ACC.
3995 -- Function: double gsl_interp_eval_integ (const gsl_interp * INTERP,
3996 const double XA[], const double YA[], double A, double B,
3997 gsl_interp_accel * ACC)
3998 -- Function: int gsl_interp_eval_integ_e (const gsl_interp * INTERP,
3999 const double XA[], const double YA[], double A, double B,
4000 gsl_interp_accel * ACC, double * RESULT)
4001 These functions return the numerical integral RESULT of an
4002 interpolated function over the range [A, B], using the
4003 interpolation object INTERP, data arrays XA and YA and the
4007 File: gsl-ref.info, Node: Higher-level Interface, Next: Interpolation Example programs, Prev: Evaluation of Interpolating Functions, Up: Interpolation
4009 26.6 Higher-level Interface
4010 ===========================
4012 The functions described in the previous sections required the user to
4013 supply pointers to the x and y arrays on each call. The following
4014 functions are equivalent to the corresponding `gsl_interp' functions
4015 but maintain a copy of this data in the `gsl_spline' object. This
4016 removes the need to pass both XA and YA as arguments on each
4017 evaluation. These functions are defined in the header file
4020 -- Function: gsl_spline * gsl_spline_alloc (const gsl_interp_type * T,
4023 -- Function: int gsl_spline_init (gsl_spline * SPLINE, const double
4024 XA[], const double YA[], size_t SIZE)
4026 -- Function: void gsl_spline_free (gsl_spline * SPLINE)
4028 -- Function: const char * gsl_spline_name (const gsl_spline * SPLINE)
4030 -- Function: unsigned int gsl_spline_min_size (const gsl_spline *
4033 -- Function: double gsl_spline_eval (const gsl_spline * SPLINE, double
4034 X, gsl_interp_accel * ACC)
4035 -- Function: int gsl_spline_eval_e (const gsl_spline * SPLINE, double
4036 X, gsl_interp_accel * ACC, double * Y)
4038 -- Function: double gsl_spline_eval_deriv (const gsl_spline * SPLINE,
4039 double X, gsl_interp_accel * ACC)
4040 -- Function: int gsl_spline_eval_deriv_e (const gsl_spline * SPLINE,
4041 double X, gsl_interp_accel * ACC, double * D)
4043 -- Function: double gsl_spline_eval_deriv2 (const gsl_spline * SPLINE,
4044 double X, gsl_interp_accel * ACC)
4045 -- Function: int gsl_spline_eval_deriv2_e (const gsl_spline * SPLINE,
4046 double X, gsl_interp_accel * ACC, double * D2)
4048 -- Function: double gsl_spline_eval_integ (const gsl_spline * SPLINE,
4049 double A, double B, gsl_interp_accel * ACC)
4050 -- Function: int gsl_spline_eval_integ_e (const gsl_spline * SPLINE,
4051 double A, double B, gsl_interp_accel * ACC, double * RESULT)
4054 File: gsl-ref.info, Node: Interpolation Example programs, Next: Interpolation References and Further Reading, Prev: Higher-level Interface, Up: Interpolation
4059 The following program demonstrates the use of the interpolation and
4060 spline functions. It computes a cubic spline interpolation of the
4061 10-point dataset (x_i, y_i) where x_i = i + \sin(i)/2 and y_i = i +
4062 \cos(i^2) for i = 0 \dots 9.
4067 #include <gsl/gsl_errno.h>
4068 #include <gsl/gsl_spline.h>
4074 double xi, yi, x[10], y[10];
4076 printf ("#m=0,S=2\n");
4078 for (i = 0; i < 10; i++)
4080 x[i] = i + 0.5 * sin (i);
4081 y[i] = i + cos (i * i);
4082 printf ("%g %g\n", x[i], y[i]);
4085 printf ("#m=1,S=0\n");
4088 gsl_interp_accel *acc
4089 = gsl_interp_accel_alloc ();
4091 = gsl_spline_alloc (gsl_interp_cspline, 10);
4093 gsl_spline_init (spline, x, y, 10);
4095 for (xi = x[0]; xi < x[9]; xi += 0.01)
4097 yi = gsl_spline_eval (spline, xi, acc);
4098 printf ("%g %g\n", xi, yi);
4100 gsl_spline_free (spline);
4101 gsl_interp_accel_free (acc);
4106 The output is designed to be used with the GNU plotutils `graph'
4109 $ ./a.out > interp.dat
4110 $ graph -T ps < interp.dat > interp.ps
4112 The result shows a smooth interpolation of the original points. The
4113 interpolation method can be changed simply by varying the first
4114 argument of `gsl_spline_alloc'.
4116 The next program demonstrates a periodic cubic spline with 4 data
4117 points. Note that the first and last points must be supplied with the
4118 same y-value for a periodic spline.
4123 #include <gsl/gsl_errno.h>
4124 #include <gsl/gsl_spline.h>
4130 double x[4] = {0.00, 0.10, 0.27, 0.30};
4131 double y[4] = {0.15, 0.70, -0.10, 0.15}; /* Note: first = last
4132 for periodic data */
4134 gsl_interp_accel *acc = gsl_interp_accel_alloc ();
4135 const gsl_interp_type *t = gsl_interp_cspline_periodic;
4136 gsl_spline *spline = gsl_spline_alloc (t, N);
4138 int i; double xi, yi;
4140 printf ("#m=0,S=5\n");
4141 for (i = 0; i < N; i++)
4143 printf ("%g %g\n", x[i], y[i]);
4146 printf ("#m=1,S=0\n");
4147 gsl_spline_init (spline, x, y, N);
4149 for (i = 0; i <= 100; i++)
4151 xi = (1 - i / 100.0) * x[0] + (i / 100.0) * x[N-1];
4152 yi = gsl_spline_eval (spline, xi, acc);
4153 printf ("%g %g\n", xi, yi);
4156 gsl_spline_free (spline);
4157 gsl_interp_accel_free (acc);
4161 The output can be plotted with GNU `graph'.
4163 $ ./a.out > interp.dat
4164 $ graph -T ps < interp.dat > interp.ps
4166 The result shows a periodic interpolation of the original points. The
4167 slope of the fitted curve is the same at the beginning and end of the
4168 data, and the second derivative is also.
4171 File: gsl-ref.info, Node: Interpolation References and Further Reading, Prev: Interpolation Example programs, Up: Interpolation
4173 26.8 References and Further Reading
4174 ===================================
4176 Descriptions of the interpolation algorithms and further references can
4177 be found in the following books:
4179 C.W. Ueberhuber, `Numerical Computation (Volume 1), Chapter 9
4180 "Interpolation"', Springer (1997), ISBN 3-540-62058-3.
4182 D.M. Young, R.T. Gregory `A Survey of Numerical Mathematics
4183 (Volume 1), Chapter 6.8', Dover (1988), ISBN 0-486-65691-8.
4187 File: gsl-ref.info, Node: Numerical Differentiation, Next: Chebyshev Approximations, Prev: Interpolation, Up: Top
4189 27 Numerical Differentiation
4190 ****************************
4192 The functions described in this chapter compute numerical derivatives by
4193 finite differencing. An adaptive algorithm is used to find the best
4194 choice of finite difference and to estimate the error in the derivative.
4195 These functions are declared in the header file `gsl_deriv.h'.
4199 * Numerical Differentiation functions::
4200 * Numerical Differentiation Examples::
4201 * Numerical Differentiation References::
4204 File: gsl-ref.info, Node: Numerical Differentiation functions, Next: Numerical Differentiation Examples, Up: Numerical Differentiation
4209 -- Function: int gsl_deriv_central (const gsl_function * F, double X,
4210 double H, double * RESULT, double * ABSERR)
4211 This function computes the numerical derivative of the function F
4212 at the point X using an adaptive central difference algorithm with
4213 a step-size of H. The derivative is returned in RESULT and an
4214 estimate of its absolute error is returned in ABSERR.
4216 The initial value of H is used to estimate an optimal step-size,
4217 based on the scaling of the truncation error and round-off error
4218 in the derivative calculation. The derivative is computed using a
4219 5-point rule for equally spaced abscissae at x-h, x-h/2, x, x+h/2,
4220 x+h, with an error estimate taken from the difference between the
4221 5-point rule and the corresponding 3-point rule x-h, x, x+h. Note
4222 that the value of the function at x does not contribute to the
4223 derivative calculation, so only 4-points are actually used.
4225 -- Function: int gsl_deriv_forward (const gsl_function * F, double X,
4226 double H, double * RESULT, double * ABSERR)
4227 This function computes the numerical derivative of the function F
4228 at the point X using an adaptive forward difference algorithm with
4229 a step-size of H. The function is evaluated only at points greater
4230 than X, and never at X itself. The derivative is returned in
4231 RESULT and an estimate of its absolute error is returned in
4232 ABSERR. This function should be used if f(x) has a discontinuity
4233 at X, or is undefined for values less than X.
4235 The initial value of H is used to estimate an optimal step-size,
4236 based on the scaling of the truncation error and round-off error
4237 in the derivative calculation. The derivative at x is computed
4238 using an "open" 4-point rule for equally spaced abscissae at x+h/4,
4239 x+h/2, x+3h/4, x+h, with an error estimate taken from the
4240 difference between the 4-point rule and the corresponding 2-point
4243 -- Function: int gsl_deriv_backward (const gsl_function * F, double X,
4244 double H, double * RESULT, double * ABSERR)
4245 This function computes the numerical derivative of the function F
4246 at the point X using an adaptive backward difference algorithm
4247 with a step-size of H. The function is evaluated only at points
4248 less than X, and never at X itself. The derivative is returned in
4249 RESULT and an estimate of its absolute error is returned in
4250 ABSERR. This function should be used if f(x) has a discontinuity
4251 at X, or is undefined for values greater than X.
4253 This function is equivalent to calling `gsl_deriv_forward' with a
4257 File: gsl-ref.info, Node: Numerical Differentiation Examples, Next: Numerical Differentiation References, Prev: Numerical Differentiation functions, Up: Numerical Differentiation
4262 The following code estimates the derivative of the function f(x) =
4263 x^{3/2} at x=2 and at x=0. The function f(x) is undefined for x<0 so
4264 the derivative at x=0 is computed using `gsl_deriv_forward'.
4267 #include <gsl/gsl_math.h>
4268 #include <gsl/gsl_deriv.h>
4270 double f (double x, void * params)
4272 return pow (x, 1.5);
4279 double result, abserr;
4284 printf ("f(x) = x^(3/2)\n");
4286 gsl_deriv_central (&F, 2.0, 1e-8, &result, &abserr);
4287 printf ("x = 2.0\n");
4288 printf ("f'(x) = %.10f +/- %.10f\n", result, abserr);
4289 printf ("exact = %.10f\n\n", 1.5 * sqrt(2.0));
4291 gsl_deriv_forward (&F, 0.0, 1e-8, &result, &abserr);
4292 printf ("x = 0.0\n");
4293 printf ("f'(x) = %.10f +/- %.10f\n", result, abserr);
4294 printf ("exact = %.10f\n", 0.0);
4299 Here is the output of the program,
4304 f'(x) = 2.1213203120 +/- 0.0000004064
4305 exact = 2.1213203436
4308 f'(x) = 0.0000000160 +/- 0.0000000339
4309 exact = 0.0000000000
4312 File: gsl-ref.info, Node: Numerical Differentiation References, Prev: Numerical Differentiation Examples, Up: Numerical Differentiation
4314 27.3 References and Further Reading
4315 ===================================
4317 The algorithms used by these functions are described in the following
4320 Abramowitz and Stegun, `Handbook of Mathematical Functions',
4321 Section 25.3.4, and Table 25.5 (Coefficients for Differentiation).
4323 S.D. Conte and Carl de Boor, `Elementary Numerical Analysis: An
4324 Algorithmic Approach', McGraw-Hill, 1972.
4327 File: gsl-ref.info, Node: Chebyshev Approximations, Next: Series Acceleration, Prev: Numerical Differentiation, Up: Top
4329 28 Chebyshev Approximations
4330 ***************************
4332 This chapter describes routines for computing Chebyshev approximations
4333 to univariate functions. A Chebyshev approximation is a truncation of
4334 the series f(x) = \sum c_n T_n(x), where the Chebyshev polynomials
4335 T_n(x) = \cos(n \arccos x) provide an orthogonal basis of polynomials
4336 on the interval [-1,1] with the weight function 1 / \sqrt{1-x^2}. The
4337 first few Chebyshev polynomials are, T_0(x) = 1, T_1(x) = x, T_2(x) = 2
4338 x^2 - 1. For further information see Abramowitz & Stegun, Chapter 22.
4340 The functions described in this chapter are declared in the header
4341 file `gsl_chebyshev.h'.
4345 * Chebyshev Definitions::
4346 * Creation and Calculation of Chebyshev Series::
4347 * Chebyshev Series Evaluation::
4348 * Derivatives and Integrals::
4349 * Chebyshev Approximation examples::
4350 * Chebyshev Approximation References and Further Reading::
4353 File: gsl-ref.info, Node: Chebyshev Definitions, Next: Creation and Calculation of Chebyshev Series, Up: Chebyshev Approximations
4358 A Chebyshev series is stored using the following structure,
4362 double * c; /* coefficients c[0] .. c[order] */
4363 int order; /* order of expansion */
4364 double a; /* lower interval point */
4365 double b; /* upper interval point */
4369 The approximation is made over the range [a,b] using ORDER+1 terms,
4370 including the coefficient c[0]. The series is computed using the
4371 following convention,
4373 f(x) = (c_0 / 2) + \sum_{n=1} c_n T_n(x)
4375 which is needed when accessing the coefficients directly.
4378 File: gsl-ref.info, Node: Creation and Calculation of Chebyshev Series, Next: Chebyshev Series Evaluation, Prev: Chebyshev Definitions, Up: Chebyshev Approximations
4380 28.2 Creation and Calculation of Chebyshev Series
4381 =================================================
4383 -- Function: gsl_cheb_series * gsl_cheb_alloc (const size_t N)
4384 This function allocates space for a Chebyshev series of order N
4385 and returns a pointer to a new `gsl_cheb_series' struct.
4387 -- Function: void gsl_cheb_free (gsl_cheb_series * CS)
4388 This function frees a previously allocated Chebyshev series CS.
4390 -- Function: int gsl_cheb_init (gsl_cheb_series * CS, const
4391 gsl_function * F, const double A, const double B)
4392 This function computes the Chebyshev approximation CS for the
4393 function F over the range (a,b) to the previously specified order.
4394 The computation of the Chebyshev approximation is an O(n^2)
4395 process, and requires n function evaluations.
4398 File: gsl-ref.info, Node: Chebyshev Series Evaluation, Next: Derivatives and Integrals, Prev: Creation and Calculation of Chebyshev Series, Up: Chebyshev Approximations
4400 28.3 Chebyshev Series Evaluation
4401 ================================
4403 -- Function: double gsl_cheb_eval (const gsl_cheb_series * CS, double
4405 This function evaluates the Chebyshev series CS at a given point X.
4407 -- Function: int gsl_cheb_eval_err (const gsl_cheb_series * CS, const
4408 double X, double * RESULT, double * ABSERR)
4409 This function computes the Chebyshev series CS at a given point X,
4410 estimating both the series RESULT and its absolute error ABSERR.
4411 The error estimate is made from the first neglected term in the
4414 -- Function: double gsl_cheb_eval_n (const gsl_cheb_series * CS,
4415 size_t ORDER, double X)
4416 This function evaluates the Chebyshev series CS at a given point
4417 N, to (at most) the given order ORDER.
4419 -- Function: int gsl_cheb_eval_n_err (const gsl_cheb_series * CS,
4420 const size_t ORDER, const double X, double * RESULT, double *
4422 This function evaluates a Chebyshev series CS at a given point X,
4423 estimating both the series RESULT and its absolute error ABSERR,
4424 to (at most) the given order ORDER. The error estimate is made
4425 from the first neglected term in the series.
4428 File: gsl-ref.info, Node: Derivatives and Integrals, Next: Chebyshev Approximation examples, Prev: Chebyshev Series Evaluation, Up: Chebyshev Approximations
4430 28.4 Derivatives and Integrals
4431 ==============================
4433 The following functions allow a Chebyshev series to be differentiated or
4434 integrated, producing a new Chebyshev series. Note that the error
4435 estimate produced by evaluating the derivative series will be
4436 underestimated due to the contribution of higher order terms being
4439 -- Function: int gsl_cheb_calc_deriv (gsl_cheb_series * DERIV, const
4440 gsl_cheb_series * CS)
4441 This function computes the derivative of the series CS, storing
4442 the derivative coefficients in the previously allocated DERIV.
4443 The two series CS and DERIV must have been allocated with the same
4446 -- Function: int gsl_cheb_calc_integ (gsl_cheb_series * INTEG, const
4447 gsl_cheb_series * CS)
4448 This function computes the integral of the series CS, storing the
4449 integral coefficients in the previously allocated INTEG. The two
4450 series CS and INTEG must have been allocated with the same order.
4451 The lower limit of the integration is taken to be the left hand
4455 File: gsl-ref.info, Node: Chebyshev Approximation examples, Next: Chebyshev Approximation References and Further Reading, Prev: Derivatives and Integrals, Up: Chebyshev Approximations
4460 The following example program computes Chebyshev approximations to a
4461 step function. This is an extremely difficult approximation to make,
4462 due to the discontinuity, and was chosen as an example where
4463 approximation error is visible. For smooth functions the Chebyshev
4464 approximation converges extremely rapidly and errors would not be
4468 #include <gsl/gsl_math.h>
4469 #include <gsl/gsl_chebyshev.h>
4472 f (double x, void *p)
4485 gsl_cheb_series *cs = gsl_cheb_alloc (40);
4492 gsl_cheb_init (cs, &F, 0.0, 1.0);
4494 for (i = 0; i < n; i++)
4496 double x = i / (double)n;
4497 double r10 = gsl_cheb_eval_n (cs, 10, x);
4498 double r40 = gsl_cheb_eval (cs, x);
4499 printf ("%g %g %g %g\n",
4500 x, GSL_FN_EVAL (&F, x), r10, r40);
4508 The output from the program gives the original function, 10-th order
4509 approximation and 40-th order approximation, all sampled at intervals of
4513 File: gsl-ref.info, Node: Chebyshev Approximation References and Further Reading, Prev: Chebyshev Approximation examples, Up: Chebyshev Approximations
4515 28.6 References and Further Reading
4516 ===================================
4518 The following paper describes the use of Chebyshev series,
4520 R. Broucke, "Ten Subroutines for the Manipulation of Chebyshev
4521 Series [C1] (Algorithm 446)". `Communications of the ACM' 16(4),
4525 File: gsl-ref.info, Node: Series Acceleration, Next: Wavelet Transforms, Prev: Chebyshev Approximations, Up: Top
4527 29 Series Acceleration
4528 **********************
4530 The functions described in this chapter accelerate the convergence of a
4531 series using the Levin u-transform. This method takes a small number of
4532 terms from the start of a series and uses a systematic approximation to
4533 compute an extrapolated value and an estimate of its error. The
4534 u-transform works for both convergent and divergent series, including
4537 These functions are declared in the header file `gsl_sum.h'.
4541 * Acceleration functions::
4542 * Acceleration functions without error estimation::
4543 * Example of accelerating a series::
4544 * Series Acceleration References::
4547 File: gsl-ref.info, Node: Acceleration functions, Next: Acceleration functions without error estimation, Up: Series Acceleration
4549 29.1 Acceleration functions
4550 ===========================
4552 The following functions compute the full Levin u-transform of a series
4553 with its error estimate. The error estimate is computed by propagating
4554 rounding errors from each term through to the final extrapolation.
4556 These functions are intended for summing analytic series where each
4557 term is known to high accuracy, and the rounding errors are assumed to
4558 originate from finite precision. They are taken to be relative errors of
4559 order `GSL_DBL_EPSILON' for each term.
4561 The calculation of the error in the extrapolated value is an O(N^2)
4562 process, which is expensive in time and memory. A faster but less
4563 reliable method which estimates the error from the convergence of the
4564 extrapolated value is described in the next section. For the method
4565 described here a full table of intermediate values and derivatives
4566 through to O(N) must be computed and stored, but this does give a
4567 reliable error estimate.
4569 -- Function: gsl_sum_levin_u_workspace * gsl_sum_levin_u_alloc (size_t
4571 This function allocates a workspace for a Levin u-transform of N
4572 terms. The size of the workspace is O(2n^2 + 3n).
4574 -- Function: void gsl_sum_levin_u_free (gsl_sum_levin_u_workspace * W)
4575 This function frees the memory associated with the workspace W.
4577 -- Function: int gsl_sum_levin_u_accel (const double * ARRAY, size_t
4578 ARRAY_SIZE, gsl_sum_levin_u_workspace * W, double *
4579 SUM_ACCEL, double * ABSERR)
4580 This function takes the terms of a series in ARRAY of size
4581 ARRAY_SIZE and computes the extrapolated limit of the series using
4582 a Levin u-transform. Additional working space must be provided in
4583 W. The extrapolated sum is stored in SUM_ACCEL, with an estimate
4584 of the absolute error stored in ABSERR. The actual term-by-term
4585 sum is returned in `w->sum_plain'. The algorithm calculates the
4586 truncation error (the difference between two successive
4587 extrapolations) and round-off error (propagated from the individual
4588 terms) to choose an optimal number of terms for the extrapolation.
4589 All the terms of the series passed in through ARRAY should be
4593 File: gsl-ref.info, Node: Acceleration functions without error estimation, Next: Example of accelerating a series, Prev: Acceleration functions, Up: Series Acceleration
4595 29.2 Acceleration functions without error estimation
4596 ====================================================
4598 The functions described in this section compute the Levin u-transform of
4599 series and attempt to estimate the error from the "truncation error" in
4600 the extrapolation, the difference between the final two approximations.
4601 Using this method avoids the need to compute an intermediate table of
4602 derivatives because the error is estimated from the behavior of the
4603 extrapolated value itself. Consequently this algorithm is an O(N)
4604 process and only requires O(N) terms of storage. If the series
4605 converges sufficiently fast then this procedure can be acceptable. It
4606 is appropriate to use this method when there is a need to compute many
4607 extrapolations of series with similar convergence properties at
4608 high-speed. For example, when numerically integrating a function
4609 defined by a parameterized series where the parameter varies only
4610 slightly. A reliable error estimate should be computed first using the
4611 full algorithm described above in order to verify the consistency of the
4614 -- Function: gsl_sum_levin_utrunc_workspace *
4615 gsl_sum_levin_utrunc_alloc (size_t N)
4616 This function allocates a workspace for a Levin u-transform of N
4617 terms, without error estimation. The size of the workspace is
4620 -- Function: void gsl_sum_levin_utrunc_free
4621 (gsl_sum_levin_utrunc_workspace * W)
4622 This function frees the memory associated with the workspace W.
4624 -- Function: int gsl_sum_levin_utrunc_accel (const double * ARRAY,
4625 size_t ARRAY_SIZE, gsl_sum_levin_utrunc_workspace * W, double
4626 * SUM_ACCEL, double * ABSERR_TRUNC)
4627 This function takes the terms of a series in ARRAY of size
4628 ARRAY_SIZE and computes the extrapolated limit of the series using
4629 a Levin u-transform. Additional working space must be provided in
4630 W. The extrapolated sum is stored in SUM_ACCEL. The actual
4631 term-by-term sum is returned in `w->sum_plain'. The algorithm
4632 terminates when the difference between two successive
4633 extrapolations reaches a minimum or is sufficiently small. The
4634 difference between these two values is used as estimate of the
4635 error and is stored in ABSERR_TRUNC. To improve the reliability
4636 of the algorithm the extrapolated values are replaced by moving
4637 averages when calculating the truncation error, smoothing out any
4641 File: gsl-ref.info, Node: Example of accelerating a series, Next: Series Acceleration References, Prev: Acceleration functions without error estimation, Up: Series Acceleration
4646 The following code calculates an estimate of \zeta(2) = \pi^2 / 6 using
4649 \zeta(2) = 1 + 1/2^2 + 1/3^2 + 1/4^2 + ...
4651 After N terms the error in the sum is O(1/N), making direct summation
4652 of the series converge slowly.
4655 #include <gsl/gsl_math.h>
4656 #include <gsl/gsl_sum.h>
4664 double sum_accel, err;
4668 gsl_sum_levin_u_workspace * w
4669 = gsl_sum_levin_u_alloc (N);
4671 const double zeta_2 = M_PI * M_PI / 6.0;
4673 /* terms for zeta(2) = \sum_{n=1}^{\infty} 1/n^2 */
4675 for (n = 0; n < N; n++)
4677 double np1 = n + 1.0;
4678 t[n] = 1.0 / (np1 * np1);
4682 gsl_sum_levin_u_accel (t, N, w, &sum_accel, &err);
4684 printf ("term-by-term sum = % .16f using %d terms\n",
4687 printf ("term-by-term sum = % .16f using %d terms\n",
4688 w->sum_plain, w->terms_used);
4690 printf ("exact value = % .16f\n", zeta_2);
4691 printf ("accelerated sum = % .16f using %d terms\n",
4692 sum_accel, w->terms_used);
4694 printf ("estimated error = % .16f\n", err);
4695 printf ("actual error = % .16f\n",
4696 sum_accel - zeta_2);
4698 gsl_sum_levin_u_free (w);
4702 The output below shows that the Levin u-transform is able to obtain an
4703 estimate of the sum to 1 part in 10^10 using the first eleven terms of
4704 the series. The error estimate returned by the function is also
4705 accurate, giving the correct number of significant digits.
4708 term-by-term sum = 1.5961632439130233 using 20 terms
4709 term-by-term sum = 1.5759958390005426 using 13 terms
4710 exact value = 1.6449340668482264
4711 accelerated sum = 1.6449340668166479 using 13 terms
4712 estimated error = 0.0000000000508580
4713 actual error = -0.0000000000315785
4715 Note that a direct summation of this series would require 10^10 terms
4716 to achieve the same precision as the accelerated sum does in 13 terms.
4719 File: gsl-ref.info, Node: Series Acceleration References, Prev: Example of accelerating a series, Up: Series Acceleration
4721 29.4 References and Further Reading
4722 ===================================
4724 The algorithms used by these functions are described in the following
4727 T. Fessler, W.F. Ford, D.A. Smith, HURRY: An acceleration
4728 algorithm for scalar sequences and series `ACM Transactions on
4729 Mathematical Software', 9(3):346-354, 1983. and Algorithm 602
4732 The theory of the u-transform was presented by Levin,
4734 D. Levin, Development of Non-Linear Transformations for Improving
4735 Convergence of Sequences, `Intern. J. Computer Math.' B3:371-388,
4738 A review paper on the Levin Transform is available online,
4739 Herbert H. H. Homeier, Scalar Levin-Type Sequence Transformations,
4740 `http://arxiv.org/abs/math/0005209'.
4743 File: gsl-ref.info, Node: Wavelet Transforms, Next: Discrete Hankel Transforms, Prev: Series Acceleration, Up: Top
4745 30 Wavelet Transforms
4746 *********************
4748 This chapter describes functions for performing Discrete Wavelet
4749 Transforms (DWTs). The library includes wavelets for real data in both
4750 one and two dimensions. The wavelet functions are declared in the
4751 header files `gsl_wavelet.h' and `gsl_wavelet2d.h'.
4756 * DWT Initialization::
4757 * DWT Transform Functions::
4762 File: gsl-ref.info, Node: DWT Definitions, Next: DWT Initialization, Up: Wavelet Transforms
4767 The continuous wavelet transform and its inverse are defined by the
4770 w(s,\tau) = \int f(t) * \psi^*_{s,\tau}(t) dt
4774 f(t) = \int \int_{-\infty}^\infty w(s, \tau) * \psi_{s,\tau}(t) d\tau ds
4776 where the basis functions \psi_{s,\tau} are obtained by scaling and
4777 translation from a single function, referred to as the "mother wavelet".
4779 The discrete version of the wavelet transform acts on equally-spaced
4780 samples, with fixed scaling and translation steps (s, \tau). The
4781 frequency and time axes are sampled "dyadically" on scales of 2^j
4782 through a level parameter j. The resulting family of functions
4783 {\psi_{j,n}} constitutes an orthonormal basis for square-integrable
4786 The discrete wavelet transform is an O(N) algorithm, and is also
4787 referred to as the "fast wavelet transform".
4790 File: gsl-ref.info, Node: DWT Initialization, Next: DWT Transform Functions, Prev: DWT Definitions, Up: Wavelet Transforms
4795 The `gsl_wavelet' structure contains the filter coefficients defining
4796 the wavelet and any associated offset parameters.
4798 -- Function: gsl_wavelet * gsl_wavelet_alloc (const gsl_wavelet_type *
4800 This function allocates and initializes a wavelet object of type
4801 T. The parameter K selects the specific member of the wavelet
4802 family. A null pointer is returned if insufficient memory is
4803 available or if a unsupported member is selected.
4805 The following wavelet types are implemented:
4807 -- Wavelet: gsl_wavelet_daubechies
4808 -- Wavelet: gsl_wavelet_daubechies_centered
4809 The is the Daubechies wavelet family of maximum phase with k/2
4810 vanishing moments. The implemented wavelets are k=4, 6, ..., 20,
4813 -- Wavelet: gsl_wavelet_haar
4814 -- Wavelet: gsl_wavelet_haar_centered
4815 This is the Haar wavelet. The only valid choice of k for the Haar
4818 -- Wavelet: gsl_wavelet_bspline
4819 -- Wavelet: gsl_wavelet_bspline_centered
4820 This is the biorthogonal B-spline wavelet family of order (i,j).
4821 The implemented values of k = 100*i + j are 103, 105, 202, 204,
4822 206, 208, 301, 303, 305 307, 309.
4824 The centered forms of the wavelets align the coefficients of the various
4825 sub-bands on edges. Thus the resulting visualization of the
4826 coefficients of the wavelet transform in the phase plane is easier to
4829 -- Function: const char * gsl_wavelet_name (const gsl_wavelet * W)
4830 This function returns a pointer to the name of the wavelet family
4833 -- Function: void gsl_wavelet_free (gsl_wavelet * W)
4834 This function frees the wavelet object W.
4836 The `gsl_wavelet_workspace' structure contains scratch space of the
4837 same size as the input data and is used to hold intermediate results
4838 during the transform.
4840 -- Function: gsl_wavelet_workspace * gsl_wavelet_workspace_alloc
4842 This function allocates a workspace for the discrete wavelet
4843 transform. To perform a one-dimensional transform on N elements,
4844 a workspace of size N must be provided. For two-dimensional
4845 transforms of N-by-N matrices it is sufficient to allocate a
4846 workspace of size N, since the transform operates on individual
4849 -- Function: void gsl_wavelet_workspace_free (gsl_wavelet_workspace *
4851 This function frees the allocated workspace WORK.
4854 File: gsl-ref.info, Node: DWT Transform Functions, Next: DWT Examples, Prev: DWT Initialization, Up: Wavelet Transforms
4856 30.3 Transform Functions
4857 ========================
4859 This sections describes the actual functions performing the discrete
4860 wavelet transform. Note that the transforms use periodic boundary
4861 conditions. If the signal is not periodic in the sample length then
4862 spurious coefficients will appear at the beginning and end of each level
4867 * DWT in one dimension::
4868 * DWT in two dimension::
4871 File: gsl-ref.info, Node: DWT in one dimension, Next: DWT in two dimension, Up: DWT Transform Functions
4873 30.3.1 Wavelet transforms in one dimension
4874 ------------------------------------------
4876 -- Function: int gsl_wavelet_transform (const gsl_wavelet * W, double
4877 * DATA, size_t STRIDE, size_t N, gsl_wavelet_direction DIR,
4878 gsl_wavelet_workspace * WORK)
4879 -- Function: int gsl_wavelet_transform_forward (const gsl_wavelet * W,
4880 double * DATA, size_t STRIDE, size_t N, gsl_wavelet_workspace
4882 -- Function: int gsl_wavelet_transform_inverse (const gsl_wavelet * W,
4883 double * DATA, size_t STRIDE, size_t N, gsl_wavelet_workspace
4885 These functions compute in-place forward and inverse discrete
4886 wavelet transforms of length N with stride STRIDE on the array
4887 DATA. The length of the transform N is restricted to powers of
4888 two. For the `transform' version of the function the argument DIR
4889 can be either `forward' (+1) or `backward' (-1). A workspace WORK
4890 of length N must be provided.
4892 For the forward transform, the elements of the original array are
4893 replaced by the discrete wavelet transform f_i -> w_{j,k} in a
4894 packed triangular storage layout, where J is the index of the level
4895 j = 0 ... J-1 and K is the index of the coefficient within each
4896 level, k = 0 ... (2^j)-1. The total number of levels is J =
4897 \log_2(n). The output data has the following form,
4899 (s_{-1,0}, d_{0,0}, d_{1,0}, d_{1,1}, d_{2,0}, ...,
4900 d_{j,k}, ..., d_{J-1,2^{J-1}-1})
4902 where the first element is the smoothing coefficient s_{-1,0},
4903 followed by the detail coefficients d_{j,k} for each level j. The
4904 backward transform inverts these coefficients to obtain the
4907 These functions return a status of `GSL_SUCCESS' upon successful
4908 completion. `GSL_EINVAL' is returned if N is not an integer power
4909 of 2 or if insufficient workspace is provided.
4912 File: gsl-ref.info, Node: DWT in two dimension, Prev: DWT in one dimension, Up: DWT Transform Functions
4914 30.3.2 Wavelet transforms in two dimension
4915 ------------------------------------------
4917 The library provides functions to perform two-dimensional discrete
4918 wavelet transforms on square matrices. The matrix dimensions must be an
4919 integer power of two. There are two possible orderings of the rows and
4920 columns in the two-dimensional wavelet transform, referred to as the
4921 "standard" and "non-standard" forms.
4923 The "standard" transform performs a complete discrete wavelet
4924 transform on the rows of the matrix, followed by a separate complete
4925 discrete wavelet transform on the columns of the resulting
4926 row-transformed matrix. This procedure uses the same ordering as a
4927 two-dimensional fourier transform.
4929 The "non-standard" transform is performed in interleaved passes on
4930 the rows and columns of the matrix for each level of the transform. The
4931 first level of the transform is applied to the matrix rows, and then to
4932 the matrix columns. This procedure is then repeated across the rows and
4933 columns of the data for the subsequent levels of the transform, until
4934 the full discrete wavelet transform is complete. The non-standard form
4935 of the discrete wavelet transform is typically used in image analysis.
4937 The functions described in this section are declared in the header
4938 file `gsl_wavelet2d.h'.
4940 -- Function: int gsl_wavelet2d_transform (const gsl_wavelet * W,
4941 double * DATA, size_t TDA, size_t SIZE1, size_t SIZE2,
4942 gsl_wavelet_direction DIR, gsl_wavelet_workspace * WORK)
4943 -- Function: int gsl_wavelet2d_transform_forward (const gsl_wavelet *
4944 W, double * DATA, size_t TDA, size_t SIZE1, size_t SIZE2,
4945 gsl_wavelet_workspace * WORK)
4946 -- Function: int gsl_wavelet2d_transform_inverse (const gsl_wavelet *
4947 W, double * DATA, size_t TDA, size_t SIZE1, size_t SIZE2,
4948 gsl_wavelet_workspace * WORK)
4949 These functions compute two-dimensional in-place forward and
4950 inverse discrete wavelet transforms in standard and non-standard
4951 forms on the array DATA stored in row-major form with dimensions
4952 SIZE1 and SIZE2 and physical row length TDA. The dimensions must
4953 be equal (square matrix) and are restricted to powers of two. For
4954 the `transform' version of the function the argument DIR can be
4955 either `forward' (+1) or `backward' (-1). A workspace WORK of the
4956 appropriate size must be provided. On exit, the appropriate
4957 elements of the array DATA are replaced by their two-dimensional
4960 The functions return a status of `GSL_SUCCESS' upon successful
4961 completion. `GSL_EINVAL' is returned if SIZE1 and SIZE2 are not
4962 equal and integer powers of 2, or if insufficient workspace is
4965 -- Function: int gsl_wavelet2d_transform_matrix (const gsl_wavelet *
4966 W, gsl_matrix * M, gsl_wavelet_direction DIR,
4967 gsl_wavelet_workspace * WORK)
4968 -- Function: int gsl_wavelet2d_transform_matrix_forward (const
4969 gsl_wavelet * W, gsl_matrix * M, gsl_wavelet_workspace * WORK)
4970 -- Function: int gsl_wavelet2d_transform_matrix_inverse (const
4971 gsl_wavelet * W, gsl_matrix * M, gsl_wavelet_workspace * WORK)
4972 These functions compute the two-dimensional in-place wavelet
4973 transform on a matrix A.
4975 -- Function: int gsl_wavelet2d_nstransform (const gsl_wavelet * W,
4976 double * DATA, size_t TDA, size_t SIZE1, size_t SIZE2,
4977 gsl_wavelet_direction DIR, gsl_wavelet_workspace * WORK)
4978 -- Function: int gsl_wavelet2d_nstransform_forward (const gsl_wavelet
4979 * W, double * DATA, size_t TDA, size_t SIZE1, size_t SIZE2,
4980 gsl_wavelet_workspace * WORK)
4981 -- Function: int gsl_wavelet2d_nstransform_inverse (const gsl_wavelet
4982 * W, double * DATA, size_t TDA, size_t SIZE1, size_t SIZE2,
4983 gsl_wavelet_workspace * WORK)
4984 These functions compute the two-dimensional wavelet transform in
4987 -- Function: int gsl_wavelet2d_nstransform_matrix (const gsl_wavelet *
4988 W, gsl_matrix * M, gsl_wavelet_direction DIR,
4989 gsl_wavelet_workspace * WORK)
4990 -- Function: int gsl_wavelet2d_nstransform_matrix_forward (const
4991 gsl_wavelet * W, gsl_matrix * M, gsl_wavelet_workspace * WORK)
4992 -- Function: int gsl_wavelet2d_nstransform_matrix_inverse (const
4993 gsl_wavelet * W, gsl_matrix * M, gsl_wavelet_workspace * WORK)
4994 These functions compute the non-standard form of the
4995 two-dimensional in-place wavelet transform on a matrix A.
4998 File: gsl-ref.info, Node: DWT Examples, Next: DWT References, Prev: DWT Transform Functions, Up: Wavelet Transforms
5003 The following program demonstrates the use of the one-dimensional
5004 wavelet transform functions. It computes an approximation to an input
5005 signal (of length 256) using the 20 largest components of the wavelet
5006 transform, while setting the others to zero.
5010 #include <gsl/gsl_sort.h>
5011 #include <gsl/gsl_wavelet.h>
5014 main (int argc, char **argv)
5016 int i, n = 256, nc = 20;
5017 double *data = malloc (n * sizeof (double));
5018 double *abscoeff = malloc (n * sizeof (double));
5019 size_t *p = malloc (n * sizeof (size_t));
5023 gsl_wavelet_workspace *work;
5025 w = gsl_wavelet_alloc (gsl_wavelet_daubechies, 4);
5026 work = gsl_wavelet_workspace_alloc (n);
5028 f = fopen (argv[1], "r");
5029 for (i = 0; i < n; i++)
5031 fscanf (f, "%lg", &data[i]);
5035 gsl_wavelet_transform_forward (w, data, 1, n, work);
5037 for (i = 0; i < n; i++)
5039 abscoeff[i] = fabs (data[i]);
5042 gsl_sort_index (p, abscoeff, 1, n);
5044 for (i = 0; (i + nc) < n; i++)
5047 gsl_wavelet_transform_inverse (w, data, 1, n, work);
5049 for (i = 0; i < n; i++)
5051 printf ("%g\n", data[i]);
5054 gsl_wavelet_free (w);
5055 gsl_wavelet_workspace_free (work);
5063 The output can be used with the GNU plotutils `graph' program,
5065 $ ./a.out ecg.dat > dwt.dat
5066 $ graph -T ps -x 0 256 32 -h 0.3 -a dwt.dat > dwt.ps
5069 File: gsl-ref.info, Node: DWT References, Prev: DWT Examples, Up: Wavelet Transforms
5071 30.5 References and Further Reading
5072 ===================================
5074 The mathematical background to wavelet transforms is covered in the
5075 original lectures by Daubechies,
5077 Ingrid Daubechies. Ten Lectures on Wavelets. `CBMS-NSF Regional
5078 Conference Series in Applied Mathematics' (1992), SIAM, ISBN
5081 An easy to read introduction to the subject with an emphasis on the
5082 application of the wavelet transform in various branches of science is,
5084 Paul S. Addison. `The Illustrated Wavelet Transform Handbook'.
5085 Institute of Physics Publishing (2002), ISBN 0750306920.
5087 For extensive coverage of signal analysis by wavelets, wavelet packets
5088 and local cosine bases see,
5090 S. G. Mallat. `A wavelet tour of signal processing' (Second
5091 edition). Academic Press (1999), ISBN 012466606X.
5093 The concept of multiresolution analysis underlying the wavelet transform
5096 S. G. Mallat. Multiresolution Approximations and Wavelet
5097 Orthonormal Bases of L^2(R). `Transactions of the American
5098 Mathematical Society', 315(1), 1989, 69-87.
5100 S. G. Mallat. A Theory for Multiresolution Signal
5101 Decomposition--The Wavelet Representation. `IEEE Transactions on
5102 Pattern Analysis and Machine Intelligence', 11, 1989, 674-693.
5104 The coefficients for the individual wavelet families implemented by the
5105 library can be found in the following papers,
5107 I. Daubechies. Orthonormal Bases of Compactly Supported Wavelets.
5108 `Communications on Pure and Applied Mathematics', 41 (1988)
5111 A. Cohen, I. Daubechies, and J.-C. Feauveau. Biorthogonal Bases
5112 of Compactly Supported Wavelets. `Communications on Pure and
5113 Applied Mathematics', 45 (1992) 485-560.
5115 The PhysioNet archive of physiological datasets can be found online at
5116 `http://www.physionet.org/' and is described in the following paper,
5118 Goldberger et al. PhysioBank, PhysioToolkit, and PhysioNet:
5119 Components of a New Research Resource for Complex Physiologic
5120 Signals. `Circulation' 101(23):e215-e220 2000.
5123 File: gsl-ref.info, Node: Discrete Hankel Transforms, Next: One dimensional Root-Finding, Prev: Wavelet Transforms, Up: Top
5125 31 Discrete Hankel Transforms
5126 *****************************
5128 This chapter describes functions for performing Discrete Hankel
5129 Transforms (DHTs). The functions are declared in the header file
5134 * Discrete Hankel Transform Definition::
5135 * Discrete Hankel Transform Functions::
5136 * Discrete Hankel Transform References::
5139 File: gsl-ref.info, Node: Discrete Hankel Transform Definition, Next: Discrete Hankel Transform Functions, Up: Discrete Hankel Transforms
5144 The discrete Hankel transform acts on a vector of sampled data, where
5145 the samples are assumed to have been taken at points related to the
5146 zeroes of a Bessel function of fixed order; compare this to the case of
5147 the discrete Fourier transform, where samples are taken at points
5148 related to the zeroes of the sine or cosine function.
5150 Specifically, let f(t) be a function on the unit interval. Then the
5151 finite \nu-Hankel transform of f(t) is defined to be the set of numbers
5153 g_m = \int_0^1 t dt J_\nu(j_(\nu,m)t) f(t),
5156 f(t) = \sum_{m=1}^\infty (2 J_\nu(j_(\nu,m)x) / J_(\nu+1)(j_(\nu,m))^2) g_m.
5158 Suppose that f is band-limited in the sense that g_m=0 for m > M. Then
5159 we have the following fundamental sampling theorem.
5160 g_m = (2 / j_(\nu,M)^2)
5161 \sum_{k=1}^{M-1} f(j_(\nu,k)/j_(\nu,M))
5162 (J_\nu(j_(\nu,m) j_(\nu,k) / j_(\nu,M)) / J_(\nu+1)(j_(\nu,k))^2).
5164 It is this discrete expression which defines the discrete Hankel
5165 transform. The kernel in the summation above defines the matrix of the
5166 \nu-Hankel transform of size M-1. The coefficients of this matrix,
5167 being dependent on \nu and M, must be precomputed and stored; the
5168 `gsl_dht' object encapsulates this data. The allocation function
5169 `gsl_dht_alloc' returns a `gsl_dht' object which must be properly
5170 initialized with `gsl_dht_init' before it can be used to perform
5171 transforms on data sample vectors, for fixed \nu and M, using the
5172 `gsl_dht_apply' function. The implementation allows a scaling of the
5173 fundamental interval, for convenience, so that one can assume the
5174 function is defined on the interval [0,X], rather than the unit
5177 Notice that by assumption f(t) vanishes at the endpoints of the
5178 interval, consistent with the inversion formula and the sampling
5179 formula given above. Therefore, this transform corresponds to an
5180 orthogonal expansion in eigenfunctions of the Dirichlet problem for the
5181 Bessel differential equation.
5184 File: gsl-ref.info, Node: Discrete Hankel Transform Functions, Next: Discrete Hankel Transform References, Prev: Discrete Hankel Transform Definition, Up: Discrete Hankel Transforms
5189 -- Function: gsl_dht * gsl_dht_alloc (size_t SIZE)
5190 This function allocates a Discrete Hankel transform object of size
5193 -- Function: int gsl_dht_init (gsl_dht * T, double NU, double XMAX)
5194 This function initializes the transform T for the given values of
5197 -- Function: gsl_dht * gsl_dht_new (size_t SIZE, double NU, double
5199 This function allocates a Discrete Hankel transform object of size
5200 SIZE and initializes it for the given values of NU and X.
5202 -- Function: void gsl_dht_free (gsl_dht * T)
5203 This function frees the transform T.
5205 -- Function: int gsl_dht_apply (const gsl_dht * T, double * F_IN,
5207 This function applies the transform T to the array F_IN whose size
5208 is equal to the size of the transform. The result is stored in
5209 the array F_OUT which must be of the same length.
5211 -- Function: double gsl_dht_x_sample (const gsl_dht * T, int N)
5212 This function returns the value of the N-th sample point in the
5213 unit interval, (j_{\nu,n+1}/j_{\nu,M}) X. These are the points
5214 where the function f(t) is assumed to be sampled.
5216 -- Function: double gsl_dht_k_sample (const gsl_dht * T, int N)
5217 This function returns the value of the N-th sample point in
5218 "k-space", j_{\nu,n+1}/X.
5221 File: gsl-ref.info, Node: Discrete Hankel Transform References, Prev: Discrete Hankel Transform Functions, Up: Discrete Hankel Transforms
5223 31.3 References and Further Reading
5224 ===================================
5226 The algorithms used by these functions are described in the following
5229 H. Fisk Johnson, Comp. Phys. Comm. 43, 181 (1987).
5231 D. Lemoine, J. Chem. Phys. 101, 3936 (1994).
5234 File: gsl-ref.info, Node: One dimensional Root-Finding, Next: One dimensional Minimization, Prev: Discrete Hankel Transforms, Up: Top
5236 32 One dimensional Root-Finding
5237 *******************************
5239 This chapter describes routines for finding roots of arbitrary
5240 one-dimensional functions. The library provides low level components
5241 for a variety of iterative solvers and convergence tests. These can be
5242 combined by the user to achieve the desired solution, with full access
5243 to the intermediate steps of the iteration. Each class of methods uses
5244 the same framework, so that you can switch between solvers at runtime
5245 without needing to recompile your program. Each instance of a solver
5246 keeps track of its own state, allowing the solvers to be used in
5247 multi-threaded programs.
5249 The header file `gsl_roots.h' contains prototypes for the root
5250 finding functions and related declarations.
5254 * Root Finding Overview::
5255 * Root Finding Caveats::
5256 * Initializing the Solver::
5257 * Providing the function to solve::
5258 * Search Bounds and Guesses::
5259 * Root Finding Iteration::
5260 * Search Stopping Parameters::
5261 * Root Bracketing Algorithms::
5262 * Root Finding Algorithms using Derivatives::
5263 * Root Finding Examples::
5264 * Root Finding References and Further Reading::
5267 File: gsl-ref.info, Node: Root Finding Overview, Next: Root Finding Caveats, Up: One dimensional Root-Finding
5272 One-dimensional root finding algorithms can be divided into two classes,
5273 "root bracketing" and "root polishing". Algorithms which proceed by
5274 bracketing a root are guaranteed to converge. Bracketing algorithms
5275 begin with a bounded region known to contain a root. The size of this
5276 bounded region is reduced, iteratively, until it encloses the root to a
5277 desired tolerance. This provides a rigorous error estimate for the
5278 location of the root.
5280 The technique of "root polishing" attempts to improve an initial
5281 guess to the root. These algorithms converge only if started "close
5282 enough" to a root, and sacrifice a rigorous error bound for speed. By
5283 approximating the behavior of a function in the vicinity of a root they
5284 attempt to find a higher order improvement of an initial guess. When
5285 the behavior of the function is compatible with the algorithm and a good
5286 initial guess is available a polishing algorithm can provide rapid
5289 In GSL both types of algorithm are available in similar frameworks.
5290 The user provides a high-level driver for the algorithms, and the
5291 library provides the individual functions necessary for each of the
5292 steps. There are three main phases of the iteration. The steps are,
5294 * initialize solver state, S, for algorithm T
5296 * update S using the iteration T
5298 * test S for convergence, and repeat iteration if necessary
5300 The state for bracketing solvers is held in a `gsl_root_fsolver'
5301 struct. The updating procedure uses only function evaluations (not
5302 derivatives). The state for root polishing solvers is held in a
5303 `gsl_root_fdfsolver' struct. The updates require both the function and
5304 its derivative (hence the name `fdf') to be supplied by the user.
5307 File: gsl-ref.info, Node: Root Finding Caveats, Next: Initializing the Solver, Prev: Root Finding Overview, Up: One dimensional Root-Finding
5312 Note that root finding functions can only search for one root at a time.
5313 When there are several roots in the search area, the first root to be
5314 found will be returned; however it is difficult to predict which of the
5315 roots this will be. _In most cases, no error will be reported if you
5316 try to find a root in an area where there is more than one._
5318 Care must be taken when a function may have a multiple root (such as
5319 f(x) = (x-x_0)^2 or f(x) = (x-x_0)^3). It is not possible to use
5320 root-bracketing algorithms on even-multiplicity roots. For these
5321 algorithms the initial interval must contain a zero-crossing, where the
5322 function is negative at one end of the interval and positive at the
5323 other end. Roots with even-multiplicity do not cross zero, but only
5324 touch it instantaneously. Algorithms based on root bracketing will
5325 still work for odd-multiplicity roots (e.g. cubic, quintic, ...). Root
5326 polishing algorithms generally work with higher multiplicity roots, but
5327 at a reduced rate of convergence. In these cases the "Steffenson
5328 algorithm" can be used to accelerate the convergence of multiple roots.
5330 While it is not absolutely required that f have a root within the
5331 search region, numerical root finding functions should not be used
5332 haphazardly to check for the _existence_ of roots. There are better
5333 ways to do this. Because it is easy to create situations where
5334 numerical root finders can fail, it is a bad idea to throw a root
5335 finder at a function you do not know much about. In general it is best
5336 to examine the function visually by plotting before searching for a
5340 File: gsl-ref.info, Node: Initializing the Solver, Next: Providing the function to solve, Prev: Root Finding Caveats, Up: One dimensional Root-Finding
5342 32.3 Initializing the Solver
5343 ============================
5345 -- Function: gsl_root_fsolver * gsl_root_fsolver_alloc (const
5346 gsl_root_fsolver_type * T)
5347 This function returns a pointer to a newly allocated instance of a
5348 solver of type T. For example, the following code creates an
5349 instance of a bisection solver,
5351 const gsl_root_fsolver_type * T
5352 = gsl_root_fsolver_bisection;
5353 gsl_root_fsolver * s
5354 = gsl_root_fsolver_alloc (T);
5356 If there is insufficient memory to create the solver then the
5357 function returns a null pointer and the error handler is invoked
5358 with an error code of `GSL_ENOMEM'.
5360 -- Function: gsl_root_fdfsolver * gsl_root_fdfsolver_alloc (const
5361 gsl_root_fdfsolver_type * T)
5362 This function returns a pointer to a newly allocated instance of a
5363 derivative-based solver of type T. For example, the following
5364 code creates an instance of a Newton-Raphson solver,
5366 const gsl_root_fdfsolver_type * T
5367 = gsl_root_fdfsolver_newton;
5368 gsl_root_fdfsolver * s
5369 = gsl_root_fdfsolver_alloc (T);
5371 If there is insufficient memory to create the solver then the
5372 function returns a null pointer and the error handler is invoked
5373 with an error code of `GSL_ENOMEM'.
5375 -- Function: int gsl_root_fsolver_set (gsl_root_fsolver * S,
5376 gsl_function * F, double X_LOWER, double X_UPPER)
5377 This function initializes, or reinitializes, an existing solver S
5378 to use the function F and the initial search interval [X_LOWER,
5381 -- Function: int gsl_root_fdfsolver_set (gsl_root_fdfsolver * S,
5382 gsl_function_fdf * FDF, double ROOT)
5383 This function initializes, or reinitializes, an existing solver S
5384 to use the function and derivative FDF and the initial guess ROOT.
5386 -- Function: void gsl_root_fsolver_free (gsl_root_fsolver * S)
5387 -- Function: void gsl_root_fdfsolver_free (gsl_root_fdfsolver * S)
5388 These functions free all the memory associated with the solver S.
5390 -- Function: const char * gsl_root_fsolver_name (const
5391 gsl_root_fsolver * S)
5392 -- Function: const char * gsl_root_fdfsolver_name (const
5393 gsl_root_fdfsolver * S)
5394 These functions return a pointer to the name of the solver. For
5397 printf ("s is a '%s' solver\n",
5398 gsl_root_fsolver_name (s));
5400 would print something like `s is a 'bisection' solver'.
5403 File: gsl-ref.info, Node: Providing the function to solve, Next: Search Bounds and Guesses, Prev: Initializing the Solver, Up: One dimensional Root-Finding
5405 32.4 Providing the function to solve
5406 ====================================
5408 You must provide a continuous function of one variable for the root
5409 finders to operate on, and, sometimes, its first derivative. In order
5410 to allow for general parameters the functions are defined by the
5411 following data types:
5413 -- Data Type: gsl_function
5414 This data type defines a general function with parameters.
5416 `double (* function) (double X, void * PARAMS)'
5417 this function should return the value f(x,params) for
5418 argument X and parameters PARAMS
5421 a pointer to the parameters of the function
5423 Here is an example for the general quadratic function,
5425 f(x) = a x^2 + b x + c
5427 with a = 3, b = 2, c = 1. The following code defines a `gsl_function'
5428 `F' which you could pass to a root finder:
5430 struct my_f_params { double a; double b; double c; };
5433 my_f (double x, void * p) {
5434 struct my_f_params * params
5435 = (struct my_f_params *)p;
5436 double a = (params->a);
5437 double b = (params->b);
5438 double c = (params->c);
5440 return (a * x + b) * x + c;
5444 struct my_f_params params = { 3.0, 2.0, 1.0 };
5449 The function f(x) can be evaluated using the following macro,
5451 #define GSL_FN_EVAL(F,x)
5452 (*((F)->function))(x,(F)->params)
5454 -- Data Type: gsl_function_fdf
5455 This data type defines a general function with parameters and its
5458 `double (* f) (double X, void * PARAMS)'
5459 this function should return the value of f(x,params) for
5460 argument X and parameters PARAMS
5462 `double (* df) (double X, void * PARAMS)'
5463 this function should return the value of the derivative of F
5464 with respect to X, f'(x,params), for argument X and
5467 `void (* fdf) (double X, void * PARAMS, double * F, double * Df)'
5468 this function should set the values of the function F to
5469 f(x,params) and its derivative DF to f'(x,params) for
5470 argument X and parameters PARAMS. This function provides an
5471 optimization of the separate functions for f(x) and f'(x)--it
5472 is always faster to compute the function and its derivative
5476 a pointer to the parameters of the function
5478 Here is an example where f(x) = 2\exp(2x):
5481 my_f (double x, void * params)
5487 my_df (double x, void * params)
5489 return 2 * exp (2 * x);
5493 my_fdf (double x, void * params,
5494 double * f, double * df)
5496 double t = exp (2 * x);
5499 *df = 2 * t; /* uses existing value */
5502 gsl_function_fdf FDF;
5509 The function f(x) can be evaluated using the following macro,
5511 #define GSL_FN_FDF_EVAL_F(FDF,x)
5512 (*((FDF)->f))(x,(FDF)->params)
5514 The derivative f'(x) can be evaluated using the following macro,
5516 #define GSL_FN_FDF_EVAL_DF(FDF,x)
5517 (*((FDF)->df))(x,(FDF)->params)
5519 and both the function y = f(x) and its derivative dy = f'(x) can be
5520 evaluated at the same time using the following macro,
5522 #define GSL_FN_FDF_EVAL_F_DF(FDF,x,y,dy)
5523 (*((FDF)->fdf))(x,(FDF)->params,(y),(dy))
5525 The macro stores f(x) in its Y argument and f'(x) in its DY
5526 argument--both of these should be pointers to `double'.
5529 File: gsl-ref.info, Node: Search Bounds and Guesses, Next: Root Finding Iteration, Prev: Providing the function to solve, Up: One dimensional Root-Finding
5531 32.5 Search Bounds and Guesses
5532 ==============================
5534 You provide either search bounds or an initial guess; this section
5535 explains how search bounds and guesses work and how function arguments
5538 A guess is simply an x value which is iterated until it is within
5539 the desired precision of a root. It takes the form of a `double'.
5541 Search bounds are the endpoints of a interval which is iterated until
5542 the length of the interval is smaller than the requested precision. The
5543 interval is defined by two values, the lower limit and the upper limit.
5544 Whether the endpoints are intended to be included in the interval or not
5545 depends on the context in which the interval is used.
5548 File: gsl-ref.info, Node: Root Finding Iteration, Next: Search Stopping Parameters, Prev: Search Bounds and Guesses, Up: One dimensional Root-Finding
5553 The following functions drive the iteration of each algorithm. Each
5554 function performs one iteration to update the state of any solver of the
5555 corresponding type. The same functions work for all solvers so that
5556 different methods can be substituted at runtime without modifications to
5559 -- Function: int gsl_root_fsolver_iterate (gsl_root_fsolver * S)
5560 -- Function: int gsl_root_fdfsolver_iterate (gsl_root_fdfsolver * S)
5561 These functions perform a single iteration of the solver S. If the
5562 iteration encounters an unexpected problem then an error code will
5566 the iteration encountered a singular point where the function
5567 or its derivative evaluated to `Inf' or `NaN'.
5570 the derivative of the function vanished at the iteration
5571 point, preventing the algorithm from continuing without a
5574 The solver maintains a current best estimate of the root at all
5575 times. The bracketing solvers also keep track of the current best
5576 interval bounding the root. This information can be accessed with the
5577 following auxiliary functions,
5579 -- Function: double gsl_root_fsolver_root (const gsl_root_fsolver * S)
5580 -- Function: double gsl_root_fdfsolver_root (const gsl_root_fdfsolver
5582 These functions return the current estimate of the root for the
5585 -- Function: double gsl_root_fsolver_x_lower (const gsl_root_fsolver *
5587 -- Function: double gsl_root_fsolver_x_upper (const gsl_root_fsolver *
5589 These functions return the current bracketing interval for the
5593 File: gsl-ref.info, Node: Search Stopping Parameters, Next: Root Bracketing Algorithms, Prev: Root Finding Iteration, Up: One dimensional Root-Finding
5595 32.7 Search Stopping Parameters
5596 ===============================
5598 A root finding procedure should stop when one of the following
5601 * A root has been found to within the user-specified precision.
5603 * A user-specified maximum number of iterations has been reached.
5605 * An error has occurred.
5607 The handling of these conditions is under user control. The functions
5608 below allow the user to test the precision of the current result in
5609 several standard ways.
5611 -- Function: int gsl_root_test_interval (double X_LOWER, double
5612 X_UPPER, double EPSABS, double EPSREL)
5613 This function tests for the convergence of the interval [X_LOWER,
5614 X_UPPER] with absolute error EPSABS and relative error EPSREL.
5615 The test returns `GSL_SUCCESS' if the following condition is
5618 |a - b| < epsabs + epsrel min(|a|,|b|)
5620 when the interval x = [a,b] does not include the origin. If the
5621 interval includes the origin then \min(|a|,|b|) is replaced by
5622 zero (which is the minimum value of |x| over the interval). This
5623 ensures that the relative error is accurately estimated for roots
5624 close to the origin.
5626 This condition on the interval also implies that any estimate of
5627 the root r in the interval satisfies the same condition with
5628 respect to the true root r^*,
5630 |r - r^*| < epsabs + epsrel r^*
5632 assuming that the true root r^* is contained within the interval.
5634 -- Function: int gsl_root_test_delta (double X1, double X0, double
5635 EPSABS, double EPSREL)
5636 This function tests for the convergence of the sequence ..., X0,
5637 X1 with absolute error EPSABS and relative error EPSREL. The test
5638 returns `GSL_SUCCESS' if the following condition is achieved,
5640 |x_1 - x_0| < epsabs + epsrel |x_1|
5642 and returns `GSL_CONTINUE' otherwise.
5644 -- Function: int gsl_root_test_residual (double F, double EPSABS)
5645 This function tests the residual value F against the absolute
5646 error bound EPSABS. The test returns `GSL_SUCCESS' if the
5647 following condition is achieved,
5651 and returns `GSL_CONTINUE' otherwise. This criterion is suitable
5652 for situations where the precise location of the root, x, is
5653 unimportant provided a value can be found where the residual,
5654 |f(x)|, is small enough.
5657 File: gsl-ref.info, Node: Root Bracketing Algorithms, Next: Root Finding Algorithms using Derivatives, Prev: Search Stopping Parameters, Up: One dimensional Root-Finding
5659 32.8 Root Bracketing Algorithms
5660 ===============================
5662 The root bracketing algorithms described in this section require an
5663 initial interval which is guaranteed to contain a root--if a and b are
5664 the endpoints of the interval then f(a) must differ in sign from f(b).
5665 This ensures that the function crosses zero at least once in the
5666 interval. If a valid initial interval is used then these algorithm
5667 cannot fail, provided the function is well-behaved.
5669 Note that a bracketing algorithm cannot find roots of even degree,
5670 since these do not cross the x-axis.
5672 -- Solver: gsl_root_fsolver_bisection
5673 The "bisection algorithm" is the simplest method of bracketing the
5674 roots of a function. It is the slowest algorithm provided by the
5675 library, with linear convergence.
5677 On each iteration, the interval is bisected and the value of the
5678 function at the midpoint is calculated. The sign of this value is
5679 used to determine which half of the interval does not contain a
5680 root. That half is discarded to give a new, smaller interval
5681 containing the root. This procedure can be continued indefinitely
5682 until the interval is sufficiently small.
5684 At any time the current estimate of the root is taken as the
5685 midpoint of the interval.
5688 -- Solver: gsl_root_fsolver_falsepos
5689 The "false position algorithm" is a method of finding roots based
5690 on linear interpolation. Its convergence is linear, but it is
5691 usually faster than bisection.
5693 On each iteration a line is drawn between the endpoints (a,f(a))
5694 and (b,f(b)) and the point where this line crosses the x-axis
5695 taken as a "midpoint". The value of the function at this point is
5696 calculated and its sign is used to determine which side of the
5697 interval does not contain a root. That side is discarded to give a
5698 new, smaller interval containing the root. This procedure can be
5699 continued indefinitely until the interval is sufficiently small.
5701 The best estimate of the root is taken from the linear
5702 interpolation of the interval on the current iteration.
5705 -- Solver: gsl_root_fsolver_brent
5706 The "Brent-Dekker method" (referred to here as "Brent's method")
5707 combines an interpolation strategy with the bisection algorithm.
5708 This produces a fast algorithm which is still robust.
5710 On each iteration Brent's method approximates the function using an
5711 interpolating curve. On the first iteration this is a linear
5712 interpolation of the two endpoints. For subsequent iterations the
5713 algorithm uses an inverse quadratic fit to the last three points,
5714 for higher accuracy. The intercept of the interpolating curve
5715 with the x-axis is taken as a guess for the root. If it lies
5716 within the bounds of the current interval then the interpolating
5717 point is accepted, and used to generate a smaller interval. If
5718 the interpolating point is not accepted then the algorithm falls
5719 back to an ordinary bisection step.
5721 The best estimate of the root is taken from the most recent
5722 interpolation or bisection.
5725 File: gsl-ref.info, Node: Root Finding Algorithms using Derivatives, Next: Root Finding Examples, Prev: Root Bracketing Algorithms, Up: One dimensional Root-Finding
5727 32.9 Root Finding Algorithms using Derivatives
5728 ==============================================
5730 The root polishing algorithms described in this section require an
5731 initial guess for the location of the root. There is no absolute
5732 guarantee of convergence--the function must be suitable for this
5733 technique and the initial guess must be sufficiently close to the root
5734 for it to work. When these conditions are satisfied then convergence is
5737 These algorithms make use of both the function and its derivative.
5739 -- Derivative Solver: gsl_root_fdfsolver_newton
5740 Newton's Method is the standard root-polishing algorithm. The
5741 algorithm begins with an initial guess for the location of the
5742 root. On each iteration, a line tangent to the function f is
5743 drawn at that position. The point where this line crosses the
5744 x-axis becomes the new guess. The iteration is defined by the
5747 x_{i+1} = x_i - f(x_i)/f'(x_i)
5749 Newton's method converges quadratically for single roots, and
5750 linearly for multiple roots.
5753 -- Derivative Solver: gsl_root_fdfsolver_secant
5754 The "secant method" is a simplified version of Newton's method
5755 which does not require the computation of the derivative on every
5758 On its first iteration the algorithm begins with Newton's method,
5759 using the derivative to compute a first step,
5761 x_1 = x_0 - f(x_0)/f'(x_0)
5763 Subsequent iterations avoid the evaluation of the derivative by
5764 replacing it with a numerical estimate, the slope of the line
5765 through the previous two points,
5767 x_{i+1} = x_i f(x_i) / f'_{est} where
5768 f'_{est} = (f(x_i) - f(x_{i-1})/(x_i - x_{i-1})
5770 When the derivative does not change significantly in the vicinity
5771 of the root the secant method gives a useful saving.
5772 Asymptotically the secant method is faster than Newton's method
5773 whenever the cost of evaluating the derivative is more than 0.44
5774 times the cost of evaluating the function itself. As with all
5775 methods of computing a numerical derivative the estimate can
5776 suffer from cancellation errors if the separation of the points
5779 On single roots, the method has a convergence of order (1 + \sqrt
5780 5)/2 (approximately 1.62). It converges linearly for multiple
5784 -- Derivative Solver: gsl_root_fdfsolver_steffenson
5785 The "Steffenson Method" provides the fastest convergence of all the
5786 routines. It combines the basic Newton algorithm with an Aitken
5787 "delta-squared" acceleration. If the Newton iterates are x_i then
5788 the acceleration procedure generates a new sequence R_i,
5790 R_i = x_i - (x_{i+1} - x_i)^2 / (x_{i+2} - 2 x_{i+1} + x_{i})
5792 which converges faster than the original sequence under reasonable
5793 conditions. The new sequence requires three terms before it can
5794 produce its first value so the method returns accelerated values
5795 on the second and subsequent iterations. On the first iteration
5796 it returns the ordinary Newton estimate. The Newton iterate is
5797 also returned if the denominator of the acceleration term ever
5800 As with all acceleration procedures this method can become
5801 unstable if the function is not well-behaved.
5804 File: gsl-ref.info, Node: Root Finding Examples, Next: Root Finding References and Further Reading, Prev: Root Finding Algorithms using Derivatives, Up: One dimensional Root-Finding
5809 For any root finding algorithm we need to prepare the function to be
5810 solved. For this example we will use the general quadratic equation
5811 described earlier. We first need a header file (`demo_fn.h') to define
5812 the function parameters,
5814 struct quadratic_params
5819 double quadratic (double x, void *params);
5820 double quadratic_deriv (double x, void *params);
5821 void quadratic_fdf (double x, void *params,
5822 double *y, double *dy);
5824 We place the function definitions in a separate file (`demo_fn.c'),
5827 quadratic (double x, void *params)
5829 struct quadratic_params *p
5830 = (struct quadratic_params *) params;
5836 return (a * x + b) * x + c;
5840 quadratic_deriv (double x, void *params)
5842 struct quadratic_params *p
5843 = (struct quadratic_params *) params;
5849 return 2.0 * a * x + b;
5853 quadratic_fdf (double x, void *params,
5854 double *y, double *dy)
5856 struct quadratic_params *p
5857 = (struct quadratic_params *) params;
5863 *y = (a * x + b) * x + c;
5864 *dy = 2.0 * a * x + b;
5867 The first program uses the function solver `gsl_root_fsolver_brent' for
5868 Brent's method and the general quadratic defined above to solve the
5873 with solution x = \sqrt 5 = 2.236068...
5876 #include <gsl/gsl_errno.h>
5877 #include <gsl/gsl_math.h>
5878 #include <gsl/gsl_roots.h>
5880 #include "demo_fn.h"
5881 #include "demo_fn.c"
5887 int iter = 0, max_iter = 100;
5888 const gsl_root_fsolver_type *T;
5889 gsl_root_fsolver *s;
5890 double r = 0, r_expected = sqrt (5.0);
5891 double x_lo = 0.0, x_hi = 5.0;
5893 struct quadratic_params params = {1.0, 0.0, -5.0};
5895 F.function = &quadratic;
5898 T = gsl_root_fsolver_brent;
5899 s = gsl_root_fsolver_alloc (T);
5900 gsl_root_fsolver_set (s, &F, x_lo, x_hi);
5902 printf ("using %s method\n",
5903 gsl_root_fsolver_name (s));
5905 printf ("%5s [%9s, %9s] %9s %10s %9s\n",
5906 "iter", "lower", "upper", "root",
5912 status = gsl_root_fsolver_iterate (s);
5913 r = gsl_root_fsolver_root (s);
5914 x_lo = gsl_root_fsolver_x_lower (s);
5915 x_hi = gsl_root_fsolver_x_upper (s);
5916 status = gsl_root_test_interval (x_lo, x_hi,
5919 if (status == GSL_SUCCESS)
5920 printf ("Converged:\n");
5922 printf ("%5d [%.7f, %.7f] %.7f %+.7f %.7f\n",
5927 while (status == GSL_CONTINUE && iter < max_iter);
5929 gsl_root_fsolver_free (s);
5934 Here are the results of the iterations,
5938 iter [ lower, upper] root err err(est)
5939 1 [1.0000000, 5.0000000] 1.0000000 -1.2360680 4.0000000
5940 2 [1.0000000, 3.0000000] 3.0000000 +0.7639320 2.0000000
5941 3 [2.0000000, 3.0000000] 2.0000000 -0.2360680 1.0000000
5942 4 [2.2000000, 3.0000000] 2.2000000 -0.0360680 0.8000000
5943 5 [2.2000000, 2.2366300] 2.2366300 +0.0005621 0.0366300
5945 6 [2.2360634, 2.2366300] 2.2360634 -0.0000046 0.0005666
5947 If the program is modified to use the bisection solver instead of
5948 Brent's method, by changing `gsl_root_fsolver_brent' to
5949 `gsl_root_fsolver_bisection' the slower convergence of the Bisection
5950 method can be observed,
5953 using bisection method
5954 iter [ lower, upper] root err err(est)
5955 1 [0.0000000, 2.5000000] 1.2500000 -0.9860680 2.5000000
5956 2 [1.2500000, 2.5000000] 1.8750000 -0.3610680 1.2500000
5957 3 [1.8750000, 2.5000000] 2.1875000 -0.0485680 0.6250000
5958 4 [2.1875000, 2.5000000] 2.3437500 +0.1076820 0.3125000
5959 5 [2.1875000, 2.3437500] 2.2656250 +0.0295570 0.1562500
5960 6 [2.1875000, 2.2656250] 2.2265625 -0.0095055 0.0781250
5961 7 [2.2265625, 2.2656250] 2.2460938 +0.0100258 0.0390625
5962 8 [2.2265625, 2.2460938] 2.2363281 +0.0002601 0.0195312
5963 9 [2.2265625, 2.2363281] 2.2314453 -0.0046227 0.0097656
5964 10 [2.2314453, 2.2363281] 2.2338867 -0.0021813 0.0048828
5965 11 [2.2338867, 2.2363281] 2.2351074 -0.0009606 0.0024414
5967 12 [2.2351074, 2.2363281] 2.2357178 -0.0003502 0.0012207
5969 The next program solves the same function using a derivative solver
5973 #include <gsl/gsl_errno.h>
5974 #include <gsl/gsl_math.h>
5975 #include <gsl/gsl_roots.h>
5977 #include "demo_fn.h"
5978 #include "demo_fn.c"
5984 int iter = 0, max_iter = 100;
5985 const gsl_root_fdfsolver_type *T;
5986 gsl_root_fdfsolver *s;
5987 double x0, x = 5.0, r_expected = sqrt (5.0);
5988 gsl_function_fdf FDF;
5989 struct quadratic_params params = {1.0, 0.0, -5.0};
5992 FDF.df = &quadratic_deriv;
5993 FDF.fdf = &quadratic_fdf;
5994 FDF.params = ¶ms;
5996 T = gsl_root_fdfsolver_newton;
5997 s = gsl_root_fdfsolver_alloc (T);
5998 gsl_root_fdfsolver_set (s, &FDF, x);
6000 printf ("using %s method\n",
6001 gsl_root_fdfsolver_name (s));
6003 printf ("%-5s %10s %10s %10s\n",
6004 "iter", "root", "err", "err(est)");
6008 status = gsl_root_fdfsolver_iterate (s);
6010 x = gsl_root_fdfsolver_root (s);
6011 status = gsl_root_test_delta (x, x0, 0, 1e-3);
6013 if (status == GSL_SUCCESS)
6014 printf ("Converged:\n");
6016 printf ("%5d %10.7f %+10.7f %10.7f\n",
6017 iter, x, x - r_expected, x - x0);
6019 while (status == GSL_CONTINUE && iter < max_iter);
6021 gsl_root_fdfsolver_free (s);
6025 Here are the results for Newton's method,
6029 iter root err err(est)
6030 1 3.0000000 +0.7639320 -2.0000000
6031 2 2.3333333 +0.0972654 -0.6666667
6032 3 2.2380952 +0.0020273 -0.0952381
6034 4 2.2360689 +0.0000009 -0.0020263
6036 Note that the error can be estimated more accurately by taking the
6037 difference between the current iterate and next iterate rather than the
6038 previous iterate. The other derivative solvers can be investigated by
6039 changing `gsl_root_fdfsolver_newton' to `gsl_root_fdfsolver_secant' or
6040 `gsl_root_fdfsolver_steffenson'.
6043 File: gsl-ref.info, Node: Root Finding References and Further Reading, Prev: Root Finding Examples, Up: One dimensional Root-Finding
6045 32.11 References and Further Reading
6046 ====================================
6048 For information on the Brent-Dekker algorithm see the following two
6051 R. P. Brent, "An algorithm with guaranteed convergence for finding
6052 a zero of a function", `Computer Journal', 14 (1971) 422-425
6054 J. C. P. Bus and T. J. Dekker, "Two Efficient Algorithms with
6055 Guaranteed Convergence for Finding a Zero of a Function", `ACM
6056 Transactions of Mathematical Software', Vol. 1 No. 4 (1975) 330-345
6059 File: gsl-ref.info, Node: One dimensional Minimization, Next: Multidimensional Root-Finding, Prev: One dimensional Root-Finding, Up: Top
6061 33 One dimensional Minimization
6062 *******************************
6064 This chapter describes routines for finding minima of arbitrary
6065 one-dimensional functions. The library provides low level components
6066 for a variety of iterative minimizers and convergence tests. These can
6067 be combined by the user to achieve the desired solution, with full
6068 access to the intermediate steps of the algorithms. Each class of
6069 methods uses the same framework, so that you can switch between
6070 minimizers at runtime without needing to recompile your program. Each
6071 instance of a minimizer keeps track of its own state, allowing the
6072 minimizers to be used in multi-threaded programs.
6074 The header file `gsl_min.h' contains prototypes for the minimization
6075 functions and related declarations. To use the minimization algorithms
6076 to find the maximum of a function simply invert its sign.
6080 * Minimization Overview::
6081 * Minimization Caveats::
6082 * Initializing the Minimizer::
6083 * Providing the function to minimize::
6084 * Minimization Iteration::
6085 * Minimization Stopping Parameters::
6086 * Minimization Algorithms::
6087 * Minimization Examples::
6088 * Minimization References and Further Reading::
6091 File: gsl-ref.info, Node: Minimization Overview, Next: Minimization Caveats, Up: One dimensional Minimization
6096 The minimization algorithms begin with a bounded region known to contain
6097 a minimum. The region is described by a lower bound a and an upper
6098 bound b, with an estimate of the location of the minimum x.
6100 The value of the function at x must be less than the value of the
6101 function at the ends of the interval,
6105 This condition guarantees that a minimum is contained somewhere within
6106 the interval. On each iteration a new point x' is selected using one
6107 of the available algorithms. If the new point is a better estimate of
6108 the minimum, i.e. where f(x') < f(x), then the current estimate of the
6109 minimum x is updated. The new point also allows the size of the
6110 bounded interval to be reduced, by choosing the most compact set of
6111 points which satisfies the constraint f(a) > f(x) < f(b). The interval
6112 is reduced until it encloses the true minimum to a desired tolerance.
6113 This provides a best estimate of the location of the minimum and a
6114 rigorous error estimate.
6116 Several bracketing algorithms are available within a single
6117 framework. The user provides a high-level driver for the algorithm,
6118 and the library provides the individual functions necessary for each of
6119 the steps. There are three main phases of the iteration. The steps
6122 * initialize minimizer state, S, for algorithm T
6124 * update S using the iteration T
6126 * test S for convergence, and repeat iteration if necessary
6128 The state for the minimizers is held in a `gsl_min_fminimizer' struct.
6129 The updating procedure uses only function evaluations (not derivatives).
6132 File: gsl-ref.info, Node: Minimization Caveats, Next: Initializing the Minimizer, Prev: Minimization Overview, Up: One dimensional Minimization
6137 Note that minimization functions can only search for one minimum at a
6138 time. When there are several minima in the search area, the first
6139 minimum to be found will be returned; however it is difficult to predict
6140 which of the minima this will be. _In most cases, no error will be
6141 reported if you try to find a minimum in an area where there is more
6144 With all minimization algorithms it can be difficult to determine the
6145 location of the minimum to full numerical precision. The behavior of
6146 the function in the region of the minimum x^* can be approximated by a
6149 y = f(x^*) + (1/2) f''(x^*) (x - x^*)^2
6151 and the second term of this expansion can be lost when added to the
6152 first term at finite precision. This magnifies the error in locating
6153 x^*, making it proportional to \sqrt \epsilon (where \epsilon is the
6154 relative accuracy of the floating point numbers). For functions with
6155 higher order minima, such as x^4, the magnification of the error is
6156 correspondingly worse. The best that can be achieved is to converge to
6157 the limit of numerical accuracy in the function values, rather than the
6158 location of the minimum itself.
6161 File: gsl-ref.info, Node: Initializing the Minimizer, Next: Providing the function to minimize, Prev: Minimization Caveats, Up: One dimensional Minimization
6163 33.3 Initializing the Minimizer
6164 ===============================
6166 -- Function: gsl_min_fminimizer * gsl_min_fminimizer_alloc (const
6167 gsl_min_fminimizer_type * T)
6168 This function returns a pointer to a newly allocated instance of a
6169 minimizer of type T. For example, the following code creates an
6170 instance of a golden section minimizer,
6172 const gsl_min_fminimizer_type * T
6173 = gsl_min_fminimizer_goldensection;
6174 gsl_min_fminimizer * s
6175 = gsl_min_fminimizer_alloc (T);
6177 If there is insufficient memory to create the minimizer then the
6178 function returns a null pointer and the error handler is invoked
6179 with an error code of `GSL_ENOMEM'.
6181 -- Function: int gsl_min_fminimizer_set (gsl_min_fminimizer * S,
6182 gsl_function * F, double X_MINIMUM, double X_LOWER, double
6184 This function sets, or resets, an existing minimizer S to use the
6185 function F and the initial search interval [X_LOWER, X_UPPER],
6186 with a guess for the location of the minimum X_MINIMUM.
6188 If the interval given does not contain a minimum, then the function
6189 returns an error code of `GSL_EINVAL'.
6191 -- Function: int gsl_min_fminimizer_set_with_values
6192 (gsl_min_fminimizer * S, gsl_function * F, double X_MINIMUM,
6193 double F_MINIMUM, double X_LOWER, double F_LOWER, double
6194 X_UPPER, double F_UPPER)
6195 This function is equivalent to `gsl_min_fminimizer_set' but uses
6196 the values F_MINIMUM, F_LOWER and F_UPPER instead of computing
6197 `f(x_minimum)', `f(x_lower)' and `f(x_upper)'.
6199 -- Function: void gsl_min_fminimizer_free (gsl_min_fminimizer * S)
6200 This function frees all the memory associated with the minimizer S.
6202 -- Function: const char * gsl_min_fminimizer_name (const
6203 gsl_min_fminimizer * S)
6204 This function returns a pointer to the name of the minimizer. For
6207 printf ("s is a '%s' minimizer\n",
6208 gsl_min_fminimizer_name (s));
6210 would print something like `s is a 'brent' minimizer'.
6213 File: gsl-ref.info, Node: Providing the function to minimize, Next: Minimization Iteration, Prev: Initializing the Minimizer, Up: One dimensional Minimization
6215 33.4 Providing the function to minimize
6216 =======================================
6218 You must provide a continuous function of one variable for the
6219 minimizers to operate on. In order to allow for general parameters the
6220 functions are defined by a `gsl_function' data type (*note Providing
6221 the function to solve::).
6224 File: gsl-ref.info, Node: Minimization Iteration, Next: Minimization Stopping Parameters, Prev: Providing the function to minimize, Up: One dimensional Minimization
6229 The following functions drive the iteration of each algorithm. Each
6230 function performs one iteration to update the state of any minimizer of
6231 the corresponding type. The same functions work for all minimizers so
6232 that different methods can be substituted at runtime without
6233 modifications to the code.
6235 -- Function: int gsl_min_fminimizer_iterate (gsl_min_fminimizer * S)
6236 This function performs a single iteration of the minimizer S. If
6237 the iteration encounters an unexpected problem then an error code
6241 the iteration encountered a singular point where the function
6242 evaluated to `Inf' or `NaN'.
6245 the algorithm could not improve the current best
6246 approximation or bounding interval.
6248 The minimizer maintains a current best estimate of the position of
6249 the minimum at all times, and the current interval bounding the minimum.
6250 This information can be accessed with the following auxiliary functions,
6252 -- Function: double gsl_min_fminimizer_x_minimum (const
6253 gsl_min_fminimizer * S)
6254 This function returns the current estimate of the position of the
6255 minimum for the minimizer S.
6257 -- Function: double gsl_min_fminimizer_x_upper (const
6258 gsl_min_fminimizer * S)
6259 -- Function: double gsl_min_fminimizer_x_lower (const
6260 gsl_min_fminimizer * S)
6261 These functions return the current upper and lower bound of the
6262 interval for the minimizer S.
6264 -- Function: double gsl_min_fminimizer_f_minimum (const
6265 gsl_min_fminimizer * S)
6266 -- Function: double gsl_min_fminimizer_f_upper (const
6267 gsl_min_fminimizer * S)
6268 -- Function: double gsl_min_fminimizer_f_lower (const
6269 gsl_min_fminimizer * S)
6270 These functions return the value of the function at the current
6271 estimate of the minimum and at the upper and lower bounds of the
6272 interval for the minimizer S.
6275 File: gsl-ref.info, Node: Minimization Stopping Parameters, Next: Minimization Algorithms, Prev: Minimization Iteration, Up: One dimensional Minimization
6277 33.6 Stopping Parameters
6278 ========================
6280 A minimization procedure should stop when one of the following
6283 * A minimum has been found to within the user-specified precision.
6285 * A user-specified maximum number of iterations has been reached.
6287 * An error has occurred.
6289 The handling of these conditions is under user control. The function
6290 below allows the user to test the precision of the current result.
6292 -- Function: int gsl_min_test_interval (double X_LOWER, double
6293 X_UPPER, double EPSABS, double EPSREL)
6294 This function tests for the convergence of the interval [X_LOWER,
6295 X_UPPER] with absolute error EPSABS and relative error EPSREL.
6296 The test returns `GSL_SUCCESS' if the following condition is
6299 |a - b| < epsabs + epsrel min(|a|,|b|)
6301 when the interval x = [a,b] does not include the origin. If the
6302 interval includes the origin then \min(|a|,|b|) is replaced by
6303 zero (which is the minimum value of |x| over the interval). This
6304 ensures that the relative error is accurately estimated for minima
6305 close to the origin.
6307 This condition on the interval also implies that any estimate of
6308 the minimum x_m in the interval satisfies the same condition with
6309 respect to the true minimum x_m^*,
6311 |x_m - x_m^*| < epsabs + epsrel x_m^*
6313 assuming that the true minimum x_m^* is contained within the
6317 File: gsl-ref.info, Node: Minimization Algorithms, Next: Minimization Examples, Prev: Minimization Stopping Parameters, Up: One dimensional Minimization
6319 33.7 Minimization Algorithms
6320 ============================
6322 The minimization algorithms described in this section require an initial
6323 interval which is guaranteed to contain a minimum--if a and b are the
6324 endpoints of the interval and x is an estimate of the minimum then f(a)
6325 > f(x) < f(b). This ensures that the function has at least one minimum
6326 somewhere in the interval. If a valid initial interval is used then
6327 these algorithm cannot fail, provided the function is well-behaved.
6329 -- Minimizer: gsl_min_fminimizer_goldensection
6330 The "golden section algorithm" is the simplest method of bracketing
6331 the minimum of a function. It is the slowest algorithm provided
6332 by the library, with linear convergence.
6334 On each iteration, the algorithm first compares the subintervals
6335 from the endpoints to the current minimum. The larger subinterval
6336 is divided in a golden section (using the famous ratio (3-\sqrt
6337 5)/2 = 0.3189660...) and the value of the function at this new
6338 point is calculated. The new value is used with the constraint
6339 f(a') > f(x') < f(b') to a select new interval containing the
6340 minimum, by discarding the least useful point. This procedure can
6341 be continued indefinitely until the interval is sufficiently
6342 small. Choosing the golden section as the bisection ratio can be
6343 shown to provide the fastest convergence for this type of
6347 -- Minimizer: gsl_min_fminimizer_brent
6348 The "Brent minimization algorithm" combines a parabolic
6349 interpolation with the golden section algorithm. This produces a
6350 fast algorithm which is still robust.
6352 The outline of the algorithm can be summarized as follows: on each
6353 iteration Brent's method approximates the function using an
6354 interpolating parabola through three existing points. The minimum
6355 of the parabola is taken as a guess for the minimum. If it lies
6356 within the bounds of the current interval then the interpolating
6357 point is accepted, and used to generate a smaller interval. If
6358 the interpolating point is not accepted then the algorithm falls
6359 back to an ordinary golden section step. The full details of
6360 Brent's method include some additional checks to improve
6364 File: gsl-ref.info, Node: Minimization Examples, Next: Minimization References and Further Reading, Prev: Minimization Algorithms, Up: One dimensional Minimization
6369 The following program uses the Brent algorithm to find the minimum of
6370 the function f(x) = \cos(x) + 1, which occurs at x = \pi. The starting
6371 interval is (0,6), with an initial guess for the minimum of 2.
6374 #include <gsl/gsl_errno.h>
6375 #include <gsl/gsl_math.h>
6376 #include <gsl/gsl_min.h>
6378 double fn1 (double x, void * params)
6380 return cos(x) + 1.0;
6387 int iter = 0, max_iter = 100;
6388 const gsl_min_fminimizer_type *T;
6389 gsl_min_fminimizer *s;
6390 double m = 2.0, m_expected = M_PI;
6391 double a = 0.0, b = 6.0;
6397 T = gsl_min_fminimizer_brent;
6398 s = gsl_min_fminimizer_alloc (T);
6399 gsl_min_fminimizer_set (s, &F, m, a, b);
6401 printf ("using %s method\n",
6402 gsl_min_fminimizer_name (s));
6404 printf ("%5s [%9s, %9s] %9s %10s %9s\n",
6405 "iter", "lower", "upper", "min",
6408 printf ("%5d [%.7f, %.7f] %.7f %+.7f %.7f\n",
6410 m, m - m_expected, b - a);
6415 status = gsl_min_fminimizer_iterate (s);
6417 m = gsl_min_fminimizer_x_minimum (s);
6418 a = gsl_min_fminimizer_x_lower (s);
6419 b = gsl_min_fminimizer_x_upper (s);
6422 = gsl_min_test_interval (a, b, 0.001, 0.0);
6424 if (status == GSL_SUCCESS)
6425 printf ("Converged:\n");
6427 printf ("%5d [%.7f, %.7f] "
6428 "%.7f %+.7f %.7f\n",
6430 m, m - m_expected, b - a);
6432 while (status == GSL_CONTINUE && iter < max_iter);
6434 gsl_min_fminimizer_free (s);
6439 Here are the results of the minimization procedure.
6442 0 [0.0000000, 6.0000000] 2.0000000 -1.1415927 6.0000000
6443 1 [2.0000000, 6.0000000] 3.2758640 +0.1342713 4.0000000
6444 2 [2.0000000, 3.2831929] 3.2758640 +0.1342713 1.2831929
6445 3 [2.8689068, 3.2831929] 3.2758640 +0.1342713 0.4142862
6446 4 [2.8689068, 3.2831929] 3.2758640 +0.1342713 0.4142862
6447 5 [2.8689068, 3.2758640] 3.1460585 +0.0044658 0.4069572
6448 6 [3.1346075, 3.2758640] 3.1460585 +0.0044658 0.1412565
6449 7 [3.1346075, 3.1874620] 3.1460585 +0.0044658 0.0528545
6450 8 [3.1346075, 3.1460585] 3.1460585 +0.0044658 0.0114510
6451 9 [3.1346075, 3.1460585] 3.1424060 +0.0008133 0.0114510
6452 10 [3.1346075, 3.1424060] 3.1415885 -0.0000041 0.0077985
6454 11 [3.1415885, 3.1424060] 3.1415927 -0.0000000 0.0008175
6457 File: gsl-ref.info, Node: Minimization References and Further Reading, Prev: Minimization Examples, Up: One dimensional Minimization
6459 33.9 References and Further Reading
6460 ===================================
6462 Further information on Brent's algorithm is available in the following
6465 Richard Brent, `Algorithms for minimization without derivatives',
6466 Prentice-Hall (1973), republished by Dover in paperback (2002),
6470 File: gsl-ref.info, Node: Multidimensional Root-Finding, Next: Multidimensional Minimization, Prev: One dimensional Minimization, Up: Top
6472 34 Multidimensional Root-Finding
6473 ********************************
6475 This chapter describes functions for multidimensional root-finding
6476 (solving nonlinear systems with n equations in n unknowns). The
6477 library provides low level components for a variety of iterative
6478 solvers and convergence tests. These can be combined by the user to
6479 achieve the desired solution, with full access to the intermediate
6480 steps of the iteration. Each class of methods uses the same framework,
6481 so that you can switch between solvers at runtime without needing to
6482 recompile your program. Each instance of a solver keeps track of its
6483 own state, allowing the solvers to be used in multi-threaded programs.
6484 The solvers are based on the original Fortran library MINPACK.
6486 The header file `gsl_multiroots.h' contains prototypes for the
6487 multidimensional root finding functions and related declarations.
6491 * Overview of Multidimensional Root Finding::
6492 * Initializing the Multidimensional Solver::
6493 * Providing the multidimensional system of equations to solve::
6494 * Iteration of the multidimensional solver::
6495 * Search Stopping Parameters for the multidimensional solver::
6496 * Algorithms using Derivatives::
6497 * Algorithms without Derivatives::
6498 * Example programs for Multidimensional Root finding::
6499 * References and Further Reading for Multidimensional Root Finding::
6502 File: gsl-ref.info, Node: Overview of Multidimensional Root Finding, Next: Initializing the Multidimensional Solver, Up: Multidimensional Root-Finding
6507 The problem of multidimensional root finding requires the simultaneous
6508 solution of n equations, f_i, in n variables, x_i,
6510 f_i (x_1, ..., x_n) = 0 for i = 1 ... n.
6512 In general there are no bracketing methods available for n dimensional
6513 systems, and no way of knowing whether any solutions exist. All
6514 algorithms proceed from an initial guess using a variant of the Newton
6517 x -> x' = x - J^{-1} f(x)
6519 where x, f are vector quantities and J is the Jacobian matrix J_{ij} =
6520 d f_i / d x_j. Additional strategies can be used to enlarge the region
6521 of convergence. These include requiring a decrease in the norm |f| on
6522 each step proposed by Newton's method, or taking steepest-descent steps
6523 in the direction of the negative gradient of |f|.
6525 Several root-finding algorithms are available within a single
6526 framework. The user provides a high-level driver for the algorithms,
6527 and the library provides the individual functions necessary for each of
6528 the steps. There are three main phases of the iteration. The steps
6531 * initialize solver state, S, for algorithm T
6533 * update S using the iteration T
6535 * test S for convergence, and repeat iteration if necessary
6537 The evaluation of the Jacobian matrix can be problematic, either because
6538 programming the derivatives is intractable or because computation of the
6539 n^2 terms of the matrix becomes too expensive. For these reasons the
6540 algorithms provided by the library are divided into two classes
6541 according to whether the derivatives are available or not.
6543 The state for solvers with an analytic Jacobian matrix is held in a
6544 `gsl_multiroot_fdfsolver' struct. The updating procedure requires both
6545 the function and its derivatives to be supplied by the user.
6547 The state for solvers which do not use an analytic Jacobian matrix is
6548 held in a `gsl_multiroot_fsolver' struct. The updating procedure uses
6549 only function evaluations (not derivatives). The algorithms estimate
6550 the matrix J or J^{-1} by approximate methods.
6553 File: gsl-ref.info, Node: Initializing the Multidimensional Solver, Next: Providing the multidimensional system of equations to solve, Prev: Overview of Multidimensional Root Finding, Up: Multidimensional Root-Finding
6555 34.2 Initializing the Solver
6556 ============================
6558 The following functions initialize a multidimensional solver, either
6559 with or without derivatives. The solver itself depends only on the
6560 dimension of the problem and the algorithm and can be reused for
6563 -- Function: gsl_multiroot_fsolver * gsl_multiroot_fsolver_alloc
6564 (const gsl_multiroot_fsolver_type * T, size_t N)
6565 This function returns a pointer to a newly allocated instance of a
6566 solver of type T for a system of N dimensions. For example, the
6567 following code creates an instance of a hybrid solver, to solve a
6568 3-dimensional system of equations.
6570 const gsl_multiroot_fsolver_type * T
6571 = gsl_multiroot_fsolver_hybrid;
6572 gsl_multiroot_fsolver * s
6573 = gsl_multiroot_fsolver_alloc (T, 3);
6575 If there is insufficient memory to create the solver then the
6576 function returns a null pointer and the error handler is invoked
6577 with an error code of `GSL_ENOMEM'.
6579 -- Function: gsl_multiroot_fdfsolver * gsl_multiroot_fdfsolver_alloc
6580 (const gsl_multiroot_fdfsolver_type * T, size_t N)
6581 This function returns a pointer to a newly allocated instance of a
6582 derivative solver of type T for a system of N dimensions. For
6583 example, the following code creates an instance of a
6584 Newton-Raphson solver, for a 2-dimensional system of equations.
6586 const gsl_multiroot_fdfsolver_type * T
6587 = gsl_multiroot_fdfsolver_newton;
6588 gsl_multiroot_fdfsolver * s =
6589 gsl_multiroot_fdfsolver_alloc (T, 2);
6591 If there is insufficient memory to create the solver then the
6592 function returns a null pointer and the error handler is invoked
6593 with an error code of `GSL_ENOMEM'.
6595 -- Function: int gsl_multiroot_fsolver_set (gsl_multiroot_fsolver * S,
6596 gsl_multiroot_function * F, const gsl_vector * X)
6597 -- Function: int gsl_multiroot_fdfsolver_set (gsl_multiroot_fdfsolver
6598 * S, gsl_multiroot_function_fdf * FDF, const gsl_vector * X)
6599 These functions set, or reset, an existing solver S to use the
6600 function F or function and derivative FDF, and the initial guess
6601 X. Note that the initial position is copied from X, this argument
6602 is not modified by subsequent iterations.
6604 -- Function: void gsl_multiroot_fsolver_free (gsl_multiroot_fsolver *
6606 -- Function: void gsl_multiroot_fdfsolver_free
6607 (gsl_multiroot_fdfsolver * S)
6608 These functions free all the memory associated with the solver S.
6610 -- Function: const char * gsl_multiroot_fsolver_name (const
6611 gsl_multiroot_fsolver * S)
6612 -- Function: const char * gsl_multiroot_fdfsolver_name (const
6613 gsl_multiroot_fdfsolver * S)
6614 These functions return a pointer to the name of the solver. For
6617 printf ("s is a '%s' solver\n",
6618 gsl_multiroot_fdfsolver_name (s));
6620 would print something like `s is a 'newton' solver'.
6623 File: gsl-ref.info, Node: Providing the multidimensional system of equations to solve, Next: Iteration of the multidimensional solver, Prev: Initializing the Multidimensional Solver, Up: Multidimensional Root-Finding
6625 34.3 Providing the function to solve
6626 ====================================
6628 You must provide n functions of n variables for the root finders to
6629 operate on. In order to allow for general parameters the functions are
6630 defined by the following data types:
6632 -- Data Type: gsl_multiroot_function
6633 This data type defines a general system of functions with
6636 `int (* f) (const gsl_vector * X, void * PARAMS, gsl_vector * F)'
6637 this function should store the vector result f(x,params) in F
6638 for argument X and parameters PARAMS, returning an
6639 appropriate error code if the function cannot be computed.
6642 the dimension of the system, i.e. the number of components of
6643 the vectors X and F.
6646 a pointer to the parameters of the function.
6648 Here is an example using Powell's test function,
6650 f_1(x) = A x_0 x_1 - 1,
6651 f_2(x) = exp(-x_0) + exp(-x_1) - (1 + 1/A)
6653 with A = 10^4. The following code defines a `gsl_multiroot_function'
6654 system `F' which you could pass to a solver:
6656 struct powell_params { double A; };
6659 powell (gsl_vector * x, void * p, gsl_vector * f) {
6660 struct powell_params * params
6661 = *(struct powell_params *)p;
6662 const double A = (params->A);
6663 const double x0 = gsl_vector_get(x,0);
6664 const double x1 = gsl_vector_get(x,1);
6666 gsl_vector_set (f, 0, A * x0 * x1 - 1);
6667 gsl_vector_set (f, 1, (exp(-x0) + exp(-x1)
6672 gsl_multiroot_function F;
6673 struct powell_params params = { 10000.0 };
6679 -- Data Type: gsl_multiroot_function_fdf
6680 This data type defines a general system of functions with
6681 parameters and the corresponding Jacobian matrix of derivatives,
6683 `int (* f) (const gsl_vector * X, void * PARAMS, gsl_vector * F)'
6684 this function should store the vector result f(x,params) in F
6685 for argument X and parameters PARAMS, returning an
6686 appropriate error code if the function cannot be computed.
6688 `int (* df) (const gsl_vector * X, void * PARAMS, gsl_matrix * J)'
6689 this function should store the N-by-N matrix result J_ij = d
6690 f_i(x,params) / d x_j in J for argument X and parameters
6691 PARAMS, returning an appropriate error code if the function
6694 `int (* fdf) (const gsl_vector * X, void * PARAMS, gsl_vector * F, gsl_matrix * J)'
6695 This function should set the values of the F and J as above,
6696 for arguments X and parameters PARAMS. This function
6697 provides an optimization of the separate functions for f(x)
6698 and J(x)--it is always faster to compute the function and its
6699 derivative at the same time.
6702 the dimension of the system, i.e. the number of components of
6703 the vectors X and F.
6706 a pointer to the parameters of the function.
6708 The example of Powell's test function defined above can be extended to
6709 include analytic derivatives using the following code,
6712 powell_df (gsl_vector * x, void * p, gsl_matrix * J)
6714 struct powell_params * params
6715 = *(struct powell_params *)p;
6716 const double A = (params->A);
6717 const double x0 = gsl_vector_get(x,0);
6718 const double x1 = gsl_vector_get(x,1);
6719 gsl_matrix_set (J, 0, 0, A * x1);
6720 gsl_matrix_set (J, 0, 1, A * x0);
6721 gsl_matrix_set (J, 1, 0, -exp(-x0));
6722 gsl_matrix_set (J, 1, 1, -exp(-x1));
6727 powell_fdf (gsl_vector * x, void * p,
6728 gsl_matrix * f, gsl_matrix * J) {
6729 struct powell_params * params
6730 = *(struct powell_params *)p;
6731 const double A = (params->A);
6732 const double x0 = gsl_vector_get(x,0);
6733 const double x1 = gsl_vector_get(x,1);
6735 const double u0 = exp(-x0);
6736 const double u1 = exp(-x1);
6738 gsl_vector_set (f, 0, A * x0 * x1 - 1);
6739 gsl_vector_set (f, 1, u0 + u1 - (1 + 1/A));
6741 gsl_matrix_set (J, 0, 0, A * x1);
6742 gsl_matrix_set (J, 0, 1, A * x0);
6743 gsl_matrix_set (J, 1, 0, -u0);
6744 gsl_matrix_set (J, 1, 1, -u1);
6748 gsl_multiroot_function_fdf FDF;
6751 FDF.df = &powell_df;
6752 FDF.fdf = &powell_fdf;
6756 Note that the function `powell_fdf' is able to reuse existing terms
6757 from the function when calculating the Jacobian, thus saving time.
6760 File: gsl-ref.info, Node: Iteration of the multidimensional solver, Next: Search Stopping Parameters for the multidimensional solver, Prev: Providing the multidimensional system of equations to solve, Up: Multidimensional Root-Finding
6765 The following functions drive the iteration of each algorithm. Each
6766 function performs one iteration to update the state of any solver of the
6767 corresponding type. The same functions work for all solvers so that
6768 different methods can be substituted at runtime without modifications to
6771 -- Function: int gsl_multiroot_fsolver_iterate (gsl_multiroot_fsolver
6773 -- Function: int gsl_multiroot_fdfsolver_iterate
6774 (gsl_multiroot_fdfsolver * S)
6775 These functions perform a single iteration of the solver S. If the
6776 iteration encounters an unexpected problem then an error code will
6780 the iteration encountered a singular point where the function
6781 or its derivative evaluated to `Inf' or `NaN'.
6784 the iteration is not making any progress, preventing the
6785 algorithm from continuing.
6787 The solver maintains a current best estimate of the root `s->x' and
6788 its function value `s->f' at all times. This information can be
6789 accessed with the following auxiliary functions,
6791 -- Function: gsl_vector * gsl_multiroot_fsolver_root (const
6792 gsl_multiroot_fsolver * S)
6793 -- Function: gsl_vector * gsl_multiroot_fdfsolver_root (const
6794 gsl_multiroot_fdfsolver * S)
6795 These functions return the current estimate of the root for the
6796 solver S, given by `s->x'.
6798 -- Function: gsl_vector * gsl_multiroot_fsolver_f (const
6799 gsl_multiroot_fsolver * S)
6800 -- Function: gsl_vector * gsl_multiroot_fdfsolver_f (const
6801 gsl_multiroot_fdfsolver * S)
6802 These functions return the function value f(x) at the current
6803 estimate of the root for the solver S, given by `s->f'.
6805 -- Function: gsl_vector * gsl_multiroot_fsolver_dx (const
6806 gsl_multiroot_fsolver * S)
6807 -- Function: gsl_vector * gsl_multiroot_fdfsolver_dx (const
6808 gsl_multiroot_fdfsolver * S)
6809 These functions return the last step dx taken by the solver S,
6813 File: gsl-ref.info, Node: Search Stopping Parameters for the multidimensional solver, Next: Algorithms using Derivatives, Prev: Iteration of the multidimensional solver, Up: Multidimensional Root-Finding
6815 34.5 Search Stopping Parameters
6816 ===============================
6818 A root finding procedure should stop when one of the following
6821 * A multidimensional root has been found to within the
6822 user-specified precision.
6824 * A user-specified maximum number of iterations has been reached.
6826 * An error has occurred.
6828 The handling of these conditions is under user control. The functions
6829 below allow the user to test the precision of the current result in
6830 several standard ways.
6832 -- Function: int gsl_multiroot_test_delta (const gsl_vector * DX,
6833 const gsl_vector * X, double EPSABS, double EPSREL)
6834 This function tests for the convergence of the sequence by
6835 comparing the last step DX with the absolute error EPSABS and
6836 relative error EPSREL to the current position X. The test returns
6837 `GSL_SUCCESS' if the following condition is achieved,
6839 |dx_i| < epsabs + epsrel |x_i|
6841 for each component of X and returns `GSL_CONTINUE' otherwise.
6843 -- Function: int gsl_multiroot_test_residual (const gsl_vector * F,
6845 This function tests the residual value F against the absolute
6846 error bound EPSABS. The test returns `GSL_SUCCESS' if the
6847 following condition is achieved,
6849 \sum_i |f_i| < epsabs
6851 and returns `GSL_CONTINUE' otherwise. This criterion is suitable
6852 for situations where the precise location of the root, x, is
6853 unimportant provided a value can be found where the residual is
6857 File: gsl-ref.info, Node: Algorithms using Derivatives, Next: Algorithms without Derivatives, Prev: Search Stopping Parameters for the multidimensional solver, Up: Multidimensional Root-Finding
6859 34.6 Algorithms using Derivatives
6860 =================================
6862 The root finding algorithms described in this section make use of both
6863 the function and its derivative. They require an initial guess for the
6864 location of the root, but there is no absolute guarantee of
6865 convergence--the function must be suitable for this technique and the
6866 initial guess must be sufficiently close to the root for it to work.
6867 When the conditions are satisfied then convergence is quadratic.
6869 -- Derivative Solver: gsl_multiroot_fdfsolver_hybridsj
6870 This is a modified version of Powell's Hybrid method as
6871 implemented in the HYBRJ algorithm in MINPACK. Minpack was
6872 written by Jorge J. More', Burton S. Garbow and Kenneth E.
6873 Hillstrom. The Hybrid algorithm retains the fast convergence of
6874 Newton's method but will also reduce the residual when Newton's
6875 method is unreliable.
6877 The algorithm uses a generalized trust region to keep each step
6878 under control. In order to be accepted a proposed new position x'
6879 must satisfy the condition |D (x' - x)| < \delta, where D is a
6880 diagonal scaling matrix and \delta is the size of the trust
6881 region. The components of D are computed internally, using the
6882 column norms of the Jacobian to estimate the sensitivity of the
6883 residual to each component of x. This improves the behavior of the
6884 algorithm for badly scaled functions.
6886 On each iteration the algorithm first determines the standard
6887 Newton step by solving the system J dx = - f. If this step falls
6888 inside the trust region it is used as a trial step in the next
6889 stage. If not, the algorithm uses the linear combination of the
6890 Newton and gradient directions which is predicted to minimize the
6891 norm of the function while staying inside the trust region,
6893 dx = - \alpha J^{-1} f(x) - \beta \nabla |f(x)|^2.
6895 This combination of Newton and gradient directions is referred to
6898 The proposed step is now tested by evaluating the function at the
6899 resulting point, x'. If the step reduces the norm of the function
6900 sufficiently then it is accepted and size of the trust region is
6901 increased. If the proposed step fails to improve the solution
6902 then the size of the trust region is decreased and another trial
6905 The speed of the algorithm is increased by computing the changes
6906 to the Jacobian approximately, using a rank-1 update. If two
6907 successive attempts fail to reduce the residual then the full
6908 Jacobian is recomputed. The algorithm also monitors the progress
6909 of the solution and returns an error if several steps fail to make
6913 the iteration is not making any progress, preventing the
6914 algorithm from continuing.
6917 re-evaluations of the Jacobian indicate that the iteration is
6918 not making any progress, preventing the algorithm from
6922 -- Derivative Solver: gsl_multiroot_fdfsolver_hybridj
6923 This algorithm is an unscaled version of `hybridsj'. The steps are
6924 controlled by a spherical trust region |x' - x| < \delta, instead
6925 of a generalized region. This can be useful if the generalized
6926 region estimated by `hybridsj' is inappropriate.
6928 -- Derivative Solver: gsl_multiroot_fdfsolver_newton
6929 Newton's Method is the standard root-polishing algorithm. The
6930 algorithm begins with an initial guess for the location of the
6931 solution. On each iteration a linear approximation to the
6932 function F is used to estimate the step which will zero all the
6933 components of the residual. The iteration is defined by the
6936 x -> x' = x - J^{-1} f(x)
6938 where the Jacobian matrix J is computed from the derivative
6939 functions provided by F. The step dx is obtained by solving the
6944 using LU decomposition.
6946 -- Derivative Solver: gsl_multiroot_fdfsolver_gnewton
6947 This is a modified version of Newton's method which attempts to
6948 improve global convergence by requiring every step to reduce the
6949 Euclidean norm of the residual, |f(x)|. If the Newton step leads
6950 to an increase in the norm then a reduced step of relative size,
6952 t = (\sqrt(1 + 6 r) - 1) / (3 r)
6954 is proposed, with r being the ratio of norms |f(x')|^2/|f(x)|^2.
6955 This procedure is repeated until a suitable step size is found.
6958 File: gsl-ref.info, Node: Algorithms without Derivatives, Next: Example programs for Multidimensional Root finding, Prev: Algorithms using Derivatives, Up: Multidimensional Root-Finding
6960 34.7 Algorithms without Derivatives
6961 ===================================
6963 The algorithms described in this section do not require any derivative
6964 information to be supplied by the user. Any derivatives needed are
6965 approximated by finite differences. Note that if the
6966 finite-differencing step size chosen by these routines is inappropriate,
6967 an explicit user-supplied numerical derivative can always be used with
6968 the algorithms described in the previous section.
6970 -- Solver: gsl_multiroot_fsolver_hybrids
6971 This is a version of the Hybrid algorithm which replaces calls to
6972 the Jacobian function by its finite difference approximation. The
6973 finite difference approximation is computed using
6974 `gsl_multiroots_fdjac' with a relative step size of
6975 `GSL_SQRT_DBL_EPSILON'. Note that this step size will not be
6976 suitable for all problems.
6978 -- Solver: gsl_multiroot_fsolver_hybrid
6979 This is a finite difference version of the Hybrid algorithm without
6982 -- Solver: gsl_multiroot_fsolver_dnewton
6983 The "discrete Newton algorithm" is the simplest method of solving a
6984 multidimensional system. It uses the Newton iteration
6986 x -> x - J^{-1} f(x)
6988 where the Jacobian matrix J is approximated by taking finite
6989 differences of the function F. The approximation scheme used by
6990 this implementation is,
6992 J_{ij} = (f_i(x + \delta_j) - f_i(x)) / \delta_j
6994 where \delta_j is a step of size \sqrt\epsilon |x_j| with \epsilon
6995 being the machine precision (\epsilon \approx 2.22 \times 10^-16).
6996 The order of convergence of Newton's algorithm is quadratic, but
6997 the finite differences require n^2 function evaluations on each
6998 iteration. The algorithm may become unstable if the finite
6999 differences are not a good approximation to the true derivatives.
7001 -- Solver: gsl_multiroot_fsolver_broyden
7002 The "Broyden algorithm" is a version of the discrete Newton
7003 algorithm which attempts to avoids the expensive update of the
7004 Jacobian matrix on each iteration. The changes to the Jacobian
7005 are also approximated, using a rank-1 update,
7007 J^{-1} \to J^{-1} - (J^{-1} df - dx) dx^T J^{-1} / dx^T J^{-1} df
7009 where the vectors dx and df are the changes in x and f. On the
7010 first iteration the inverse Jacobian is estimated using finite
7011 differences, as in the discrete Newton algorithm.
7013 This approximation gives a fast update but is unreliable if the
7014 changes are not small, and the estimate of the inverse Jacobian
7015 becomes worse as time passes. The algorithm has a tendency to
7016 become unstable unless it starts close to the root. The Jacobian
7017 is refreshed if this instability is detected (consult the source
7020 This algorithm is included only for demonstration purposes, and is
7021 not recommended for serious use.
7024 File: gsl-ref.info, Node: Example programs for Multidimensional Root finding, Next: References and Further Reading for Multidimensional Root Finding, Prev: Algorithms without Derivatives, Up: Multidimensional Root-Finding
7029 The multidimensional solvers are used in a similar way to the
7030 one-dimensional root finding algorithms. This first example
7031 demonstrates the `hybrids' scaled-hybrid algorithm, which does not
7032 require derivatives. The program solves the Rosenbrock system of
7035 f_1 (x, y) = a (1 - x)
7036 f_2 (x, y) = b (y - x^2)
7038 with a = 1, b = 10. The solution of this system lies at (x,y) = (1,1)
7041 The first stage of the program is to define the system of equations,
7045 #include <gsl/gsl_vector.h>
7046 #include <gsl/gsl_multiroots.h>
7055 rosenbrock_f (const gsl_vector * x, void *params,
7058 double a = ((struct rparams *) params)->a;
7059 double b = ((struct rparams *) params)->b;
7061 const double x0 = gsl_vector_get (x, 0);
7062 const double x1 = gsl_vector_get (x, 1);
7064 const double y0 = a * (1 - x0);
7065 const double y1 = b * (x1 - x0 * x0);
7067 gsl_vector_set (f, 0, y0);
7068 gsl_vector_set (f, 1, y1);
7073 The main program begins by creating the function object `f', with the
7074 arguments `(x,y)' and parameters `(a,b)'. The solver `s' is initialized
7075 to use this function, with the `hybrids' method.
7080 const gsl_multiroot_fsolver_type *T;
7081 gsl_multiroot_fsolver *s;
7087 struct rparams p = {1.0, 10.0};
7088 gsl_multiroot_function f = {&rosenbrock_f, n, &p};
7090 double x_init[2] = {-10.0, -5.0};
7091 gsl_vector *x = gsl_vector_alloc (n);
7093 gsl_vector_set (x, 0, x_init[0]);
7094 gsl_vector_set (x, 1, x_init[1]);
7096 T = gsl_multiroot_fsolver_hybrids;
7097 s = gsl_multiroot_fsolver_alloc (T, 2);
7098 gsl_multiroot_fsolver_set (s, &f, x);
7100 print_state (iter, s);
7105 status = gsl_multiroot_fsolver_iterate (s);
7107 print_state (iter, s);
7109 if (status) /* check if solver is stuck */
7113 gsl_multiroot_test_residual (s->f, 1e-7);
7115 while (status == GSL_CONTINUE && iter < 1000);
7117 printf ("status = %s\n", gsl_strerror (status));
7119 gsl_multiroot_fsolver_free (s);
7120 gsl_vector_free (x);
7124 Note that it is important to check the return status of each solver
7125 step, in case the algorithm becomes stuck. If an error condition is
7126 detected, indicating that the algorithm cannot proceed, then the error
7127 can be reported to the user, a new starting point chosen or a different
7130 The intermediate state of the solution is displayed by the following
7131 function. The solver state contains the vector `s->x' which is the
7132 current position, and the vector `s->f' with corresponding function
7136 print_state (size_t iter, gsl_multiroot_fsolver * s)
7138 printf ("iter = %3u x = % .3f % .3f "
7139 "f(x) = % .3e % .3e\n",
7141 gsl_vector_get (s->x, 0),
7142 gsl_vector_get (s->x, 1),
7143 gsl_vector_get (s->f, 0),
7144 gsl_vector_get (s->f, 1));
7147 Here are the results of running the program. The algorithm is started at
7148 (-10,-5) far from the solution. Since the solution is hidden in a
7149 narrow valley the earliest steps follow the gradient of the function
7150 downhill, in an attempt to reduce the large value of the residual. Once
7151 the root has been approximately located, on iteration 8, the Newton
7152 behavior takes over and convergence is very rapid.
7154 iter = 0 x = -10.000 -5.000 f(x) = 1.100e+01 -1.050e+03
7155 iter = 1 x = -10.000 -5.000 f(x) = 1.100e+01 -1.050e+03
7156 iter = 2 x = -3.976 24.827 f(x) = 4.976e+00 9.020e+01
7157 iter = 3 x = -3.976 24.827 f(x) = 4.976e+00 9.020e+01
7158 iter = 4 x = -3.976 24.827 f(x) = 4.976e+00 9.020e+01
7159 iter = 5 x = -1.274 -5.680 f(x) = 2.274e+00 -7.302e+01
7160 iter = 6 x = -1.274 -5.680 f(x) = 2.274e+00 -7.302e+01
7161 iter = 7 x = 0.249 0.298 f(x) = 7.511e-01 2.359e+00
7162 iter = 8 x = 0.249 0.298 f(x) = 7.511e-01 2.359e+00
7163 iter = 9 x = 1.000 0.878 f(x) = 1.268e-10 -1.218e+00
7164 iter = 10 x = 1.000 0.989 f(x) = 1.124e-11 -1.080e-01
7165 iter = 11 x = 1.000 1.000 f(x) = 0.000e+00 0.000e+00
7168 Note that the algorithm does not update the location on every
7169 iteration. Some iterations are used to adjust the trust-region
7170 parameter, after trying a step which was found to be divergent, or to
7171 recompute the Jacobian, when poor convergence behavior is detected.
7173 The next example program adds derivative information, in order to
7174 accelerate the solution. There are two derivative functions
7175 `rosenbrock_df' and `rosenbrock_fdf'. The latter computes both the
7176 function and its derivative simultaneously. This allows the
7177 optimization of any common terms. For simplicity we substitute calls to
7178 the separate `f' and `df' functions at this point in the code below.
7181 rosenbrock_df (const gsl_vector * x, void *params,
7184 const double a = ((struct rparams *) params)->a;
7185 const double b = ((struct rparams *) params)->b;
7187 const double x0 = gsl_vector_get (x, 0);
7189 const double df00 = -a;
7190 const double df01 = 0;
7191 const double df10 = -2 * b * x0;
7192 const double df11 = b;
7194 gsl_matrix_set (J, 0, 0, df00);
7195 gsl_matrix_set (J, 0, 1, df01);
7196 gsl_matrix_set (J, 1, 0, df10);
7197 gsl_matrix_set (J, 1, 1, df11);
7203 rosenbrock_fdf (const gsl_vector * x, void *params,
7204 gsl_vector * f, gsl_matrix * J)
7206 rosenbrock_f (x, params, f);
7207 rosenbrock_df (x, params, J);
7212 The main program now makes calls to the corresponding `fdfsolver'
7213 versions of the functions,
7218 const gsl_multiroot_fdfsolver_type *T;
7219 gsl_multiroot_fdfsolver *s;
7225 struct rparams p = {1.0, 10.0};
7226 gsl_multiroot_function_fdf f = {&rosenbrock_f,
7231 double x_init[2] = {-10.0, -5.0};
7232 gsl_vector *x = gsl_vector_alloc (n);
7234 gsl_vector_set (x, 0, x_init[0]);
7235 gsl_vector_set (x, 1, x_init[1]);
7237 T = gsl_multiroot_fdfsolver_gnewton;
7238 s = gsl_multiroot_fdfsolver_alloc (T, n);
7239 gsl_multiroot_fdfsolver_set (s, &f, x);
7241 print_state (iter, s);
7247 status = gsl_multiroot_fdfsolver_iterate (s);
7249 print_state (iter, s);
7254 status = gsl_multiroot_test_residual (s->f, 1e-7);
7256 while (status == GSL_CONTINUE && iter < 1000);
7258 printf ("status = %s\n", gsl_strerror (status));
7260 gsl_multiroot_fdfsolver_free (s);
7261 gsl_vector_free (x);
7265 The addition of derivative information to the `hybrids' solver does not
7266 make any significant difference to its behavior, since it able to
7267 approximate the Jacobian numerically with sufficient accuracy. To
7268 illustrate the behavior of a different derivative solver we switch to
7269 `gnewton'. This is a traditional Newton solver with the constraint that
7270 it scales back its step if the full step would lead "uphill". Here is
7271 the output for the `gnewton' algorithm,
7273 iter = 0 x = -10.000 -5.000 f(x) = 1.100e+01 -1.050e+03
7274 iter = 1 x = -4.231 -65.317 f(x) = 5.231e+00 -8.321e+02
7275 iter = 2 x = 1.000 -26.358 f(x) = -8.882e-16 -2.736e+02
7276 iter = 3 x = 1.000 1.000 f(x) = -2.220e-16 -4.441e-15
7279 The convergence is much more rapid, but takes a wide excursion out to
7280 the point (-4.23,-65.3). This could cause the algorithm to go astray in
7281 a realistic application. The hybrid algorithm follows the downhill
7282 path to the solution more reliably.
7285 File: gsl-ref.info, Node: References and Further Reading for Multidimensional Root Finding, Prev: Example programs for Multidimensional Root finding, Up: Multidimensional Root-Finding
7287 34.9 References and Further Reading
7288 ===================================
7290 The original version of the Hybrid method is described in the following
7293 M.J.D. Powell, "A Hybrid Method for Nonlinear Equations" (Chap 6, p
7294 87-114) and "A Fortran Subroutine for Solving systems of Nonlinear
7295 Algebraic Equations" (Chap 7, p 115-161), in `Numerical Methods for
7296 Nonlinear Algebraic Equations', P. Rabinowitz, editor. Gordon and
7299 The following papers are also relevant to the algorithms described in
7302 J.J. More', M.Y. Cosnard, "Numerical Solution of Nonlinear
7303 Equations", `ACM Transactions on Mathematical Software', Vol 5, No
7306 C.G. Broyden, "A Class of Methods for Solving Nonlinear
7307 Simultaneous Equations", `Mathematics of Computation', Vol 19
7310 J.J. More', B.S. Garbow, K.E. Hillstrom, "Testing Unconstrained
7311 Optimization Software", ACM Transactions on Mathematical Software,
7312 Vol 7, No 1 (1981), p 17-41
7315 File: gsl-ref.info, Node: Multidimensional Minimization, Next: Least-Squares Fitting, Prev: Multidimensional Root-Finding, Up: Top
7317 35 Multidimensional Minimization
7318 ********************************
7320 This chapter describes routines for finding minima of arbitrary
7321 multidimensional functions. The library provides low level components
7322 for a variety of iterative minimizers and convergence tests. These can
7323 be combined by the user to achieve the desired solution, while providing
7324 full access to the intermediate steps of the algorithms. Each class of
7325 methods uses the same framework, so that you can switch between
7326 minimizers at runtime without needing to recompile your program. Each
7327 instance of a minimizer keeps track of its own state, allowing the
7328 minimizers to be used in multi-threaded programs. The minimization
7329 algorithms can be used to maximize a function by inverting its sign.
7331 The header file `gsl_multimin.h' contains prototypes for the
7332 minimization functions and related declarations.
7336 * Multimin Overview::
7337 * Multimin Caveats::
7338 * Initializing the Multidimensional Minimizer::
7339 * Providing a function to minimize::
7340 * Multimin Iteration::
7341 * Multimin Stopping Criteria::
7342 * Multimin Algorithms::
7343 * Multimin Examples::
7344 * Multimin References and Further Reading::
7347 File: gsl-ref.info, Node: Multimin Overview, Next: Multimin Caveats, Up: Multidimensional Minimization
7352 The problem of multidimensional minimization requires finding a point x
7353 such that the scalar function,
7357 takes a value which is lower than at any neighboring point. For smooth
7358 functions the gradient g = \nabla f vanishes at the minimum. In general
7359 there are no bracketing methods available for the minimization of
7360 n-dimensional functions. The algorithms proceed from an initial guess
7361 using a search algorithm which attempts to move in a downhill direction.
7363 Algorithms making use of the gradient of the function perform a
7364 one-dimensional line minimisation along this direction until the lowest
7365 point is found to a suitable tolerance. The search direction is then
7366 updated with local information from the function and its derivatives,
7367 and the whole process repeated until the true n-dimensional minimum is
7370 The Nelder-Mead Simplex algorithm applies a different strategy. It
7371 maintains n+1 trial parameter vectors as the vertices of a
7372 n-dimensional simplex. In each iteration step it tries to improve the
7373 worst vertex by a simple geometrical transformation until the size of
7374 the simplex falls below a given tolerance.
7376 Both types of algorithms use a standard framework. The user provides
7377 a high-level driver for the algorithms, and the library provides the
7378 individual functions necessary for each of the steps. There are three
7379 main phases of the iteration. The steps are,
7381 * initialize minimizer state, S, for algorithm T
7383 * update S using the iteration T
7385 * test S for convergence, and repeat iteration if necessary
7387 Each iteration step consists either of an improvement to the
7388 line-minimisation in the current direction or an update to the search
7389 direction itself. The state for the minimizers is held in a
7390 `gsl_multimin_fdfminimizer' struct or a `gsl_multimin_fminimizer'
7394 File: gsl-ref.info, Node: Multimin Caveats, Next: Initializing the Multidimensional Minimizer, Prev: Multimin Overview, Up: Multidimensional Minimization
7399 Note that the minimization algorithms can only search for one local
7400 minimum at a time. When there are several local minima in the search
7401 area, the first minimum to be found will be returned; however it is
7402 difficult to predict which of the minima this will be. In most cases,
7403 no error will be reported if you try to find a local minimum in an area
7404 where there is more than one.
7406 It is also important to note that the minimization algorithms find
7407 local minima; there is no way to determine whether a minimum is a global
7408 minimum of the function in question.
7411 File: gsl-ref.info, Node: Initializing the Multidimensional Minimizer, Next: Providing a function to minimize, Prev: Multimin Caveats, Up: Multidimensional Minimization
7413 35.3 Initializing the Multidimensional Minimizer
7414 ================================================
7416 The following function initializes a multidimensional minimizer. The
7417 minimizer itself depends only on the dimension of the problem and the
7418 algorithm and can be reused for different problems.
7420 -- Function: gsl_multimin_fdfminimizer *
7421 gsl_multimin_fdfminimizer_alloc (const gsl_multimin_fdfminimizer_type *
7423 -- Function: gsl_multimin_fminimizer * gsl_multimin_fminimizer_alloc
7424 (const gsl_multimin_fminimizer_type * T, size_t N)
7425 This function returns a pointer to a newly allocated instance of a
7426 minimizer of type T for an N-dimension function. If there is
7427 insufficient memory to create the minimizer then the function
7428 returns a null pointer and the error handler is invoked with an
7429 error code of `GSL_ENOMEM'.
7431 -- Function: int gsl_multimin_fdfminimizer_set
7432 (gsl_multimin_fdfminimizer * S, gsl_multimin_function_fdf *
7433 FDF, const gsl_vector * X, double STEP_SIZE, double TOL)
7434 This function initializes the minimizer S to minimize the function
7435 FDF starting from the initial point X. The size of the first
7436 trial step is given by STEP_SIZE. The accuracy of the line
7437 minimization is specified by TOL. The precise meaning of this
7438 parameter depends on the method used. Typically the line
7439 minimization is considered successful if the gradient of the
7440 function g is orthogonal to the current search direction p to a
7441 relative accuracy of TOL, where dot(p,g) < tol |p| |g|. A TOL
7442 value of 0.1 is suitable for most purposes, since line
7443 minimization only needs to be carried out approximately. Note
7444 that setting TOL to zero will force the use of "exact"
7445 line-searches, which are extremely expensive.
7447 -- Function: int gsl_multimin_fminimizer_set (gsl_multimin_fminimizer
7448 * S, gsl_multimin_function * F, const gsl_vector * X, const
7449 gsl_vector * STEP_SIZE)
7450 This function initializes the minimizer S to minimize the function
7451 F, starting from the initial point X. The size of the initial
7452 trial steps is given in vector STEP_SIZE. The precise meaning of
7453 this parameter depends on the method used.
7455 -- Function: void gsl_multimin_fdfminimizer_free
7456 (gsl_multimin_fdfminimizer * S)
7457 -- Function: void gsl_multimin_fminimizer_free
7458 (gsl_multimin_fminimizer * S)
7459 This function frees all the memory associated with the minimizer S.
7461 -- Function: const char * gsl_multimin_fdfminimizer_name (const
7462 gsl_multimin_fdfminimizer * S)
7463 -- Function: const char * gsl_multimin_fminimizer_name (const
7464 gsl_multimin_fminimizer * S)
7465 This function returns a pointer to the name of the minimizer. For
7468 printf ("s is a '%s' minimizer\n",
7469 gsl_multimin_fdfminimizer_name (s));
7471 would print something like `s is a 'conjugate_pr' minimizer'.