This is gsl-ref.info, produced by makeinfo version 4.6 from
gsl-ref.texi.

INFO-DIR-SECTION Scientific software
START-INFO-DIR-ENTRY
* gsl-ref: (gsl-ref).                   GNU Scientific Library - Reference
END-INFO-DIR-ENTRY


File: gsl-ref.info,  Node: Updating and accessing 2D histogram elements,  Next: Searching 2D histogram ranges,  Prev: Copying 2D Histograms,  Up: Histograms

Updating and accessing 2D histogram elements
============================================

You can access the bins of a two-dimensional histogram either by
specifying a pair of (x,y) coordinates or by using the bin indices
(i,j) directly.  The functions for accessing the histogram through
(x,y) coordinates use binary searches in the x and y directions to
identify the bin which covers the appropriate range.

 - Function: int gsl_histogram2d_increment (gsl_histogram2d * H, double
          X, double Y)
     This function updates the histogram H by adding one (1.0) to the
     bin whose x and y ranges contain the coordinates (X,Y).

     If the point (x,y) lies inside the valid ranges of the histogram
     then the function returns zero to indicate success.  If (x,y) lies
     outside the limits of the histogram then the function returns
     `GSL_EDOM', and none of the bins are modified.  The error handler
     is not called, since it is often necessary to compute histograms
     for a small range of a larger dataset, ignoring any coordinates
     outside the range of interest.

 - Function: int gsl_histogram2d_accumulate (gsl_histogram2d * H,
          double X, double Y, double WEIGHT)
     This function is similar to `gsl_histogram2d_increment' but
     increases the value of the appropriate bin in the histogram H by
     the floating-point number WEIGHT.

 - Function: double gsl_histogram2d_get (const gsl_histogram2d * H,
          size_t I, size_t J)
     This function returns the contents of the (I,J)th bin of the
     histogram H.  If (I,J) lies outside the valid range of indices for
     the histogram then the error handler is called with an error code
     of `GSL_EDOM' and the function returns 0.

 - Function: int gsl_histogram2d_get_xrange (const gsl_histogram2d * H,
          size_t I, double * XLOWER, double * XUPPER)
 - Function: int gsl_histogram2d_get_yrange (const gsl_histogram2d * H,
          size_t J, double * YLOWER, double * YUPPER)
     These functions find the upper and lower range limits of the I-th
     and J-th bins in the x and y directions of the histogram H.  The
     range limits are stored in XLOWER and XUPPER or YLOWER and YUPPER.
     The lower limits are inclusive (i.e. events with these
     coordinates are included in the bin) and the upper limits are
     exclusive (i.e. events with the value of the upper limit are not
     included and fall in the neighboring higher bin, if it exists).
     The functions return 0 to indicate success.  If I or J lies
     outside the valid range of indices for the histogram then the
     error handler is called with an error code of `GSL_EDOM'.

 - Function: double gsl_histogram2d_xmax (const gsl_histogram2d * H)
 - Function: double gsl_histogram2d_xmin (const gsl_histogram2d * H)
 - Function: size_t gsl_histogram2d_nx (const gsl_histogram2d * H)
 - Function: double gsl_histogram2d_ymax (const gsl_histogram2d * H)
 - Function: double gsl_histogram2d_ymin (const gsl_histogram2d * H)
 - Function: size_t gsl_histogram2d_ny (const gsl_histogram2d * H)
     These functions return the maximum upper and minimum lower range
     limits and the number of bins for the x and y directions of the
     histogram H.  They provide a way of determining these values
     without accessing the `gsl_histogram2d' struct directly.

 - Function: void gsl_histogram2d_reset (gsl_histogram2d * H)
     This function resets all the bins of the histogram H to zero.


File: gsl-ref.info,  Node: Searching 2D histogram ranges,  Next: 2D Histogram Statistics,  Prev: Updating and accessing 2D histogram elements,  Up: Histograms

Searching 2D histogram ranges
=============================

The following functions are used by the access and update routines to
locate the bin which corresponds to a given (x,y) coordinate.

 - Function: int gsl_histogram2d_find (const gsl_histogram2d * H,
          double X, double Y, size_t * I, size_t * J)
     This function finds and sets the indices I and J to the to the bin
     which covers the coordinates (X,Y). The bin is located using a
     binary search.  The search includes an optimization for histogram
     with uniform ranges, and will return the correct bin immediately
     in this case. If (x,y) is found then the function sets the indices
     (I,J) and returns `GSL_SUCCESS'.  If (x,y) lies outside the valid
     range of the histogram then the function returns `GSL_EDOM' and
     the error handler is invoked.


File: gsl-ref.info,  Node: 2D Histogram Statistics,  Next: 2D Histogram Operations,  Prev: Searching 2D histogram ranges,  Up: Histograms

2D Histogram Statistics
=======================

 - Function: double gsl_histogram2d_max_val (const gsl_histogram2d * H)
     This function returns the maximum value contained in the histogram
     bins.

 - Function: void gsl_histogram2d_max_bin (const gsl_histogram2d * H,
          size_t * I, size_t * J)
     This function returns the indices (I,J) of the bin containing the
     maximum value in the histogram H. In the case where several bins
     contain the same maximum value the first bin found is returned.

 - Function: double gsl_histogram2d_min_val (const gsl_histogram2d * H)
     This function returns the minimum value contained in the histogram
     bins.

 - Function: void gsl_histogram2d_min_bin (const gsl_histogram2d * H,
          size_t * I, size_t * J)
     This function returns the indices (I,J) of the bin containing the
     minimum value in the histogram H. In the case where several bins
     contain the same maximum value the first bin found is returned.

 - Function: double gsl_histogram2d_xmean (const gsl_histogram2d * H)
     This function returns the mean of the histogrammed x variable,
     where the histogram is regarded as a probability distribution.
     Negative bin values are ignored for the purposes of this
     calculation.

 - Function: double gsl_histogram2d_ymean (const gsl_histogram2d * H)
     This function returns the mean of the histogrammed y variable,
     where the histogram is regarded as a probability distribution.
     Negative bin values are ignored for the purposes of this
     calculation.

 - Function: double gsl_histogram2d_xsigma (const gsl_histogram2d * H)
     This function returns the standard deviation of the histogrammed x
     variable, where the histogram is regarded as a probability
     distribution. Negative bin values are ignored for the purposes of
     this calculation.

 - Function: double gsl_histogram2d_ysigma (const gsl_histogram2d * H)
     This function returns the standard deviation of the histogrammed y
     variable, where the histogram is regarded as a probability
     distribution. Negative bin values are ignored for the purposes of
     this calculation.

 - Function: double gsl_histogram2d_cov (const gsl_histogram2d * H)
     This function returns the covariance of the histogrammed x and y
     variables, where the histogram is regarded as a probability
     distribution. Negative bin values are ignored for the purposes of
     this calculation.

 - Function: double gsl_histogram2d_sum (const gsl_histogram2d * H)
     This function returns the sum of all bin values. Negative bin
     values are included in the sum.


File: gsl-ref.info,  Node: 2D Histogram Operations,  Next: Reading and writing 2D histograms,  Prev: 2D Histogram Statistics,  Up: Histograms

2D Histogram Operations
=======================

 - Function: int gsl_histogram2d_equal_bins_p (const gsl_histogram2d
          *H1, const gsl_histogram2d *H2)
     This function returns 1 if the all of the individual bin ranges of
     the two histograms are identical, and 0 otherwise.

 - Function: int gsl_histogram2d_add (gsl_histogram2d *H1, const
          gsl_histogram2d *H2)
     This function adds the contents of the bins in histogram H2 to the
     corresponding bins of histogram H1, i.e. h'_1(i,j) = h_1(i,j) +
     h_2(i,j).  The two histograms must have identical bin ranges.

 - Function: int gsl_histogram2d_sub (gsl_histogram2d *H1, const
          gsl_histogram2d *H2)
     This function subtracts the contents of the bins in histogram H2
     from the corresponding bins of histogram H1, i.e. h'_1(i,j) =
     h_1(i,j) - h_2(i,j).  The two histograms must have identical bin
     ranges.

 - Function: int gsl_histogram2d_mul (gsl_histogram2d *H1, const
          gsl_histogram2d *H2)
     This function multiplies the contents of the bins of histogram H1
     by the contents of the corresponding bins in histogram H2, i.e.
     h'_1(i,j) = h_1(i,j) * h_2(i,j).  The two histograms must have
     identical bin ranges.

 - Function: int gsl_histogram2d_div (gsl_histogram2d *H1, const
          gsl_histogram2d *H2)
     This function divides the contents of the bins of histogram H1 by
     the contents of the corresponding bins in histogram H2, i.e.
     h'_1(i,j) = h_1(i,j) / h_2(i,j).  The two histograms must have
     identical bin ranges.

 - Function: int gsl_histogram2d_scale (gsl_histogram2d *H, double
          SCALE)
     This function multiplies the contents of the bins of histogram H
     by the constant SCALE, i.e. h'_1(i,j) = h_1(i,j) scale.

 - Function: int gsl_histogram2d_shift (gsl_histogram2d *H, double
          OFFSET)
     This function shifts the contents of the bins of histogram H by
     the constant OFFSET, i.e. h'_1(i,j) = h_1(i,j) + offset.


File: gsl-ref.info,  Node: Reading and writing 2D histograms,  Next: Resampling from 2D histograms,  Prev: 2D Histogram Operations,  Up: Histograms

Reading and writing 2D histograms
=================================

The library provides functions for reading and writing two dimensional
histograms to a file as binary data or formatted text.

 - Function: int gsl_histogram2d_fwrite (FILE * STREAM, const
          gsl_histogram2d * H)
     This function writes the ranges and bins of the histogram H to the
     stream STREAM in binary format.  The return value is 0 for success
     and `GSL_EFAILED' if there was a problem writing to the file.
     Since the data is written in the native binary format it may not
     be portable between different architectures.

 - Function: int gsl_histogram2d_fread (FILE * STREAM, gsl_histogram2d
          * H)
     This function reads into the histogram H from the stream STREAM in
     binary format.  The histogram H must be preallocated with the
     correct size since the function uses the number of x and y bins in
     H to determine how many bytes to read.  The return value is 0 for
     success and `GSL_EFAILED' if there was a problem reading from the
     file.  The data is assumed to have been written in the native
     binary format on the same architecture.

 - Function: int gsl_histogram2d_fprintf (FILE * STREAM, const
          gsl_histogram2d * H, const char * RANGE_FORMAT, const char *
          BIN_FORMAT)
     This function writes the ranges and bins of the histogram H
     line-by-line to the stream STREAM using the format specifiers
     RANGE_FORMAT and BIN_FORMAT.  These should be one of the `%g',
     `%e' or `%f' formats for floating point numbers.  The function
     returns 0 for success and `GSL_EFAILED' if there was a problem
     writing to the file.  The histogram output is formatted in five
     columns, and the columns are separated by spaces, like this,

          xrange[0] xrange[1] yrange[0] yrange[1] bin(0,0)
          xrange[0] xrange[1] yrange[1] yrange[2] bin(0,1)
          xrange[0] xrange[1] yrange[2] yrange[3] bin(0,2)
          ....
          xrange[0] xrange[1] yrange[ny-1] yrange[ny] bin(0,ny-1)
          
          xrange[1] xrange[2] yrange[0] yrange[1] bin(1,0)
          xrange[1] xrange[2] yrange[1] yrange[2] bin(1,1)
          xrange[1] xrange[2] yrange[1] yrange[2] bin(1,2)
          ....
          xrange[1] xrange[2] yrange[ny-1] yrange[ny] bin(1,ny-1)
          
          ....
          
          xrange[nx-1] xrange[nx] yrange[0] yrange[1] bin(nx-1,0)
          xrange[nx-1] xrange[nx] yrange[1] yrange[2] bin(nx-1,1)
          xrange[nx-1] xrange[nx] yrange[1] yrange[2] bin(nx-1,2)
          ....
          xrange[nx-1] xrange[nx] yrange[ny-1] yrange[ny] bin(nx-1,ny-1)

     Each line contains the lower and upper limits of the bin and the
     contents of the bin.  Since the upper limits of the each bin are
     the lower limits of the neighboring bins there is duplication of
     these values but this allows the histogram to be manipulated with
     line-oriented tools.

 - Function: int gsl_histogram2d_fscanf (FILE * STREAM, gsl_histogram2d
          * H)
     This function reads formatted data from the stream STREAM into the
     histogram H.  The data is assumed to be in the five-column format
     used by `gsl_histogram_fprintf'.  The histogram H must be
     preallocated with the correct lengths since the function uses the
     sizes of H to determine how many numbers to read.  The function
     returns 0 for success and `GSL_EFAILED' if there was a problem
     reading from the file.


File: gsl-ref.info,  Node: Resampling from 2D histograms,  Next: Example programs for 2D histograms,  Prev: Reading and writing 2D histograms,  Up: Histograms

Resampling from 2D histograms
=============================

As in the one-dimensional case, a two-dimensional histogram made by
counting events can be regarded as a measurement of a probability
distribution.  Allowing for statistical error, the height of each bin
represents the probability of an event where (x,y) falls in the range
of that bin.  For a two-dimensional histogram the probability
distribution takes the form p(x,y) dx dy where,

     p(x,y) = n_{ij}/ (N A_{ij})

In this equation n_{ij} is the number of events in the bin which
contains (x,y), A_{ij} is the area of the bin and N is the total number
of events.  The distribution of events within each bin is assumed to be
uniform.

 - Data Type: gsl_histogram2d_pdf
    `size_t nx, ny'
          This is the number of histogram bins used to approximate the
          probability distribution function in the x and y directions.

    `double * xrange'
          The ranges of the bins in the x-direction are stored in an
          array of NX + 1 elements pointed to by XRANGE.

    `double * yrange'
          The ranges of the bins in the y-direction are stored in an
          array of NY + 1 pointed to by YRANGE.

    `double * sum'
          The cumulative probability for the bins is stored in an array
          of NX*NY elements pointed to by SUM.

The following functions allow you to create a `gsl_histogram2d_pdf'
struct which represents a two dimensional probability distribution and
generate random samples from it.

 - Function: gsl_histogram2d_pdf * gsl_histogram2d_pdf_alloc (size_t
          NX, size_t NY)
     This function allocates memory for a two-dimensional probability
     distribution of size NX-by-NY and returns a pointer to a newly
     initialized `gsl_histogram2d_pdf' struct. If insufficient memory
     is available a null pointer is returned and the error handler is
     invoked with an error code of `GSL_ENOMEM'.

 - Function: int gsl_histogram2d_pdf_init (gsl_histogram2d_pdf * P,
          const gsl_histogram2d * H)
     This function initializes the two-dimensional probability
     distribution calculated P from the histogram H.  If any of the
     bins of H are negative then the error handler is invoked with an
     error code of `GSL_EDOM' because a probability distribution cannot
     contain negative values.

 - Function: void gsl_histogram2d_pdf_free (gsl_histogram2d_pdf * P)
     This function frees the two-dimensional probability distribution
     function P and all of the memory associated with it.

 - Function: int gsl_histogram2d_pdf_sample (const gsl_histogram2d_pdf
          * P, double R1, double R2, double * X, double * Y)
     This function uses two uniform random numbers between zero and one,
     R1 and R2, to compute a single random sample from the
     two-dimensional probability distribution P.


File: gsl-ref.info,  Node: Example programs for 2D histograms,  Prev: Resampling from 2D histograms,  Up: Histograms

Example programs for 2D histograms
==================================

This program demonstrates two features of two-dimensional histograms.
First a 10-by-10 two-dimensional histogram is created with x and y
running from 0 to 1.  Then a few sample points are added to the
histogram, at (0.3,0.3) with a height of 1, at (0.8,0.1) with a height
of 5 and at (0.7,0.9) with a height of 0.5.  This histogram with three
events is used to generate a random sample of 1000 simulated events,
which are printed out.

     #include <stdio.h>
     #include <gsl/gsl_rng.h>
     #include <gsl/gsl_histogram2d.h>
     
     int
     main (void)
     {
       const gsl_rng_type * T;
       gsl_rng * r;
     
       gsl_histogram2d * h = gsl_histogram2d_alloc (10, 10);
     
       gsl_histogram2d_set_ranges_uniform (h, 
                                           0.0, 1.0,
                                           0.0, 1.0);
     
       gsl_histogram2d_accumulate (h, 0.3, 0.3, 1);
       gsl_histogram2d_accumulate (h, 0.8, 0.1, 5);
       gsl_histogram2d_accumulate (h, 0.7, 0.9, 0.5);
     
       gsl_rng_env_setup ();
       
       T = gsl_rng_default;
       r = gsl_rng_alloc (T);
     
       {
         int i;
         gsl_histogram2d_pdf * p 
           = gsl_histogram2d_pdf_alloc (h->nx, h->ny);
         
         gsl_histogram2d_pdf_init (p, h);
     
         for (i = 0; i < 1000; i++) {
           double x, y;
           double u = gsl_rng_uniform (r);
           double v = gsl_rng_uniform (r);
            
           gsl_histogram2d_pdf_sample (p, u, v, &x, &y);
           
           printf ("%g %g\n", x, y);
         }
       }
       
      return 0;
     }

The following plot shows the distribution of the simulated events.
Using a higher resolution grid we can see the original underlying
histogram and also the statistical fluctuations caused by the events
being uniformly distributed over the area of the original bins.


File: gsl-ref.info,  Node: N-tuples,  Next: Monte Carlo Integration,  Prev: Histograms,  Up: Top

N-tuples
********

This chapter describes functions for creating and manipulating
"ntuples", sets of values associated with events.  The ntuples are
stored in files. Their values can be extracted in any combination and
"booked" in a histogram using a selection function.

   The values to be stored are held in a user-defined data structure,
and an ntuple is created associating this data structure with a file.
The values are then written to the file (normally inside a loop) using
the ntuple functions described below.

   A histogram can be created from ntuple data by providing a selection
function and a value function.  The selection function specifies whether
an event should be included in the subset to be analyzed or not. The
value function computes the entry to be added to the histogram for each
event.

   All the ntuple functions are defined in the header file
`gsl_ntuple.h'

* Menu:

* The ntuple struct::
* Creating ntuples::
* Opening an existing ntuple file::
* Writing ntuples::
* Reading ntuples ::
* Closing an ntuple file::
* Histogramming ntuple values::
* Example ntuple programs::
* Ntuple References and Further Reading::


File: gsl-ref.info,  Node: The ntuple struct,  Next: Creating ntuples,  Up: N-tuples

The ntuple struct
=================

Ntuples are manipulated using the `gsl_ntuple' struct. This struct
contains information on the file where the ntuple data is stored, a
pointer to the current ntuple data row and the size of the user-defined
ntuple data struct.

     typedef struct {
         FILE * file;
         void * ntuple_data;
         size_t size;
     } gsl_ntuple;


File: gsl-ref.info,  Node: Creating ntuples,  Next: Opening an existing ntuple file,  Prev: The ntuple struct,  Up: N-tuples

Creating ntuples
================

 - Function: gsl_ntuple * gsl_ntuple_create (char * FILENAME, void *
          NTUPLE_DATA, size_t SIZE)
     This function creates a new write-only ntuple file FILENAME for
     ntuples of size SIZE and returns a pointer to the newly created
     ntuple struct.  Any existing file with the same name is truncated
     to zero length and overwritten.  A pointer to memory for the
     current ntuple row NTUPLE_DATA must be supplied--this is used to
     copy ntuples in and out of the file.


File: gsl-ref.info,  Node: Opening an existing ntuple file,  Next: Writing ntuples,  Prev: Creating ntuples,  Up: N-tuples

Opening an existing ntuple file
===============================

 - Function: gsl_ntuple * gsl_ntuple_open (char * FILENAME, void *
          NTUPLE_DATA, size_t SIZE)
     This function opens an existing ntuple file FILENAME for reading
     and returns a pointer to a corresponding ntuple struct. The
     ntuples in the file must have size SIZE.  A pointer to memory for
     the current ntuple row NTUPLE_DATA must be supplied--this is used
     to copy ntuples in and out of the file.


File: gsl-ref.info,  Node: Writing ntuples,  Next: Reading ntuples,  Prev: Opening an existing ntuple file,  Up: N-tuples

Writing ntuples
===============

 - Function: int gsl_ntuple_write (gsl_ntuple * NTUPLE)
     This function writes the current ntuple NTUPLE->NTUPLE_DATA of
     size NTUPLE->SIZE to the corresponding file.

 - Function: int gsl_ntuple_bookdata (gsl_ntuple * NTUPLE)
     This function is a synonym for `gsl_ntuple_write'.


File: gsl-ref.info,  Node: Reading ntuples,  Next: Closing an ntuple file,  Prev: Writing ntuples,  Up: N-tuples

Reading ntuples
===============

 - Function: int gsl_ntuple_read (gsl_ntuple * NTUPLE)
     This function reads the current row of the ntuple file for NTUPLE
     and stores the values in NTUPLE->DATA.


File: gsl-ref.info,  Node: Closing an ntuple file,  Next: Histogramming ntuple values,  Prev: Reading ntuples,  Up: N-tuples

Closing an ntuple file
======================

 - Function: int gsl_ntuple_close (gsl_ntuple * NTUPLE)
     This function closes the ntuple file NTUPLE and frees its
     associated allocated memory.


File: gsl-ref.info,  Node: Histogramming ntuple values,  Next: Example ntuple programs,  Prev: Closing an ntuple file,  Up: N-tuples

Histogramming ntuple values
===========================

Once an ntuple has been created its contents can be histogrammed in
various ways using the function `gsl_ntuple_project'.  Two user-defined
functions must be provided, a function to select events and a function
to compute scalar values. The selection function and the value function
both accept the ntuple row as a first argument and other parameters as
a second argument.

   The "selection function" determines which ntuple rows are selected
for histogramming.  It is defined by the following struct,
     typedef struct {
       int (* function) (void * ntuple_data, void * params);
       void * params;
     } gsl_ntuple_select_fn;

The struct component FUNCTION should return a non-zero value for each
ntuple row that is to be included in the histogram.

   The "value function" computes scalar values for those ntuple rows
selected by the selection function,
     typedef struct {
       double (* function) (void * ntuple_data, void * params);
       void * params;
     } gsl_ntuple_value_fn;

In this case the struct component FUNCTION should return the value to
be added to the histogram for the ntuple row.

 - Function: int gsl_ntuple_project (gsl_histogram * H, gsl_ntuple *
          NTUPLE, gsl_ntuple_value_fn *VALUE_FUNC, gsl_ntuple_select_fn
          *SELECT_FUNC)
     This function updates the histogram H from the ntuple NTUPLE using
     the functions VALUE_FUNC and SELECT_FUNC. For each ntuple row
     where the selection function SELECT_FUNC is non-zero the
     corresponding value of that row is computed using the function
     VALUE_FUNC and added to the histogram.  Those ntuple rows where
     SELECT_FUNC returns zero are ignored.  New entries are added to
     the histogram, so subsequent calls can be used to accumulate
     further data in the same histogram.


File: gsl-ref.info,  Node: Example ntuple programs,  Next: Ntuple References and Further Reading,  Prev: Histogramming ntuple values,  Up: N-tuples

Examples
========

The following example programs demonstrate the use of ntuples in
managing a large dataset.  The first program creates a set of 10,000
simulated "events", each with 3 associated values (x,y,z).  These are
generated from a gaussian distribution with unit variance, for
demonstration purposes, and written to the ntuple file `test.dat'.

     #include <gsl/gsl_ntuple.h>
     #include <gsl/gsl_rng.h>
     #include <gsl/gsl_randist.h>
     
     struct data
     {
       double x;
       double y;
       double z;
     };
     
     int
     main (void)
     {
       const gsl_rng_type * T;
       gsl_rng * r;
     
       struct data ntuple_row;
       int i;
     
       gsl_ntuple *ntuple 
         = gsl_ntuple_create ("test.dat", &ntuple_row, 
                              sizeof (ntuple_row));
     
       gsl_rng_env_setup ();
     
       T = gsl_rng_default; 
       r = gsl_rng_alloc (T);
     
       for (i = 0; i < 10000; i++)
         {
           ntuple_row.x = gsl_ran_ugaussian (r);
           ntuple_row.y = gsl_ran_ugaussian (r);
           ntuple_row.z = gsl_ran_ugaussian (r);
           
           gsl_ntuple_write (ntuple);
         }
       
       gsl_ntuple_close (ntuple);
       return 0;
     }

The next program analyses the ntuple data in the file `test.dat'.  The
analysis procedure is to compute the squared-magnitude of each event,
E^2=x^2+y^2+z^2, and select only those which exceed a lower limit of
1.5.  The selected events are then histogrammed using their E^2 values.

     #include <math.h>
     #include <gsl/gsl_ntuple.h>
     #include <gsl/gsl_histogram.h>
     
     struct data
     {
       double x;
       double y;
       double z;
     };
     
     int sel_func (void *ntuple_data, void *params);
     double val_func (void *ntuple_data, void *params);
     
     int
     main (void)
     {
       struct data ntuple_row;
     
       gsl_ntuple *ntuple 
         = gsl_ntuple_open ("test.dat", &ntuple_row,
                            sizeof (ntuple_row));
       double lower = 1.5;
     
       gsl_ntuple_select_fn S;
       gsl_ntuple_value_fn V;
     
       gsl_histogram *h = gsl_histogram_alloc (100);
       gsl_histogram_set_ranges_uniform(h, 0.0, 10.0);
     
       S.function = &sel_func;
       S.params = &lower;
     
       V.function = &val_func;
       V.params = 0;
     
       gsl_ntuple_project (h, ntuple, &V, &S);
       gsl_histogram_fprintf (stdout, h, "%f", "%f");
       gsl_histogram_free (h);
       gsl_ntuple_close (ntuple);
     
       return 0;
     }
     
     int
     sel_func (void *ntuple_data, void *params)
     {
       struct data * data = (struct data *) ntuple_data;  
       double x, y, z, E2, scale;
       scale = *(double *) params;
       
       x = data->x;
       y = data->y;
       z = data->z;
     
       E2 = x * x + y * y + z * z;
     
       return E2 > scale;
     }
     
     double
     val_func (void *ntuple_data, void *params)
     {
       struct data * data = (struct data *) ntuple_data;  
       double x, y, z;
     
       x = data->x;
       y = data->y;
       z = data->z;
     
       return x * x + y * y + z * z;
     }

   The following plot shows the distribution of the selected events.
Note the cut-off at the lower bound.


File: gsl-ref.info,  Node: Ntuple References and Further Reading,  Prev: Example ntuple programs,  Up: N-tuples

References and Further Reading
==============================

Further information on the use of ntuples can be found in the
documentation for the CERN packages PAW and HBOOK (available online).


File: gsl-ref.info,  Node: Monte Carlo Integration,  Next: Simulated Annealing,  Prev: N-tuples,  Up: Top

Monte Carlo Integration
***********************

This chapter describes routines for multidimensional Monte Carlo
integration.  These include the traditional Monte Carlo method and
adaptive algorithms such as VEGAS and MISER which use importance
sampling and stratified sampling techniques. Each algorithm computes an
estimate of a multidimensional definite integral of the form,

     I = \int_xl^xu dx \int_yl^yu  dy ...  f(x, y, ...)

over a hypercubic region ((x_l,x_u), (y_l,y_u), ...) using a fixed
number of function calls.  The routines also provide a statistical
estimate of the error on the result.  This error estimate should be
taken as a guide rather than as a strict error bound -- random sampling
of the region may not uncover all the important features of the
function, resulting in an underestimate of the error.

   The functions are defined in separate header files for each routine,
`gsl_monte_plain.h', `gsl_monte_miser.h' and `gsl_monte_vegas.h'.

* Menu:

* Monte Carlo Interface::
* PLAIN Monte Carlo::
* MISER::
* VEGAS::
* Monte Carlo Examples::
* Monte Carlo Integration References and Further Reading::


File: gsl-ref.info,  Node: Monte Carlo Interface,  Next: PLAIN Monte Carlo,  Up: Monte Carlo Integration

Interface
=========

All of the Monte Carlo integration routines use the same general form of
interface.  There is an allocator to allocate memory for control
variables and workspace, a routine to initialize those control
variables, the integrator itself, and a function to free the space when
done.

   Each integration function requires a random number generator to be
supplied, and returns an estimate of the integral and its standard
deviation.  The accuracy of the result is determined by the number of
function calls specified by the user.  If a known level of accuracy is
required this can be achieved by calling the integrator several times
and averaging the individual results until the desired accuracy is
obtained.

   Random sample points used within the Monte Carlo routines are always
chosen strictly within the integration region, so that endpoint
singularities are automatically avoided.

   The function to be integrated has its own datatype, defined in the
header file `gsl_monte.h'.

 - Data Type: gsl_monte_function
     This data type defines a general function with parameters for Monte
     Carlo integration.

    `double (* f) (double * X, size_t DIM, void * PARAMS)'
          this function should return the value f(x,params) for
          argument X and parameters PARAMS, where X is an array of size
          DIM giving the coordinates of the point where the function is
          to be evaluated.

    `size_t dim'
          the number of dimensions for X.

    `void * params'
          a pointer to the parameters of the function.

Here is an example for a quadratic function in two dimensions,

     f(x,y) = a x^2 + b x y + c y^2

with a = 3, b = 2, c = 1.  The following code defines a
`gsl_monte_function' `F' which you could pass to an integrator:

     struct my_f_params { double a; double b; double c; };
     
     double
     my_f (double x[], size_t dim, void * p) {
        struct my_f_params * fp = (struct my_f_params *)p;
     
        if (dim != 2)
           {
             fprintf (stderr, "error: dim != 2");
             abort ();
           }
     
        return  fp->a * x[0] * x[0]
                  + fp->b * x[0] * x[1]
                    + fp->c * x[1] * x[1];
     }
     
     gsl_monte_function F;
     struct my_f_params params = { 3.0, 2.0, 1.0 };
     
     F.f = &my_f;
     F.dim = 2;
     F.params = &params;

The function f(x) can be evaluated using the following macro,

     #define GSL_MONTE_FN_EVAL(F,x)
         (*((F)->f))(x,(F)->dim,(F)->params)


File: gsl-ref.info,  Node: PLAIN Monte Carlo,  Next: MISER,  Prev: Monte Carlo Interface,  Up: Monte Carlo Integration

PLAIN Monte Carlo
=================

The plain Monte Carlo algorithm samples points randomly from the
integration region to estimate the integral and its error.  Using this
algorithm the estimate of the integral E(f; N) for N randomly
distributed points x_i is given by,

     E(f; N) = =  V <f> = (V / N) \sum_i^N f(x_i)

where V is the volume of the integration region.  The error on this
estimate \sigma(E;N) is calculated from the estimated variance of the
mean,

     \sigma^2 (E; N) = (V / N) \sum_i^N (f(x_i) -  <f>)^2.

For large N this variance decreases asymptotically as var(f)/N, where
var(f) is the true variance of the function over the integration
region.  The error estimate itself should decrease as
\sigma(f)/\sqrt{N}.  The familiar law of errors decreasing as
1/\sqrt{N} applies -- to reduce the error by a factor of 10 requires a
100-fold increase in the number of sample points.

   The functions described in this section are declared in the header
file `gsl_monte_plain.h'.

 - Function: gsl_monte_plain_state * gsl_monte_plain_alloc (size_t DIM)
     This function allocates and initializes a workspace for Monte Carlo
     integration in DIM dimensions.

 - Function: int gsl_monte_plain_init (gsl_monte_plain_state* S)
     This function initializes a previously allocated integration state.
     This allows an existing workspace to be reused for different
     integrations.

 - Function: int gsl_monte_plain_integrate (gsl_monte_function * F,
          double * XL, double * XU, size_t DIM, size_t CALLS, gsl_rng *
          R, gsl_monte_plain_state * S, double * RESULT, double *
          ABSERR)
     This routines uses the plain Monte Carlo algorithm to integrate the
     function F over the DIM-dimensional hypercubic region defined by
     the lower and upper limits in the arrays XL and XU, each of size
     DIM.  The integration uses a fixed number of function calls CALLS,
     and obtains random sampling points using the random number
     generator R. A previously allocated workspace S must be supplied.
     The result of the integration is returned in RESULT, with an
     estimated absolute error ABSERR.

 - Function: void gsl_monte_plain_free (gsl_monte_plain_state * S)
     This function frees the memory associated with the integrator state
     S.


File: gsl-ref.info,  Node: MISER,  Next: VEGAS,  Prev: PLAIN Monte Carlo,  Up: Monte Carlo Integration

MISER
=====

The MISER algorithm of Press and Farrar is based on recursive
stratified sampling.  This technique aims to reduce the overall
integration error by concentrating integration points in the regions of
highest variance.

   The idea of stratified sampling begins with the observation that for
two disjoint regions a and b with Monte Carlo estimates of the integral
E_a(f) and E_b(f) and variances \sigma_a^2(f) and \sigma_b^2(f), the
variance Var(f) of the combined estimate E(f) = (1/2) (E_a(f) + E_b(f))
is given by,

     Var(f) = (\sigma_a^2(f) / 4 N_a) + (\sigma_b^2(f) / 4 N_b).

It can be shown that this variance is minimized by distributing the
points such that,

     N_a / (N_a + N_b) = \sigma_a / (\sigma_a + \sigma_b).

Hence the smallest error estimate is obtained by allocating sample
points in proportion to the standard deviation of the function in each
sub-region.

   The MISER algorithm proceeds by bisecting the integration region
along one coordinate axis to give two sub-regions at each step.  The
direction is chosen by examining all d possible bisections and
selecting the one which will minimize the combined variance of the two
sub-regions.  The variance in the sub-regions is estimated by sampling
with a fraction of the total number of points available to the current
step.  The same procedure is then repeated recursively for each of the
two half-spaces from the best bisection. The remaining sample points are
allocated to the sub-regions using the formula for N_a and N_b.  This
recursive allocation of integration points continues down to a
user-specified depth where each sub-region is integrated using a plain
Monte Carlo estimate.  These individual values and their error
estimates are then combined upwards to give an overall result and an
estimate of its error.

   The functions described in this section are declared in the header
file `gsl_monte_miser.h'.

 - Function: gsl_monte_miser_state * gsl_monte_miser_alloc (size_t DIM)
     This function allocates and initializes a workspace for Monte Carlo
     integration in DIM dimensions.  The workspace is used to maintain
     the state of the integration.

 - Function: int gsl_monte_miser_init (gsl_monte_miser_state* S)
     This function initializes a previously allocated integration state.
     This allows an existing workspace to be reused for different
     integrations.

 - Function: int gsl_monte_miser_integrate (gsl_monte_function * F,
          double * XL, double * XU, size_t DIM, size_t CALLS, gsl_rng *
          R, gsl_monte_miser_state * S, double * RESULT, double *
          ABSERR)
     This routines uses the MISER Monte Carlo algorithm to integrate the
     function F over the DIM-dimensional hypercubic region defined by
     the lower and upper limits in the arrays XL and XU, each of size
     DIM.  The integration uses a fixed number of function calls CALLS,
     and obtains random sampling points using the random number
     generator R. A previously allocated workspace S must be supplied.
     The result of the integration is returned in RESULT, with an
     estimated absolute error ABSERR.

 - Function: void gsl_monte_miser_free (gsl_monte_miser_state * S)
     This function frees the memory associated with the integrator state
     S.

   The MISER algorithm has several configurable parameters. The
following variables can be accessed through the `gsl_monte_miser_state'
struct,

 - Variable: double estimate_frac
     This parameter specifies the fraction of the currently available
     number of function calls which are allocated to estimating the
     variance at each recursive step. The default value is 0.1.

 - Variable: size_t min_calls
     This parameter specifies the minimum number of function calls
     required for each estimate of the variance. If the number of
     function calls allocated to the estimate using ESTIMATE_FRAC falls
     below MIN_CALLS then MIN_CALLS are used instead.  This ensures
     that each estimate maintains a reasonable level of accuracy.  The
     default value of MIN_CALLS is `16 * dim'.

 - Variable: size_t min_calls_per_bisection
     This parameter specifies the minimum number of function calls
     required to proceed with a bisection step.  When a recursive step
     has fewer calls available than MIN_CALLS_PER_BISECTION it performs
     a plain Monte Carlo estimate of the current sub-region and
     terminates its branch of the recursion.  The default value of this
     parameter is `32 * min_calls'.

 - Variable: double alpha
     This parameter controls how the estimated variances for the two
     sub-regions of a bisection are combined when allocating points.
     With recursive sampling the overall variance should scale better
     than 1/N, since the values from the sub-regions will be obtained
     using a procedure which explicitly minimizes their variance.  To
     accommodate this behavior the MISER algorithm allows the total
     variance to depend on a scaling parameter \alpha,

          Var(f) = {\sigma_a \over N_a^\alpha} + {\sigma_b \over N_b^\alpha}.

     The authors of the original paper describing MISER recommend the
     value \alpha = 2 as a good choice, obtained from numerical
     experiments, and this is used as the default value in this
     implementation.

 - Variable: double dither
     This parameter introduces a random fractional variation of size
     DITHER into each bisection, which can be used to break the
     symmetry of integrands which are concentrated near the exact
     center of the hypercubic integration region.  The default value of
     dither is zero, so no variation is introduced. If needed, a
     typical value of DITHER is around 0.1.


File: gsl-ref.info,  Node: VEGAS,  Next: Monte Carlo Examples,  Prev: MISER,  Up: Monte Carlo Integration

VEGAS
=====

The VEGAS algorithm of Lepage is based on importance sampling.  It
samples points from the probability distribution described by the
function |f|, so that the points are concentrated in the regions that
make the largest contribution to the integral.

   In general, if the Monte Carlo integral of f is sampled with points
distributed according to a probability distribution described by the
function g, we obtain an estimate E_g(f; N),

     E_g(f; N) = E(f/g; N)

with a corresponding variance,

     Var_g(f; N) = Var(f/g; N).

If the probability distribution is chosen as g = |f|/I(|f|) then it can
be shown that the variance V_g(f; N) vanishes, and the error in the
estimate will be zero.  In practice it is not possible to sample from
the exact distribution g for an arbitrary function, so importance
sampling algorithms aim to produce efficient approximations to the
desired distribution.

   The VEGAS algorithm approximates the exact distribution by making a
number of passes over the integration region while histogramming the
function f. Each histogram is used to define a sampling distribution
for the next pass.  Asymptotically this procedure converges to the
desired distribution. In order to avoid the number of histogram bins
growing like K^d the probability distribution is approximated by a
separable function: g(x_1, x_2, ...) = g_1(x_1) g_2(x_2) ...  so that
the number of bins required is only Kd.  This is equivalent to locating
the peaks of the function from the projections of the integrand onto
the coordinate axes.  The efficiency of VEGAS depends on the validity
of this assumption.  It is most efficient when the peaks of the
integrand are well-localized.  If an integrand can be rewritten in a
form which is approximately separable this will increase the efficiency
of integration with VEGAS.

   VEGAS incorporates a number of additional features, and combines both
stratified sampling and importance sampling.  The integration region is
divided into a number of "boxes", with each box getting a fixed number
of points (the goal is 2).  Each box can then have a fractional number
of bins, but if bins/box is less than two, Vegas switches to a kind
variance reduction (rather than importance sampling).

 - Function: gsl_monte_vegas_state * gsl_monte_vegas_alloc (size_t DIM)
     This function allocates and initializes a workspace for Monte Carlo
     integration in DIM dimensions.  The workspace is used to maintain
     the state of the integration.

 - Function: int gsl_monte_vegas_init (gsl_monte_vegas_state* S)
     This function initializes a previously allocated integration state.
     This allows an existing workspace to be reused for different
     integrations.

 - Function: int gsl_monte_vegas_integrate (gsl_monte_function * F,
          double * XL, double * XU, size_t DIM, size_t CALLS, gsl_rng *
          R, gsl_monte_vegas_state * S, double * RESULT, double *
          ABSERR)
     This routines uses the VEGAS Monte Carlo algorithm to integrate the
     function F over the DIM-dimensional hypercubic region defined by
     the lower and upper limits in the arrays XL and XU, each of size
     DIM.  The integration uses a fixed number of function calls CALLS,
     and obtains random sampling points using the random number
     generator R. A previously allocated workspace S must be supplied.
     The result of the integration is returned in RESULT, with an
     estimated absolute error ABSERR.  The result and its error
     estimate are based on a weighted average of independent samples.
     The chi-squared per degree of freedom for the weighted average is
     returned via the state struct component, S->CHISQ, and must be
     consistent with 1 for the weighted average to be reliable.

 - Function: void gsl_monte_vegas_free (gsl_monte_vegas_state * S)
     This function frees the memory associated with the integrator state
     S.

   The VEGAS algorithm computes a number of independent estimates of the
integral internally, according to the `iterations' parameter described
below, and returns their weighted average.  Random sampling of the
integrand can occasionally produce an estimate where the error is zero,
particularly if the function is constant in some regions. An estimate
with zero error causes the weighted average to break down and must be
handled separately. In the original Fortran implementations of VEGAS
the error estimate is made non-zero by substituting a small value
(typically `1e-30').  The implementation in GSL differs from this and
avoids the use of an arbitrary constant--it either assigns the value a
weight which is the average weight of the preceding estimates or
discards it according to the following procedure,

current estimate has zero error, weighted average has finite error
     The current estimate is assigned a weight which is the average
     weight of the preceding estimates.

current estimate has finite error, previous estimates had zero error
     The previous estimates are discarded and the weighted averaging
     procedure begins with the current estimate.

current estimate has zero error, previous estimates had zero error
     The estimates are averaged using the arithmetic mean, but no error
     is computed.

   The VEGAS algorithm is highly configurable. The following variables
can be accessed through the `gsl_monte_vegas_state' struct,

 - Variable: double result
 - Variable: double sigma
     These parameters contain the raw value of the integral RESULT and
     its error SIGMA from the last iteration of the algorithm.

 - Variable: double chisq
     This parameter gives the chi-squared per degree of freedom for the
     weighted estimate of the integral.  The value of CHISQ should be
     close to 1.  A value of CHISQ which differs significantly from 1
     indicates that the values from different iterations are
     inconsistent.  In this case the weighted error will be
     under-estimated, and further iterations of the algorithm are
     needed to obtain reliable results.

 - Variable: double alpha
     The parameter `alpha' controls the stiffness of the rebinning
     algorithm.  It is typically set between one and two. A value of
     zero prevents rebinning of the grid.  The default value is 1.5.

 - Variable: size_t iterations
     The number of iterations to perform for each call to the routine.
     The default value is 5 iterations.

 - Variable: int stage
     Setting this determines the "stage" of the calculation.  Normally,
     `stage = 0' which begins with a new uniform grid and empty weighted
     average.  Calling vegas with `stage = 1' retains the grid from the
     previous run but discards the weighted average, so that one can
     "tune" the grid using a relatively small number of points and then
     do a large run with `stage = 1' on the optimized grid.  Setting
     `stage = 2' keeps the grid and the weighted average from the
     previous run, but may increase (or decrease) the number of
     histogram bins in the grid depending on the number of calls
     available.  Choosing `stage = 3' enters at the main loop, so that
     nothing is changed, and is equivalent to performing additional
     iterations in a previous call.

 - Variable: int mode
     The possible choices are `GSL_VEGAS_MODE_IMPORTANCE',
     `GSL_VEGAS_MODE_STRATIFIED', `GSL_VEGAS_MODE_IMPORTANCE_ONLY'.
     This determines whether VEGAS will use importance sampling or
     stratified sampling, or whether it can pick on its own.  In low
     dimensions VEGAS uses strict stratified sampling (more precisely,
     stratified sampling is chosen if there are fewer than 2 bins per
     box).

 - Variable: int verbose
 - Variable: FILE * ostream
     These parameters set the level of information printed by VEGAS. All
     information is written to the stream OSTREAM.  The default setting
     of VERBOSE is `-1', which turns off all output.  A VERBOSE value
     of `0' prints summary information about the weighted average and
     final result, while a value of `1' also displays the grid
     coordinates.  A value of `2' prints information from the rebinning
     procedure for each iteration.


File: gsl-ref.info,  Node: Monte Carlo Examples,  Next: Monte Carlo Integration References and Further Reading,  Prev: VEGAS,  Up: Monte Carlo Integration

Examples
========

The example program below uses the Monte Carlo routines to estimate the
value of the following 3-dimensional integral from the theory of random
walks,

     I = \int_{-pi}^{+pi} {dk_x/(2 pi)}
         \int_{-pi}^{+pi} {dk_y/(2 pi)}
         \int_{-pi}^{+pi} {dk_z/(2 pi)}
          1 / (1 - cos(k_x)cos(k_y)cos(k_z)).

The analytic value of this integral can be shown to be I =
\Gamma(1/4)^4/(4 \pi^3) = 1.393203929685676859....  The integral gives
the mean time spent at the origin by a random walk on a body-centered
cubic lattice in three dimensions.

   For simplicity we will compute the integral over the region (0,0,0)
to (\pi,\pi,\pi) and multiply by 8 to obtain the full result.  The
integral is slowly varying in the middle of the region but has
integrable singularities at the corners (0,0,0), (0,\pi,\pi),
(\pi,0,\pi) and (\pi,\pi,0).  The Monte Carlo routines only select
points which are strictly within the integration region and so no
special measures are needed to avoid these singularities.

     #include <stdlib.h>
     #include <gsl/gsl_math.h>
     #include <gsl/gsl_monte.h>
     #include <gsl/gsl_monte_plain.h>
     #include <gsl/gsl_monte_miser.h>
     #include <gsl/gsl_monte_vegas.h>
     
     /* Computation of the integral,
     
           I = int (dx dy dz)/(2pi)^3  1/(1-cos(x)cos(y)cos(z))
     
        over (-pi,-pi,-pi) to (+pi, +pi, +pi).  The exact answer
        is Gamma(1/4)^4/(4 pi^3).  This example is taken from
        C.Itzykson, J.M.Drouffe, "Statistical Field Theory -
        Volume 1", Section 1.1, p21, which cites the original
        paper M.L.Glasser, I.J.Zucker, Proc.Natl.Acad.Sci.USA 74
        1800 (1977) */
     
     /* For simplicity we compute the integral over the region 
        (0,0,0) -> (pi,pi,pi) and multiply by 8 */
     
     double exact = 1.3932039296856768591842462603255;
     
     double
     g (double *k, size_t dim, void *params)
     {
       double A = 1.0 / (M_PI * M_PI * M_PI);
       return A / (1.0 - cos (k[0]) * cos (k[1]) * cos (k[2]));
     }
     
     void
     display_results (char *title, double result, double error)
     {
       printf ("%s ==================\n", title);
       printf ("result = % .6f\n", result);
       printf ("sigma  = % .6f\n", error);
       printf ("exact  = % .6f\n", exact);
       printf ("error  = % .6f = %.1g sigma\n", result - exact,
               fabs (result - exact) / error);
     }
     
     int
     main (void)
     {
       double res, err;
     
       double xl[3] = { 0, 0, 0 };
       double xu[3] = { M_PI, M_PI, M_PI };
     
       const gsl_rng_type *T;
       gsl_rng *r;
     
       gsl_monte_function G = { &g, 3, 0 };
     
       size_t calls = 500000;
     
       gsl_rng_env_setup ();
     
       T = gsl_rng_default;
       r = gsl_rng_alloc (T);
     
       {
         gsl_monte_plain_state *s = gsl_monte_plain_alloc (3);
         gsl_monte_plain_integrate (&G, xl, xu, 3, calls, r, s, 
                                    &res, &err);
         gsl_monte_plain_free (s);
     
         display_results ("plain", res, err);
       }
     
       {
         gsl_monte_miser_state *s = gsl_monte_miser_alloc (3);
         gsl_monte_miser_integrate (&G, xl, xu, 3, calls, r, s,
                                    &res, &err);
         gsl_monte_miser_free (s);
     
         display_results ("miser", res, err);
       }
     
       {
         gsl_monte_vegas_state *s = gsl_monte_vegas_alloc (3);
     
         gsl_monte_vegas_integrate (&G, xl, xu, 3, 10000, r, s,
                                    &res, &err);
         display_results ("vegas warm-up", res, err);
     
         printf ("converging...\n");
     
         do
           {
             gsl_monte_vegas_integrate (&G, xl, xu, 3, calls/5, r, s,
                                        &res, &err);
             printf ("result = % .6f sigma = % .6f "
                     "chisq/dof = %.1f\n", res, err, s->chisq);
           }
         while (fabs (s->chisq - 1.0) > 0.5);
     
         display_results ("vegas final", res, err);
     
         gsl_monte_vegas_free (s);
       }
       return 0;
     }

With 500,000 function calls the plain Monte Carlo algorithm achieves a
fractional error of 0.6%.  The estimated error `sigma' is consistent
with the actual error, and the computed result differs from the true
result by about one standard deviation,

     plain ==================
     result =  1.385867
     sigma  =  0.007938
     exact  =  1.393204
     error  = -0.007337 = 0.9 sigma

The MISER algorithm reduces the error by a factor of two, and also
correctly estimates the error,

     miser ==================
     result =  1.390656
     sigma  =  0.003743
     exact  =  1.393204
     error  = -0.002548 = 0.7 sigma

In the case of the VEGAS algorithm the program uses an initial warm-up
run of 10,000 function calls to prepare, or "warm up", the grid.  This
is followed by a main run with five iterations of 100,000 function
calls. The chi-squared per degree of freedom for the five iterations are
checked for consistency with 1, and the run is repeated if the results
have not converged. In this case the estimates are consistent on the
first pass.

     vegas warm-up ==================
     result =  1.386925
     sigma  =  0.002651
     exact  =  1.393204
     error  = -0.006278 = 2 sigma
     converging...
     result =  1.392957 sigma =  0.000452 chisq/dof = 1.1
     vegas final ==================
     result =  1.392957
     sigma  =  0.000452
     exact  =  1.393204
     error  = -0.000247 = 0.5 sigma

If the value of `chisq' had differed significantly from 1 it would
indicate inconsistent results, with a correspondingly underestimated
error.  The final estimate from VEGAS (using a similar number of
function calls) is significantly more accurate than the other two
algorithms.


File: gsl-ref.info,  Node: Monte Carlo Integration References and Further Reading,  Prev: Monte Carlo Examples,  Up: Monte Carlo Integration

References and Further Reading
==============================

The MISER algorithm is described in the following article,

     W.H. Press, G.R. Farrar, `Recursive Stratified Sampling for
     Multidimensional Monte Carlo Integration', Computers in Physics,
     v4 (1990), pp190-195.

The VEGAS algorithm is described in the following papers,

     G.P. Lepage, `A New Algorithm for Adaptive Multidimensional
     Integration', Journal of Computational Physics 27, 192-203, (1978)

     G.P. Lepage, `VEGAS: An Adaptive Multi-dimensional Integration
     Program', Cornell preprint CLNS 80-447, March 1980


File: gsl-ref.info,  Node: Simulated Annealing,  Next: Ordinary Differential Equations,  Prev: Monte Carlo Integration,  Up: Top

Simulated Annealing
*******************

Stochastic search techniques are used when the structure of a space is
not well understood or is not smooth, so that techniques like Newton's
method (which requires calculating Jacobian derivative matrices) cannot
be used. In particular, these techniques are frequently used to solve
combinatorial optimization problems, such as the traveling salesman
problem.

   The goal is to find a point in the space at which a real valued
"energy function" (or "cost function") is minimized.  Simulated
annealing is a minimization technique which has given good results in
avoiding local minima; it is based on the idea of taking a random walk
through the space at successively lower temperatures, where the
probability of taking a step is given by a Boltzmann distribution.

   The functions described in this chapter are declared in the header
file `gsl_siman.h'.

* Menu:

* Simulated Annealing algorithm::
* Simulated Annealing functions::
* Examples with Simulated Annealing::


File: gsl-ref.info,  Node: Simulated Annealing algorithm,  Next: Simulated Annealing functions,  Prev: Simulated Annealing,  Up: Simulated Annealing

Simulated Annealing algorithm
=============================

The simulated annealing algorithm takes random walks through the problem
space, looking for points with low energies; in these random walks, the
probability of taking a step is determined by the Boltzmann
distribution,

     p = e^{-(E_{i+1} - E_i)/(kT)}

if E_{i+1} > E_i, and p = 1 when E_{i+1} <= E_i.

   In other words, a step will occur if the new energy is lower.  If
the new energy is higher, the transition can still occur, and its
likelihood is proportional to the temperature T and inversely
proportional to the energy difference E_{i+1} - E_i.

   The temperature T is initially set to a high value, and a random
walk is carried out at that temperature.  Then the temperature is
lowered very slightly according to a "cooling schedule", for example: T
-> T/mu_T where \mu_T is slightly greater than 1.

   The slight probability of taking a step that gives higher energy is
what allows simulated annealing to frequently get out of local minima.


File: gsl-ref.info,  Node: Simulated Annealing functions,  Next: Examples with Simulated Annealing,  Prev: Simulated Annealing algorithm,  Up: Simulated Annealing

Simulated Annealing functions
=============================

 - Function: void gsl_siman_solve (const gsl_rng * R, void *X0_P,
          gsl_siman_Efunc_t EF, gsl_siman_step_t TAKE_STEP,
          gsl_siman_metric_t DISTANCE, gsl_siman_print_t
          PRINT_POSITION, gsl_siman_copy_t COPYFUNC,
          gsl_siman_copy_construct_t COPY_CONSTRUCTOR,
          gsl_siman_destroy_t DESTRUCTOR, size_t ELEMENT_SIZE,
          gsl_siman_params_t PARAMS)
     This function performs a simulated annealing search through a given
     space.  The space is specified by providing the functions EF and
     DISTANCE.  The simulated annealing steps are generated using the
     random number generator R and the function TAKE_STEP.

     The starting configuration of the system should be given by X0_P.
     The routine offers two modes for updating configurations, a
     fixed-size mode and a variable-size mode.  In the fixed-size mode
     the configuration is stored as a single block of memory of size
     ELEMENT_SIZE.  Copies of this configuration are created, copied
     and destroyed internally using the standard library functions
     `malloc', `memcpy' and `free'.  The function pointers COPYFUNC,
     COPY_CONSTRUCTOR and DESTRUCTOR should be null pointers in
     fixed-size mode.  In the variable-size mode the functions
     COPYFUNC, COPY_CONSTRUCTOR and DESTRUCTOR are used to create, copy
     and destroy configurations internally.  The variable ELEMENT_SIZE
     should be zero in the variable-size mode.

     The PARAMS structure (described below) controls the run by
     providing the temperature schedule and other tunable parameters to
     the algorithm.

     On exit the best result achieved during the search is placed in
     `*X0_P'.  If the annealing process has been successful this should
     be a good approximation to the optimal point in the space.

     If the function pointer PRINT_POSITION is not null, a debugging
     log will be printed to `stdout' with the following columns:

          number_of_iterations temperature x x-(*x0_p) Ef(x)

     and the output of the function PRINT_POSITION itself.  If
     PRINT_POSITION is null then no information is printed.

The simulated annealing routines require several user-specified
functions to define the configuration space and energy function.  The
prototypes for these functions are given below.

 - Data Type: gsl_siman_Efunc_t
     This function type should return the energy of a configuration XP.
          double (*gsl_siman_Efunc_t) (void *xp)

 - Data Type: gsl_siman_step_t
     This function type should modify the configuration XP using a
     random step taken from the generator R, up to a maximum distance of
     STEP_SIZE.
          void (*gsl_siman_step_t) (const gsl_rng *r, void *xp,
                                    double step_size)

 - Data Type: gsl_siman_metric_t
     This function type should return the distance between two
     configurations XP and YP.
          double (*gsl_siman_metric_t) (void *xp, void *yp)

 - Data Type: gsl_siman_print_t
     This function type should print the contents of the configuration
     XP.
          void (*gsl_siman_print_t) (void *xp)

 - Data Type: gsl_siman_copy_t
     This function type should copy the configuration DEST into SOURCE.
          void (*gsl_siman_copy_t) (void *source, void *dest)

 - Data Type: gsl_siman_copy_construct_t
     This function type should create a new copy of the configuration
     XP.
          void * (*gsl_siman_copy_construct_t) (void *xp)

 - Data Type: gsl_siman_destroy_t
     This function type should destroy the configuration XP, freeing its
     memory.
          void (*gsl_siman_destroy_t) (void *xp)

 - Data Type: gsl_siman_params_t
     These are the parameters that control a run of `gsl_siman_solve'.
     This structure contains all the information needed to control the
     search, beyond the energy function, the step function and the
     initial guess.

    `int n_tries'
          The number of points to try for each step.

    `int iters_fixed_T'
          The number of iterations at each temperature.

    `double step_size'
          The maximum step size in the random walk.

    `double k, t_initial, mu_t, t_min'
          The parameters of the Boltzmann distribution and cooling
          schedule.


File: gsl-ref.info,  Node: Examples with Simulated Annealing,  Prev: Simulated Annealing functions,  Up: Simulated Annealing

Examples
========

The simulated annealing package is clumsy, and it has to be because it
is written in C, for C callers, and tries to be polymorphic at the same
time.  But here we provide some examples which can be pasted into your
application with little change and should make things easier.

* Menu:

* Trivial example::
* Traveling Salesman Problem::


File: gsl-ref.info,  Node: Trivial example,  Next: Traveling Salesman Problem,  Prev: Examples with Simulated Annealing,  Up: Examples with Simulated Annealing

Trivial example
---------------

The first example, in one dimensional cartesian space, sets up an energy
function which is a damped sine wave; this has many local minima, but
only one global minimum, somewhere between 1.0 and 1.5.  The initial
guess given is 15.5, which is several local minima away from the global
minimum.

     #include <math.h>
     #include <stdlib.h>
     #include <gsl/gsl_siman.h>
     
     /* set up parameters for this simulated annealing run */
     
     /* how many points do we try before stepping */
     #define N_TRIES 200             
     
     /* how many iterations for each T? */
     #define ITERS_FIXED_T 10        
     
     /* max step size in random walk */
     #define STEP_SIZE 10            
     
     /* Boltzmann constant */
     #define K 1.0                   
     
     /* initial temperature */
     #define T_INITIAL 0.002         
     
     /* damping factor for temperature */
     #define MU_T 1.005              
     #define T_MIN 2.0e-6
     
     gsl_siman_params_t params 
       = {N_TRIES, ITERS_FIXED_T, STEP_SIZE,
          K, T_INITIAL, MU_T, T_MIN};
     
     /* now some functions to test in one dimension */
     double E1(void *xp)
     {
       double x = * ((double *) xp);
     
       return exp(-pow((x-1.0),2.0))*sin(8*x);
     }
     
     double M1(void *xp, void *yp)
     {
       double x = *((double *) xp);
       double y = *((double *) yp);
     
       return fabs(x - y);
     }
     
     void S1(const gsl_rng * r, void *xp, double step_size)
     {
       double old_x = *((double *) xp);
       double new_x;
     
       double u = gsl_rng_uniform(r);
       new_x = u * 2 * step_size - step_size + old_x;
     
       memcpy(xp, &new_x, sizeof(new_x));
     }
     
     void P1(void *xp)
     {
       printf ("%12g", *((double *) xp));
     }
     
     int
     main(int argc, char *argv[])
     {
       const gsl_rng_type * T;
       gsl_rng * r;
     
       double x_initial = 15.5;
     
       gsl_rng_env_setup();
     
       T = gsl_rng_default;
       r = gsl_rng_alloc(T);
     
       gsl_siman_solve(r, &x_initial, E1, S1, M1, P1,
                       NULL, NULL, NULL, 
                       sizeof(double), params);
       return 0;
     }

   Here are a couple of plots that are generated by running
`siman_test' in the following way:
     ./siman_test | grep -v "^#"
       | xyplot -xyil -y -0.88 -0.83 -d "x...y"
       | xyps -d > siman-test.eps
     ./siman_test | grep -v "^#"
       | xyplot -xyil -xl "generation" -yl "energy" -d "x..y"
       | xyps -d > siman-energy.eps


File: gsl-ref.info,  Node: Traveling Salesman Problem,  Prev: Trivial example,  Up: Examples with Simulated Annealing

Traveling Salesman Problem
--------------------------

The TSP ("Traveling Salesman Problem") is the classic combinatorial
optimization problem.  I have provided a very simple version of it,
based on the coordinates of twelve cities in the southwestern United
States.  This should maybe be called the "Flying Salesman Problem",
since I am using the great-circle distance between cities, rather than
the driving distance.  Also: I assume the earth is a sphere, so I don't
use geoid distances.

   The `gsl_siman_solve()' routine finds a route which is 3490.62
Kilometers long; this is confirmed by an exhaustive search of all
possible routes with the same initial city.

   The full code can be found in `siman/siman_tsp.c', but I include
here some plots generated in the following way:
     ./siman_tsp > tsp.output
     grep -v "^#" tsp.output
       | xyplot -xyil -d "x................y"
                -lx "generation" -ly "distance"
                -lt "TSP -- 12 southwest cities"
       | xyps -d > 12-cities.eps
     grep initial_city_coord tsp.output
       | awk '{print $2, $3, $4, $5}'
       | xyplot -xyil -lb0 -cs 0.8
                -lx "longitude (- means west)"
                -ly "latitude"
                -lt "TSP -- initial-order"
       | xyps -d > initial-route.eps
     grep final_city_coord tsp.output
       | awk '{print $2, $3, $4, $5}'
       | xyplot -xyil -lb0 -cs 0.8
                -lx "longitude (- means west)"
                -ly "latitude"
                -lt "TSP -- final-order"
       | xyps -d > final-route.eps

   This is the output showing the initial order of the cities;
longitude is negative, since it is west and I want the plot to look
like a map.
     # initial coordinates of cities (longitude and latitude)
     ###initial_city_coord: -105.95 35.68 Santa Fe
     ###initial_city_coord: -112.07 33.54 Phoenix
     ###initial_city_coord: -106.62 35.12 Albuquerque
     ###initial_city_coord: -103.2 34.41 Clovis
     ###initial_city_coord: -107.87 37.29 Durango
     ###initial_city_coord: -96.77 32.79 Dallas
     ###initial_city_coord: -105.92 35.77 Tesuque
     ###initial_city_coord: -107.84 35.15 Grants
     ###initial_city_coord: -106.28 35.89 Los Alamos
     ###initial_city_coord: -106.76 32.34 Las Cruces
     ###initial_city_coord: -108.58 37.35 Cortez
     ###initial_city_coord: -108.74 35.52 Gallup
     ###initial_city_coord: -105.95 35.68 Santa Fe

   The optimal route turns out to be:
     # final coordinates of cities (longitude and latitude)
     ###final_city_coord: -105.95 35.68 Santa Fe
     ###final_city_coord: -106.28 35.89 Los Alamos
     ###final_city_coord: -106.62 35.12 Albuquerque
     ###final_city_coord: -107.84 35.15 Grants
     ###final_city_coord: -107.87 37.29 Durango
     ###final_city_coord: -108.58 37.35 Cortez
     ###final_city_coord: -108.74 35.52 Gallup
     ###final_city_coord: -112.07 33.54 Phoenix
     ###final_city_coord: -106.76 32.34 Las Cruces
     ###final_city_coord: -96.77 32.79 Dallas
     ###final_city_coord: -103.2 34.41 Clovis
     ###final_city_coord: -105.92 35.77 Tesuque
     ###final_city_coord: -105.95 35.68 Santa Fe

Here's a plot of the cost function (energy) versus generation (point in
the calculation at which a new temperature is set) for this problem:


File: gsl-ref.info,  Node: Ordinary Differential Equations,  Next: Interpolation,  Prev: Simulated Annealing,  Up: Top

Ordinary Differential Equations
*******************************

This chapter describes functions for solving ordinary differential
equation (ODE) initial value problems.  The library provides a variety
of low-level methods, such as Runge-Kutta and Bulirsch-Stoer routines,
and higher-level components for adaptive step-size control.  The
components can be combined by the user to achieve the desired solution,
with full access to any intermediate steps.

   These functions are declared in the header file `gsl_odeiv.h'.

* Menu:

* Defining the ODE System::
* Stepping Functions::
* Adaptive Step-size Control::
* Evolution::
* ODE Example programs::
* ODE References and Further Reading::


File: gsl-ref.info,  Node: Defining the ODE System,  Next: Stepping Functions,  Up: Ordinary Differential Equations

Defining the ODE System
=======================

The routines solve the general n-dimensional first-order system,

     dy_i(t)/dt = f_i(t, y_1(t), ..., y_n(t))

for i = 1, \dots, n.  The stepping functions rely on the vector of
derivatives f_i and the Jacobian matrix, J_{ij} = df_i(t,y(t)) / dy_j.
A system of equations is defined using the `gsl_odeiv_system' datatype.

 - Data Type: gsl_odeiv_system
     This data type defines a general ODE system with arbitrary
     parameters.

    `int (* function) (double t, const double y[], double dydt[], void * params)'
          This function should store the vector elements
          f_i(t,y,params) in the array DYDT, for arguments (T,Y) and
          parameters PARAMS

    `int (* jacobian) (double t, const double y[], double * dfdy, double dfdt[], void * params);'
          This function should store the vector elements
          df_i(t,y,params)/dt in the array DFDT and the Jacobian matrix
          J_{ij} in the array DFDY regarded as a row-ordered matrix
          `J(i,j) = dfdy[i * dimension + j]' where `dimension' is the
          dimension of the system.

          Some of the simpler solver algorithms do not make use of the
          Jacobian matrix, so it is not always strictly necessary to
          provide it (the `jacobian' element of the struct can be
          replaced by a null pointer for those algorithms).  However,
          it is useful to provide the Jacobian to allow the solver
          algorithms to be interchanged--the best algorithms make use
          of the Jacobian.

    `size_t dimension;'
          This is the dimension of the system of equations.

    `void * params'
          This is a pointer to the arbitrary parameters of the system.


File: gsl-ref.info,  Node: Stepping Functions,  Next: Adaptive Step-size Control,  Prev: Defining the ODE System,  Up: Ordinary Differential Equations

Stepping Functions
==================

The lowest level components are the "stepping functions" which advance
a solution from time t to t+h for a fixed step-size h and estimate the
resulting local error.

 - Function: gsl_odeiv_step * gsl_odeiv_step_alloc (const
          gsl_odeiv_step_type * T, size_t DIM)
     This function returns a pointer to a newly allocated instance of a
     stepping function of type T for a system of DIM dimensions.

 - Function: int gsl_odeiv_step_reset (gsl_odeiv_step * S)
     This function resets the stepping function S.  It should be used
     whenever the next use of S will not be a continuation of a
     previous step.

 - Function: void gsl_odeiv_step_free (gsl_odeiv_step * S)
     This function frees all the memory associated with the stepping
     function S.

 - Function: const char * gsl_odeiv_step_name (const gsl_odeiv_step * S)
     This function returns a pointer to the name of the stepping
     function.  For example,

          printf ("step method is '%s'\n",
                   gsl_odeiv_step_name (s));

     would print something like `step method is 'rk4''.

 - Function: unsigned int gsl_odeiv_step_order (const gsl_odeiv_step *
          S)
     This function returns the order of the stepping function on the
     previous step.  This order can vary if the stepping function
     itself is adaptive.

 - Function: int gsl_odeiv_step_apply (gsl_odeiv_step * S, double T,
          double H, double Y[], double YERR[], const double DYDT_IN[],
          double DYDT_OUT[], const gsl_odeiv_system * DYDT)
     This function applies the stepping function S to the system of
     equations defined by DYDT, using the step size H to advance the
     system from time T and state Y to time T+H.  The new state of the
     system is stored in Y on output, with an estimate of the absolute
     error in each component stored in YERR.  If the argument DYDT_IN
     is not null it should point an array containing the derivatives
     for the system at time T on input. This is optional as the
     derivatives will be computed internally if they are not provided,
     but allows the reuse of existing derivative information.  On
     output the new derivatives of the system at time T+H will be
     stored in DYDT_OUT if it is not null.

   The following algorithms are available,

 - Step Type: gsl_odeiv_step_rk2
     Embedded Runge-Kutta (2, 3) method.

 - Step Type: gsl_odeiv_step_rk4
     4th order (classical) Runge-Kutta.

 - Step Type: gsl_odeiv_step_rkf45
     Embedded Runge-Kutta-Fehlberg (4, 5) method.  This method is a good
     general-purpose integrator.

 - Step Type: gsl_odeiv_step_rkck
     Embedded Runge-Kutta Cash-Karp (4, 5) method.

 - Step Type: gsl_odeiv_step_rk8pd
     Embedded Runge-Kutta Prince-Dormand (8,9) method.

 - Step Type: gsl_odeiv_step_rk2imp
     Implicit 2nd order Runge-Kutta at Gaussian points.

 - Step Type: gsl_odeiv_step_rk4imp
     Implicit 4th order Runge-Kutta at Gaussian points.

 - Step Type: gsl_odeiv_step_bsimp
     Implicit Bulirsch-Stoer method of Bader and Deuflhard.  This
     algorithm requires the Jacobian.

 - Step Type: gsl_odeiv_step_gear1
     M=1 implicit Gear method.

 - Step Type: gsl_odeiv_step_gear2
     M=2 implicit Gear method.


File: gsl-ref.info,  Node: Adaptive Step-size Control,  Next: Evolution,  Prev: Stepping Functions,  Up: Ordinary Differential Equations

Adaptive Step-size Control
==========================

The control function examines the proposed change to the solution and
its error estimate produced by a stepping function and attempts to
determine the optimal step-size for a user-specified level of error.

 - Function: gsl_odeiv_control * gsl_odeiv_control_standard_new (double
          EPS_ABS, double EPS_REL, double A_Y, double A_DYDT)
     The standard control object is a four parameter heuristic based on
     absolute and relative errors EPS_ABS and EPS_REL, and scaling
     factors A_Y and A_DYDT for the system state y(t) and derivatives
     y'(t) respectively.

     The step-size adjustment procedure for this method begins by
     computing the desired error level D_i for each component,

          D_i = eps_abs + eps_rel * (a_y |y_i| + a_dydt h |y'_i|)

     and comparing it with the observed error E_i = |yerr_i|.  If the
     observed error E exceeds the desired error level D by more than
     10% for any component then the method reduces the step-size by an
     appropriate factor,

          h_new = h_old * S * (E/D)^(-1/q)

     where q is the consistency order of method (e.g. q=4 for 4(5)
     embedded RK), and S is a safety factor of 0.9. The ratio E/D is
     taken to be the maximum of the ratios E_i/D_i.

     If the observed error E is less than 50% of the desired error
     level D for the maximum ratio E_i/D_i then the algorithm takes the
     opportunity to increase the step-size to bring the error in line
     with the desired level,

          h_new = h_old * S * (E/D)^(-1/(q+1))

     This encompasses all the standard error scaling methods. To avoid
     uncontrolled changes in the stepsize, the overall scaling factor is
     limited to the range 1/5 to 5.

 - Function: gsl_odeiv_control * gsl_odeiv_control_y_new (double
          EPS_ABS, double EPS_REL)
     This function creates a new control object which will keep the
     local error on each step within an absolute error of EPS_ABS and
     relative error of EPS_REL with respect to the solution y_i(t).
     This is equivalent to the standard control object with A_Y=1 and
     A_DYDT=0.

 - Function: gsl_odeiv_control * gsl_odeiv_control_yp_new (double
          EPS_ABS, double EPS_REL)
     This function creates a new control object which will keep the
     local error on each step within an absolute error of EPS_ABS and
     relative error of EPS_REL with respect to the derivatives of the
     solution y'_i(t).  This is equivalent to the standard control
     object with A_Y=0 and A_DYDT=1.

 - Function: gsl_odeiv_control * gsl_odeiv_control_scaled_new (double
          EPS_ABS, double EPS_REL, double A_Y, double A_DYDT, const
          double SCALE_ABS[], size_t DIM)
     This function creates a new control object which uses the same
     algorithm as `gsl_odeiv_control_standard_new' but with an absolute
     error which is scaled for each component by the array SCALE_ABS.
     The formula for D_i for this control object is,

          D_i = eps_abs * s_i + eps_rel * (a_y |y_i| + a_dydt h |y'_i|)

     where s_i is the i-th component of the array SCALE_ABS.  The same
     error control heuristic is used by the Matlab ODE suite.

 - Function: gsl_odeiv_control * gsl_odeiv_control_alloc (const
          gsl_odeiv_control_type * T)
     This function returns a pointer to a newly allocated instance of a
     control function of type T.  This function is only needed for
     defining new types of control functions.  For most purposes the
     standard control functions described above should be sufficient.

 - Function: int gsl_odeiv_control_init (gsl_odeiv_control * C, double
          EPS_ABS, double EPS_REL, double A_Y, double A_DYDT)
     This function initializes the control function C with the
     parameters EPS_ABS (absolute error), EPS_REL (relative error), A_Y
     (scaling factor for y) and A_DYDT (scaling factor for derivatives).

 - Function: void gsl_odeiv_control_free (gsl_odeiv_control * C)
     This function frees all the memory associated with the control
     function C.

 - Function: int gsl_odeiv_control_hadjust (gsl_odeiv_control * C,
          gsl_odeiv_step * S, const double Y0[], const double YERR[],
          const double DYDT[], double * H)
     This function adjusts the step-size H using the control function
     C, and the current values of Y, YERR and DYDT.  The stepping
     function STEP is also needed to determine the order of the method.
     If the error in the y-values YERR is found to be too large then
     the step-size H is reduced and the function returns
     `GSL_ODEIV_HADJ_DEC'.  If the error is sufficiently small then H
     may be increased and `GSL_ODEIV_HADJ_INC' is returned.  The
     function returns `GSL_ODEIV_HADJ_NIL' if the step-size is
     unchanged.  The goal of the function is to estimate the largest
     step-size which satisfies the user-specified accuracy requirements
     for the current point.

 - Function: const char * gsl_odeiv_control_name (const
          gsl_odeiv_control * C)
     This function returns a pointer to the name of the control
     function.  For example,

          printf ("control method is '%s'\n",
                  gsl_odeiv_control_name (c));

     would print something like `control method is 'standard''


File: gsl-ref.info,  Node: Evolution,  Next: ODE Example programs,  Prev: Adaptive Step-size Control,  Up: Ordinary Differential Equations

Evolution
=========

The highest level of the system is the evolution function which combines
the results of a stepping function and control function to reliably
advance the solution forward over an interval (t_0, t_1).  If the
control function signals that the step-size should be decreased the
evolution function backs out of the current step and tries the proposed
smaller step-size.  This process is continued until an acceptable
step-size is found.

 - Function: gsl_odeiv_evolve * gsl_odeiv_evolve_alloc (size_t DIM)
     This function returns a pointer to a newly allocated instance of an
     evolution function for a system of DIM dimensions.

 - Function: int gsl_odeiv_evolve_apply (gsl_odeiv_evolve * E,
          gsl_odeiv_control * CON, gsl_odeiv_step * STEP, const
          gsl_odeiv_system * DYDT, double * T, double T1, double * H,
          double Y[])
     This function advances the system (E, DYDT) from time T and
     position Y using the stepping function STEP.  The new time and
     position are stored in T and Y on output.  The initial step-size
     is taken as H, but this will be modified using the control
     function C to achieve the appropriate error bound if necessary.
     The routine may make several calls to STEP in order to determine
     the optimum step-size. If the step-size has been changed the value
     of H will be modified on output.  The maximum time T1 is
     guaranteed not to be exceeded by the time-step.  On the final
     time-step the value of T will be set to T1 exactly.

 - Function: int gsl_odeiv_evolve_reset (gsl_odeiv_evolve * E)
     This function resets the evolution function E.  It should be used
     whenever the next use of E will not be a continuation of a
     previous step.

 - Function: void gsl_odeiv_evolve_free (gsl_odeiv_evolve * E)
     This function frees all the memory associated with the evolution
     function E.


File: gsl-ref.info,  Node: ODE Example programs,  Next: ODE References and Further Reading,  Prev: Evolution,  Up: Ordinary Differential Equations

Examples
========

The following program solves the second-order nonlinear Van der Pol
oscillator equation,

     x''(t) + \mu x'(t) (x(t)^2 - 1) + x(t) = 0

This can be converted into a first order system suitable for use with
the routines described in this chapter by introducing a separate
variable for the velocity, y = x'(t),

     x' = y
     y' = -x + \mu y (1-x^2)

The program begins by defining functions for these derivatives and
their Jacobian,

     #include <stdio.h>
     #include <gsl/gsl_errno.h>
     #include <gsl/gsl_matrix.h>
     #include <gsl/gsl_odeiv.h>
     
     int
     func (double t, const double y[], double f[],
           void *params)
     {
       double mu = *(double *)params;
       f[0] = y[1];
       f[1] = -y[0] - mu*y[1]*(y[0]*y[0] - 1);
       return GSL_SUCCESS;
     }
     
     int
     jac (double t, const double y[], double *dfdy, 
          double dfdt[], void *params)
     {
       double mu = *(double *)params;
       gsl_matrix_view dfdy_mat 
         = gsl_matrix_view_array (dfdy, 2, 2);
       gsl_matrix * m = &dfdy_mat.matrix; 
       gsl_matrix_set (m, 0, 0, 0.0);
       gsl_matrix_set (m, 0, 1, 1.0);
       gsl_matrix_set (m, 1, 0, -2.0*mu*y[0]*y[1] - 1.0);
       gsl_matrix_set (m, 1, 1, -mu*(y[0]*y[0] - 1.0));
       dfdt[0] = 0.0;
       dfdt[1] = 0.0;
       return GSL_SUCCESS;
     }
     
     int
     main (void)
     {
       const gsl_odeiv_step_type * T 
         = gsl_odeiv_step_rk8pd;
     
       gsl_odeiv_step * s 
         = gsl_odeiv_step_alloc (T, 2);
       gsl_odeiv_control * c 
         = gsl_odeiv_control_y_new (1e-6, 0.0);
       gsl_odeiv_evolve * e 
         = gsl_odeiv_evolve_alloc (2);
     
       double mu = 10;
       gsl_odeiv_system sys = {func, jac, 2, &mu};
     
       double t = 0.0, t1 = 100.0;
       double h = 1e-6;
       double y[2] = { 1.0, 0.0 };
     
       while (t < t1)
         {
           int status = gsl_odeiv_evolve_apply (e, c, s,
                                                &sys, 
                                                &t, t1,
                                                &h, y);
     
           if (status != GSL_SUCCESS)
               break;
     
           printf ("%.5e %.5e %.5e\n", t, y[0], y[1]);
         }
     
       gsl_odeiv_evolve_free (e);
       gsl_odeiv_control_free (c);
       gsl_odeiv_step_free (s);
       return 0;
     }

The main loop of the program evolves the solution from (y, y') = (1, 0)
at t=0 to t=100.  The step-size h is automatically adjusted by the
controller to maintain an absolute accuracy of 10^{-6} in the function
values Y.

To obtain the values at regular intervals, rather than the variable
spacings chosen by the control function, the main loop can be modified
to advance the solution from one point to the next.  For example, the
following main loop prints the solution at the fixed points t = 0, 1,
2, \dots, 100,

       for (i = 1; i <= 100; i++)
         {
           double ti = i * t1 / 100.0;
     
           while (t < ti)
             {
               gsl_odeiv_evolve_apply (e, c, s,
                                       &sys,
                                       &t, ti, &h,
                                       y);
             }
     
           printf ("%.5e %.5e %.5e\n", t, y[0], y[1]);
         }

It is also possible to work with a non-adaptive integrator, using only
the stepping function itself.  The following program uses the `rk4'
fourth-order Runge-Kutta stepping function with a fixed stepsize of
0.01,

     int
     main (void)
     {
       const gsl_odeiv_step_type * T 
         = gsl_odeiv_step_rk4;
     
       gsl_odeiv_step * s 
         = gsl_odeiv_step_alloc (T, 2);
     
       double mu = 10;
       gsl_odeiv_system sys = {func, jac, 2, &mu};
     
       double t = 0.0, t1 = 100.0;
       double h = 1e-2;
       double y[2] = { 1.0, 0.0 }, y_err[2];
       double dydt_in[2], dydt_out[2];
     
       /* initialise dydt_in from system parameters */
       GSL_ODEIV_FN_EVAL(&sys, t, y, dydt_in);
     
       while (t < t1)
         {
           int status = gsl_odeiv_step_apply (s, t, h, 
                                              y, y_err, 
                                              dydt_in, 
                                              dydt_out, 
                                              &sys);
     
           if (status != GSL_SUCCESS)
               break;
     
           dydt_in[0] = dydt_out[0];
           dydt_in[1] = dydt_out[1];
     
           t += h;
     
           printf ("%.5e %.5e %.5e\n", t, y[0], y[1]);
         }
     
       gsl_odeiv_step_free (s);
       return 0;
     }

The derivatives must be initialised for the starting point t=0 before
the first step is taken.  Subsequent steps use the output derivatives
DYDT_OUT as inputs to the next step by copying their values into
DYDT_IN.


File: gsl-ref.info,  Node: ODE References and Further Reading,  Prev: ODE Example programs,  Up: Ordinary Differential Equations

References and Further Reading
==============================

Many of the basic Runge-Kutta formulas can be found in the Handbook of
Mathematical Functions,

     Abramowitz & Stegun (eds.), `Handbook of Mathematical Functions',
     Section 25.5.

The implicit Bulirsch-Stoer algorithm `bsimp' is described in the
following paper,

     G. Bader and P. Deuflhard, "A Semi-Implicit Mid-Point Rule for
     Stiff Systems of Ordinary Differential Equations.", Numer. Math.
     41, 373-398, 1983.


File: gsl-ref.info,  Node: Interpolation,  Next: Numerical Differentiation,  Prev: Ordinary Differential Equations,  Up: Top

Interpolation
*************

This chapter describes functions for performing interpolation.  The
library provides a variety of interpolation methods, including Cubic
splines and Akima splines.  The interpolation types are interchangeable,
allowing different methods to be used without recompiling.
Interpolations can be defined for both normal and periodic boundary
conditions.  Additional functions are available for computing
derivatives and integrals of interpolating functions.

   The functions described in this section are declared in the header
files `gsl_interp.h' and `gsl_spline.h'.

* Menu:

* Introduction to Interpolation::
* Interpolation Functions::
* Interpolation Types::
* Index Look-up and Acceleration::
* Evaluation of Interpolating Functions::
* Higher-level Interface::
* Interpolation Example programs::
* Interpolation References and Further Reading::


File: gsl-ref.info,  Node: Introduction to Interpolation,  Next: Interpolation Functions,  Up: Interpolation

Introduction
============

Given a set of data points (x_1, y_1) \dots (x_n, y_n) the routines
described in this section compute a continuous interpolating function
y(x) such that y_i = y(x_i).  The interpolation is piecewise smooth,
and its behavior at the end-points is determined by the type of
interpolation used.


File: gsl-ref.info,  Node: Interpolation Functions,  Next: Interpolation Types,  Prev: Introduction to Interpolation,  Up: Interpolation

Interpolation Functions
=======================

The interpolation function for a given dataset is stored in a
`gsl_interp' object.  These are created by the following functions.

 - Function: gsl_interp * gsl_interp_alloc (const gsl_interp_type * T,
          size_t SIZE)
     This function returns a pointer to a newly allocated interpolation
     object of type T for SIZE data-points.

 - Function: int gsl_interp_init (gsl_interp * INTERP, const double
          XA[], const double YA[], size_t SIZE)
     This function initializes the interpolation object INTERP for the
     data (XA,YA) where XA and YA are arrays of size SIZE.  The
     interpolation object (`gsl_interp') does not save the data arrays
     XA and YA and only stores the static state computed from the data.
     The XA data array is always assumed to be strictly ordered; the
     behavior for other arrangements is not defined.

 - Function: void gsl_interp_free (gsl_interp * INTERP)
     This function frees the interpolation object INTERP.


File: gsl-ref.info,  Node: Interpolation Types,  Next: Index Look-up and Acceleration,  Prev: Interpolation Functions,  Up: Interpolation

Interpolation Types
===================

The interpolation library provides five interpolation types:

 - Interpolation Type: gsl_interp_linear
     Linear interpolation.  This interpolation method does not require
     any additional memory.

 - Interpolation Type: gsl_interp_polynomial
     Polynomial interpolation.  This method should only be used for
     interpolating small numbers of points because polynomial
     interpolation introduces large oscillations, even for well-behaved
     datasets.  The number of terms in the interpolating polynomial is
     equal to the number of points.

 - Interpolation Type: gsl_interp_cspline
     Cubic spline with natural boundary conditions.

 - Interpolation Type: gsl_interp_cspline_periodic
     Cubic spline with periodic boundary conditions

 - Interpolation Type: gsl_interp_akima
     Non-rounded Akima spline with natural boundary conditions.  This
     method uses the non-rounded corner algorithm of Wodicka.

 - Interpolation Type: gsl_interp_akima_periodic
     Non-rounded Akima spline with periodic boundary conditions.  This
     method uses the non-rounded corner algorithm of Wodicka.

   The following related functions are available:

 - Function: const char * gsl_interp_name (const gsl_interp * INTERP)
     This function returns the name of the interpolation type used by
     INTERP.  For example,

          printf ("interp uses '%s' interpolation.\n",
                  gsl_interp_name (interp));

     would print something like,
          interp uses 'cspline' interpolation.

 - Function: unsigned int gsl_interp_min_size (const gsl_interp *
          INTERP)
     This function returns the minimum number of points required by the
     interpolation type of INTERP.  For example, Akima spline
     interpolation requires a minimum of 5 points.


File: gsl-ref.info,  Node: Index Look-up and Acceleration,  Next: Evaluation of Interpolating Functions,  Prev: Interpolation Types,  Up: Interpolation

Index Look-up and Acceleration
==============================

The state of searches can be stored in a `gsl_interp_accel' object,
which is a kind of iterator for interpolation lookups.  It caches the
previous value of an index lookup.  When the subsequent interpolation
point falls in the same interval its index value can be returned
immediately.

 - Function: size_t gsl_interp_bsearch (const double X_ARRAY[], double
          X, size_t INDEX_LO, size_t INDEX_HI)
     This function returns the index i of the array X_ARRAY such that
     `x_array[i] <= x < x_array[i+1]'.  The index is searched for in
     the range [INDEX_LO,INDEX_HI].

 - Function: gsl_interp_accel * gsl_interp_accel_alloc (void)
     This function returns a pointer to an accelerator object, which is
     a kind of iterator for interpolation lookups.  It tracks the state
     of lookups, thus allowing for application of various acceleration
     strategies.

 - Function: size_t gsl_interp_accel_find (gsl_interp_accel * A, const
          double X_ARRAY[], size_t SIZE, double X)
     This function performs a lookup action on the data array X_ARRAY
     of size SIZE, using the given accelerator A.  This is how lookups
     are performed during evaluation of an interpolation.  The function
     returns an index i such that `x_array[i] <= x < x_array[i+1]'.

 - Function: void gsl_interp_accel_free (gsl_interp_accel* ACC)
     This function frees the accelerator object ACC.


File: gsl-ref.info,  Node: Evaluation of Interpolating Functions,  Next: Higher-level Interface,  Prev: Index Look-up and Acceleration,  Up: Interpolation

Evaluation of Interpolating Functions
=====================================

 - Function: double gsl_interp_eval (const gsl_interp * INTERP, const
          double XA[], const double YA[], double X, gsl_interp_accel *
          ACC)
 - Function: int gsl_interp_eval_e (const gsl_interp * INTERP, const
          double XA[], const double YA[], double X, gsl_interp_accel *
          ACC, double * Y)
     These functions return the interpolated value of Y for a given
     point X, using the interpolation object INTERP, data arrays XA and
     YA and the accelerator ACC.

 - Function: double gsl_interp_eval_deriv (const gsl_interp * INTERP,
          const double XA[], const double YA[], double X,
          gsl_interp_accel * ACC)
 - Function: int gsl_interp_eval_deriv_e (const gsl_interp * INTERP,
          const double XA[], const double YA[], double X,
          gsl_interp_accel * ACC, double * D)
     These functions return the derivative D of an interpolated
     function for a given point X, using the interpolation object
     INTERP, data arrays XA and YA and the accelerator ACC.

 - Function: double gsl_interp_eval_deriv2 (const gsl_interp * INTERP,
          const double XA[], const double YA[], double X,
          gsl_interp_accel * ACC)
 - Function: int gsl_interp_eval_deriv2_e (const gsl_interp * INTERP,
          const double XA[], const double YA[], double X,
          gsl_interp_accel * ACC, double * D2)
     These functions return the second derivative D2 of an interpolated
     function for a given point X, using the interpolation object
     INTERP, data arrays XA and YA and the accelerator ACC.

 - Function: double gsl_interp_eval_integ (const gsl_interp * INTERP,
          const double XA[], const double YA[], double A, double B,
          gsl_interp_accel * ACC)
 - Function: int gsl_interp_eval_integ_e (const gsl_interp * INTERP,
          const double XA[], const double YA[], double A, double B,
          gsl_interp_accel * ACC, double * RESULT)
     These functions return the numerical integral RESULT of an
     interpolated function over the range [A, B], using the
     interpolation object INTERP, data arrays XA and YA and the
     accelerator ACC.


File: gsl-ref.info,  Node: Higher-level Interface,  Next: Interpolation Example programs,  Prev: Evaluation of Interpolating Functions,  Up: Interpolation

Higher-level Interface
======================

The functions described in the previous sections required the user to
supply pointers to the x and y arrays on each call.  The following
functions are equivalent to the corresponding `gsl_interp' functions
but maintain a copy of this data in the `gsl_spline' object.  This
removes the need to pass both XA and YA as arguments on each
evaluation. These functions are defined in the header file
`gsl_spline.h'.

 - Function: gsl_spline * gsl_spline_alloc (const gsl_interp_type * T,
          size_t SIZE)

 - Function: int gsl_spline_init (gsl_spline * SPLINE, const double
          XA[], const double YA[], size_t SIZE)

 - Function: void gsl_spline_free (gsl_spline * SPLINE)

 - Function: double gsl_spline_eval (const gsl_spline * SPLINE, double
          X, gsl_interp_accel * ACC)
 - Function: int gsl_spline_eval_e (const gsl_spline * SPLINE, double
          X, gsl_interp_accel * ACC, double * Y)

 - Function: double gsl_spline_eval_deriv (const gsl_spline * SPLINE,
          double X, gsl_interp_accel * ACC)
 - Function: int gsl_spline_eval_deriv_e (const gsl_spline * SPLINE,
          double X, gsl_interp_accel * ACC, double * D)

 - Function: double gsl_spline_eval_deriv2 (const gsl_spline * SPLINE,
          double X, gsl_interp_accel * ACC)
 - Function: int gsl_spline_eval_deriv2_e (const gsl_spline * SPLINE,
          double X, gsl_interp_accel * ACC, double * D2)

 - Function: double gsl_spline_eval_integ (const gsl_spline * SPLINE,
          double A, double B, gsl_interp_accel * ACC)
 - Function: int gsl_spline_eval_integ_e (const gsl_spline * SPLINE,
          double A, double B, gsl_interp_accel * ACC, double * RESULT)


File: gsl-ref.info,  Node: Interpolation Example programs,  Next: Interpolation References and Further Reading,  Prev: Higher-level Interface,  Up: Interpolation

Examples
========

The following program demonstrates the use of the interpolation and
spline functions.  It computes a cubic spline interpolation of the
10-point dataset (x_i, y_i) where x_i = i + \sin(i)/2 and y_i = i +
\cos(i^2) for i = 0 \dots 9.

     #include <stdlib.h>
     #include <stdio.h>
     #include <math.h>
     #include <gsl/gsl_errno.h>
     #include <gsl/gsl_spline.h>
     
     int
     main (void)
     {
       int i;
       double xi, yi, x[10], y[10];
     
       printf ("#m=0,S=2\n");
     
       for (i = 0; i < 10; i++)
         {
           x[i] = i + 0.5 * sin (i);
           y[i] = i + cos (i * i);
           printf ("%g %g\n", x[i], y[i]);
         }
     
       printf ("#m=1,S=0\n");
     
       {
         gsl_interp_accel *acc 
           = gsl_interp_accel_alloc ();
         gsl_spline *spline 
           = gsl_spline_alloc (gsl_interp_cspline, 10);
     
         gsl_spline_init (spline, x, y, 10);
     
         for (xi = x[0]; xi < x[9]; xi += 0.01)
           {
             yi = gsl_spline_eval (spline, xi, acc);
             printf ("%g %g\n", xi, yi);
           }
         gsl_spline_free (spline);
         gsl_interp_accel_free (acc);
       }
       return 0;
     }

The output is designed to be used with the GNU plotutils `graph'
program,

     $ ./a.out > interp.dat
     $ graph -T ps < interp.dat > interp.ps

The result shows a smooth interpolation of the original points.  The
interpolation method can changed simply by varying the first argument of
`gsl_spline_alloc'.


File: gsl-ref.info,  Node: Interpolation References and Further Reading,  Prev: Interpolation Example programs,  Up: Interpolation

References and Further Reading
==============================

Descriptions of the interpolation algorithms and further references can
be found in the following book,

     C.W. Ueberhuber, `Numerical Computation (Volume 1), Chapter 9
     "Interpolation"', Springer (1997), ISBN 3-540-62058-3.


File: gsl-ref.info,  Node: Numerical Differentiation,  Next: Chebyshev Approximations,  Prev: Interpolation,  Up: Top

Numerical Differentiation
*************************

The functions described in this chapter compute numerical derivatives by
finite differencing.  An adaptive algorithm is used to find the best
choice of finite difference and to estimate the error in the derivative.
These functions are declared in the header file `gsl_deriv.h'.

* Menu:

* Numerical Differentiation functions::
* Numerical Differentiation Examples::
* Numerical Differentiation References::


File: gsl-ref.info,  Node: Numerical Differentiation functions,  Next: Numerical Differentiation Examples,  Up: Numerical Differentiation

Functions
=========

 - Function: int gsl_deriv_central (const gsl_function *F, double X,
          double H, double *RESULT, double *ABSERR)
     This function computes the numerical derivative of the function F
     at the point X using an adaptive central difference algorithm with
     a step-size of H.   The derivative is returned in RESULT and an
     estimate of its absolute error is returned in ABSERR.

     The initial value of H is used to estimate an optimal step-size,
     based on the scaling of the truncation error and round-off error
     in the derivative calculation.  The derivative is computed using a
     5-point rule for equally spaced abscissae at x-h, x-h/2, x, x+h/2,
     x, with an error estimate taken from the difference between the
     5-point rule and the corresponding 3-point rule x-h, x, x+h.  Note
     that the value of the function at x does not contribute to the
     derivative calculation, so only 4-points are actually used.

 - Function: int gsl_deriv_forward (const gsl_function *F, double X,
          double H, double *RESULT, double *ABSERR)
     This function computes the numerical derivative of the function F
     at the point X using an adaptive forward difference algorithm with
     a step-size of H. The function is evaluated only at points greater
     than X, and never at X itself.  The derivative is returned in
     RESULT and an estimate of its absolute error is returned in
     ABSERR.  This function should be used if f(x) has a discontinuity
     at X, or is undefined for values less than X.

     The initial value of H is used to estimate an optimal step-size,
     based on the scaling of the truncation error and round-off error
     in the derivative calculation.  The derivative at x is computed
     using an "open" 4-point rule for equally spaced abscissae at x+h/4,
     x+h/2, x+3h/4, x+h, with an error estimate taken from the
     difference between the 4-point rule and the corresponding 2-point
     rule x+h/2, x+h.

 - Function: int gsl_deriv_backward (const gsl_function *F, double X,
          double H, double *RESULT, double *ABSERR)
     This function computes the numerical derivative of the function F
     at the point X using an adaptive backward difference algorithm
     with a step-size of H. The function is evaluated only at points
     less than X, and never at X itself.  The derivative is returned in
     RESULT and an estimate of its absolute error is returned in
     ABSERR.  This function should be used if f(x) has a discontinuity
     at X, or is undefined for values greater than X.

     This function is equivalent to calling `gsl_deriv_forward' with a
     negative step-size.


File: gsl-ref.info,  Node: Numerical Differentiation Examples,  Next: Numerical Differentiation References,  Prev: Numerical Differentiation functions,  Up: Numerical Differentiation

Examples
========

The following code estimates the derivative of the function f(x) =
x^{3/2} at x=2 and at x=0.  The function f(x) is undefined for x<0 so
the derivative at x=0 is computed using `gsl_deriv_forward'.

     #include <stdio.h>
     #include <gsl/gsl_math.h>
     #include <gsl/gsl_deriv.h>
     
     double f (double x, void * params)
     {
       return pow (x, 1.5);
     }
     
     int
     main (void)
     {
       gsl_function F;
       double result, abserr;
     
       F.function = &f;
       F.params = 0;
     
       printf ("f(x) = x^(3/2)\n");
     
       gsl_deriv_central (&F, 2.0, 1e-8, &result, &abserr);
       printf ("x = 2.0\n");
       printf ("f'(x) = %.10f +/- %.10f\n", result, abserr);
       printf ("exact = %.10f\n\n", 1.5 * sqrt(2.0));
     
       gsl_deriv_forward (&F, 0.0, 1e-8, &result, &abserr);
       printf ("x = 0.0\n");
       printf ("f'(x) = %.10f +/- %.10f\n", result, abserr);
       printf ("exact = %.10f\n", 0.0);
     
       return 0;
     }

Here is the output of the program,

     $ ./a.out

     f(x) = x^(3/2)
     x = 2.0
     f'(x) = 2.1213203120 +/- 0.0000004064
     exact = 2.1213203436
     
     x = 0.0
     f'(x) = 0.0000000160 +/- 0.0000000339
     exact = 0.0000000000


File: gsl-ref.info,  Node: Numerical Differentiation References,  Prev: Numerical Differentiation Examples,  Up: Numerical Differentiation

References and Further Reading
==============================

The algorithms used by these functions are described in the following
sources:

     Abramowitz and Stegun, `Handbook of Mathematical Functions',
     Section 25.3.4, and Table 25.5 (Coefficients for Differentiation).

     S.D. Conte and Carl de Boor, `Elementary Numerical Analysis: An
     Algorithmic Approach', McGraw-Hill, 1972.


File: gsl-ref.info,  Node: Chebyshev Approximations,  Next: Series Acceleration,  Prev: Numerical Differentiation,  Up: Top

Chebyshev Approximations
************************

This chapter describes routines for computing Chebyshev approximations
to univariate functions.  A Chebyshev approximation is a truncation of
the series f(x) = \sum c_n T_n(x), where the Chebyshev polynomials
T_n(x) = \cos(n \arccos x) provide an orthogonal basis of polynomials
on the interval [-1,1] with the weight function 1 / \sqrt{1-x^2}.  The
first few Chebyshev polynomials are, T_0(x) = 1, T_1(x) = x, T_2(x) = 2
x^2 - 1.  For further information see Abramowitz & Stegun, Chapter 22.

   The functions described in this chapter are declared in the header
file `gsl_chebyshev.h'.

* Menu:

* The gsl_cheb_series struct::
* Creation and Calculation of Chebyshev Series::
* Chebyshev Series Evaluation::
* Derivatives and Integrals::
* Chebyshev Approximation examples::
* Chebyshev Approximation References and Further Reading::


File: gsl-ref.info,  Node: The gsl_cheb_series struct,  Next: Creation and Calculation of Chebyshev Series,  Up: Chebyshev Approximations

The gsl_cheb_series struct
==========================

A Chebyshev series  is stored using the following structure,

     typedef struct
     {
       double * c;   /* coefficients  c[0] .. c[order] */
       int order;    /* order of expansion             */
       double a;     /* lower interval point           */
       double b;     /* upper interval point           */
     } gsl_cheb_struct

The approximation is made over the range [a,b] using ORDER+1 terms,
including the coefficient c[0].  The series is computed using the
following convention,

     f(x) = (c_0 / 2) + \sum_{n=1} c_n T_n(x)

which is needed when accessing the coefficients directly.


File: gsl-ref.info,  Node: Creation and Calculation of Chebyshev Series,  Next: Chebyshev Series Evaluation,  Prev: The gsl_cheb_series struct,  Up: Chebyshev Approximations

Creation and Calculation of Chebyshev Series
============================================

 - Function: gsl_cheb_series * gsl_cheb_alloc (const size_t N)
     This function allocates space for a Chebyshev series of order N
     and returns a pointer to a new `gsl_cheb_series' struct.

 - Function: void gsl_cheb_free (gsl_cheb_series * CS)
     This function frees a previously allocated Chebyshev series CS.

 - Function: int gsl_cheb_init (gsl_cheb_series * CS, const
          gsl_function * F, const double A, const double B)
     This function computes the Chebyshev approximation CS for the
     function F over the range (a,b) to the previously specified order.
     The computation of the Chebyshev approximation is an O(n^2)
     process, and requires n function evaluations.


File: gsl-ref.info,  Node: Chebyshev Series Evaluation,  Next: Derivatives and Integrals,  Prev: Creation and Calculation of Chebyshev Series,  Up: Chebyshev Approximations

Chebyshev Series Evaluation
===========================

 - Function: double gsl_cheb_eval (const gsl_cheb_series * CS, double X)
     This function evaluates the Chebyshev series CS at a given point X.

 - Function: int gsl_cheb_eval_err (const gsl_cheb_series * CS, const
          double X, double * RESULT, double * ABSERR)
     This function computes the Chebyshev series CS at a given point X,
     estimating both the series RESULT and its absolute error ABSERR.
     The error estimate is made from the first neglected term in the
     series.

 - Function: double gsl_cheb_eval_n (const gsl_cheb_series * CS, size_t
          ORDER, double X)
     This function evaluates the Chebyshev series CS at a given point
     N, to (at most) the given order ORDER.

 - Function: int gsl_cheb_eval_n_err (const gsl_cheb_series * CS, const
          size_t ORDER, const double X, double * RESULT, double *
          ABSERR)
     This function evaluates a Chebyshev series CS at a given point X,
     estimating both the series RESULT and its absolute error ABSERR,
     to (at most) the given order ORDER.  The error estimate is made
     from the first neglected term in the series.


File: gsl-ref.info,  Node: Derivatives and Integrals,  Next: Chebyshev Approximation examples,  Prev: Chebyshev Series Evaluation,  Up: Chebyshev Approximations

Derivatives and Integrals
=========================

The following functions allow a Chebyshev series to be differentiated or
integrated, producing a new Chebyshev series.  Note that the error
estimate produced by evaluating the derivative series will be
underestimated due to the contribution of higher order terms being
neglected.

 - Function: int gsl_cheb_calc_deriv (gsl_cheb_series * DERIV, const
          gsl_cheb_series * CS)
     This function computes the derivative of the series CS, storing
     the derivative coefficients in the previously allocated DERIV.
     The two series CS and DERIV must have been allocated with the same
     order.

 - Function: int gsl_cheb_calc_integ (gsl_cheb_series * INTEG, const
          gsl_cheb_series * CS)
     This function computes the integral of the series CS, storing the
     integral coefficients in the previously allocated INTEG.  The two
     series CS and INTEG must have been allocated with the same order.
     The lower limit of the integration is taken to be the left hand
     end of the range A.


File: gsl-ref.info,  Node: Chebyshev Approximation examples,  Next: Chebyshev Approximation References and Further Reading,  Prev: Derivatives and Integrals,  Up: Chebyshev Approximations

Examples
========

The following example program computes Chebyshev approximations to a
step function.  This is an extremely difficult approximation to make,
due to the discontinuity, and was chosen as an example where
approximation error is visible.  For smooth functions the Chebyshev
approximation converges extremely rapidly and errors would not be
visible.

     #include <stdio.h>
     #include <gsl/gsl_math.h>
     #include <gsl/gsl_chebyshev.h>
     
     double
     f (double x, void *p)
     {
       if (x < 0.5)
         return 0.25;
       else
         return 0.75;
     }
     
     int
     main (void)
     {
       int i, n = 10000; 
     
       gsl_cheb_series *cs = gsl_cheb_alloc (40);
     
       gsl_function F;
     
       F.function = f;
       F.params = 0;
     
       gsl_cheb_init (cs, &F, 0.0, 1.0);
     
       for (i = 0; i < n; i++)
         {
           double x = i / (double)n;
           double r10 = gsl_cheb_eval_n (cs, 10, x);
           double r40 = gsl_cheb_eval (cs, x);
           printf ("%g %g %g %g\n", 
                   x, GSL_FN_EVAL (&F, x), r10, r40);
         }
     
       gsl_cheb_free (cs);
     
       return 0;
     }

The output from the program gives the original function, 10-th order
approximation and 40-th order approximation, all sampled at intervals of
0.001 in x.


File: gsl-ref.info,  Node: Chebyshev Approximation References and Further Reading,  Prev: Chebyshev Approximation examples,  Up: Chebyshev Approximations

References and Further Reading
==============================

The following paper describes the use of Chebyshev series,

     R. Broucke, "Ten Subroutines for the Manipulation of Chebyshev
     Series [C1] (Algorithm 446)". `Communications of the ACM' 16(4),
     254-256 (1973)


File: gsl-ref.info,  Node: Series Acceleration,  Next: Wavelet Transforms,  Prev: Chebyshev Approximations,  Up: Top

Series Acceleration
*******************

The functions described in this chapter accelerate the convergence of a
series using the Levin u-transform.  This method takes a small number of
terms from the start of a series and uses a systematic approximation to
compute an extrapolated value and an estimate of its error.  The
u-transform works for both convergent and divergent series, including
asymptotic series.

   These functions are declared in the header file `gsl_sum.h'.

* Menu:

* Acceleration functions::
* Acceleration functions without error estimation::
* Example of accelerating a series::
* Series Acceleration References::


File: gsl-ref.info,  Node: Acceleration functions,  Next: Acceleration functions without error estimation,  Up: Series Acceleration

Acceleration functions
======================

The following functions compute the full Levin u-transform of a series
with its error estimate.  The error estimate is computed by propagating
rounding errors from each term through to the final extrapolation.

   These functions are intended for summing analytic series where each
term is known to high accuracy, and the rounding errors are assumed to
originate from finite precision. They are taken to be relative errors of
order `GSL_DBL_EPSILON' for each term.

   The calculation of the error in the extrapolated value is an O(N^2)
process, which is expensive in time and memory.  A faster but less
reliable method which estimates the error from the convergence of the
extrapolated value is described in the next section.  For the method
described here a full table of intermediate values and derivatives
through to O(N) must be computed and stored, but this does give a
reliable error estimate.

 - Function: gsl_sum_levin_u_workspace * gsl_sum_levin_u_alloc (size_t
          N)
     This function allocates a workspace for a Levin u-transform of N
     terms.  The size of the workspace is O(2n^2 + 3n).

 - Function: int gsl_sum_levin_u_free (gsl_sum_levin_u_workspace * W)
     This function frees the memory associated with the workspace W.

 - Function: int gsl_sum_levin_u_accel (const double * ARRAY, size_t
          ARRAY_SIZE, gsl_sum_levin_u_workspace * W, double *
          SUM_ACCEL, double * ABSERR)
     This function takes the terms of a series in ARRAY of size
     ARRAY_SIZE and computes the extrapolated limit of the series using
     a Levin u-transform.  Additional working space must be provided in
     W.  The extrapolated sum is stored in SUM_ACCEL, with an estimate
     of the absolute error stored in ABSERR.  The actual term-by-term
     sum is returned in `w->sum_plain'. The algorithm calculates the
     truncation error (the difference between two successive
     extrapolations) and round-off error (propagated from the individual
     terms) to choose an optimal number of terms for the extrapolation.


File: gsl-ref.info,  Node: Acceleration functions without error estimation,  Next: Example of accelerating a series,  Prev: Acceleration functions,  Up: Series Acceleration

Acceleration functions without error estimation
===============================================

The functions described in this section compute the Levin u-transform of
series and attempt to estimate the error from the "truncation error" in
the extrapolation, the difference between the final two approximations.
Using this method avoids the need to compute an intermediate table of
derivatives because the error is estimated from the behavior of the
extrapolated value itself. Consequently this algorithm is an O(N)
process and only requires O(N) terms of storage.  If the series
converges sufficiently fast then this procedure can be acceptable.  It
is appropriate to use this method when there is a need to compute many
extrapolations of series with similar convergence properties at
high-speed.  For example, when numerically integrating a function
defined by a parameterized series where the parameter varies only
slightly. A reliable error estimate should be computed first using the
full algorithm described above in order to verify the consistency of the
results.

 - Function: gsl_sum_levin_utrunc_workspace *
          gsl_sum_levin_utrunc_alloc (size_t N)
     This function allocates a workspace for a Levin u-transform of N
     terms, without error estimation.  The size of the workspace is
     O(3n).

 - Function: int gsl_sum_levin_utrunc_free
          (gsl_sum_levin_utrunc_workspace * W)
     This function frees the memory associated with the workspace W.

 - Function: int gsl_sum_levin_utrunc_accel (const double * ARRAY,
          size_t ARRAY_SIZE, gsl_sum_levin_utrunc_workspace * W, double
          * SUM_ACCEL, double * ABSERR_TRUNC)
     This function takes the terms of a series in ARRAY of size
     ARRAY_SIZE and computes the extrapolated limit of the series using
     a Levin u-transform.  Additional working space must be provided in
     W.  The extrapolated sum is stored in SUM_ACCEL.  The actual
     term-by-term sum is returned in `w->sum_plain'. The algorithm
     terminates when the difference between two successive
     extrapolations reaches a minimum or is sufficiently small. The
     difference between these two values is used as estimate of the
     error and is stored in ABSERR_TRUNC.  To improve the reliability
     of the algorithm the extrapolated values are replaced by moving
     averages when calculating the truncation error, smoothing out any
     fluctuations.


File: gsl-ref.info,  Node: Example of accelerating a series,  Next: Series Acceleration References,  Prev: Acceleration functions without error estimation,  Up: Series Acceleration

Examples
========

The following code calculates an estimate of \zeta(2) = \pi^2 / 6 using
the series,

     \zeta(2) = 1 + 1/2^2 + 1/3^2 + 1/4^2 + ...

After N terms the error in the sum is O(1/N), making direct summation
of the series converge slowly.

     #include <stdio.h>
     #include <gsl/gsl_math.h>
     #include <gsl/gsl_sum.h>
     
     #define N 20
     
     int
     main (void)
     {
       double t[N];
       double sum_accel, err;
       double sum = 0;
       int n;
       
       gsl_sum_levin_u_workspace * w 
         = gsl_sum_levin_u_alloc (N);
     
       const double zeta_2 = M_PI * M_PI / 6.0;
       
       /* terms for zeta(2) = \sum_{n=1}^{\infty} 1/n^2 */
     
       for (n = 0; n < N; n++)
         {
           double np1 = n + 1.0;
           t[n] = 1.0 / (np1 * np1);
           sum += t[n];
         }
       
       gsl_sum_levin_u_accel (t, N, w, &sum_accel, &err);
     
       printf ("term-by-term sum = % .16f using %d terms\n", 
               sum, N);
     
       printf ("term-by-term sum = % .16f using %d terms\n", 
               w->sum_plain, w->terms_used);
     
       printf ("exact value      = % .16f\n", zeta_2);
       printf ("accelerated sum  = % .16f using %d terms\n", 
               sum_accel, w->terms_used);
     
       printf ("estimated error  = % .16f\n", err);
       printf ("actual error     = % .16f\n", 
               sum_accel - zeta_2);
     
       gsl_sum_levin_u_free (w);
       return 0;
     }

The output below shows that the Levin u-transform is able to obtain an
estimate of the sum to 1 part in 10^10 using the first eleven terms of
the series.  The error estimate returned by the function is also
accurate, giving the correct number of significant digits.

     bash$ ./a.out

     term-by-term sum =  1.5961632439130233 using 20 terms
     term-by-term sum =  1.5759958390005426 using 13 terms
     exact value      =  1.6449340668482264
     accelerated sum  =  1.6449340668166479 using 13 terms
     estimated error  =  0.0000000000508580
     actual error     = -0.0000000000315785

Note that a direct summation of this series would require 10^10 terms
to achieve the same precision as the accelerated sum does in 13 terms.


File: gsl-ref.info,  Node: Series Acceleration References,  Prev: Example of accelerating a series,  Up: Series Acceleration

References and Further Reading
==============================

The algorithms used by these functions are described in the following
papers,

     T. Fessler, W.F. Ford, D.A. Smith, HURRY: An acceleration
     algorithm for scalar sequences and series `ACM Transactions on
     Mathematical Software', 9(3):346-354, 1983.  and Algorithm 602
     9(3):355-357, 1983.

The theory of the u-transform was presented by Levin,

     D. Levin, Development of Non-Linear Transformations for Improving
     Convergence of Sequences, `Intern. J. Computer Math.' B3:371-388,
     1973.

A review paper on the Levin Transform is available online,
     Herbert H. H. Homeier, Scalar Levin-Type Sequence Transformations,
     `http://xxx.lanl.gov/abs/math/0005209'.


File: gsl-ref.info,  Node: Wavelet Transforms,  Next: Discrete Hankel Transforms,  Prev: Series Acceleration,  Up: Top

Wavelet Transforms
******************

This chapter describes functions for performing Discrete Wavelet
Transforms (DWTs).  The library includes wavelets for real data in both
one and two dimensions.  The wavelet functions are declared in the
header files `gsl_wavelet.h' and `gsl_wavelet2d.h'.

* Menu:

* DWT Definitions::
* DWT Initialization::
* DWT Transform Functions::
* DWT Examples::
* DWT References::


File: gsl-ref.info,  Node: DWT Definitions,  Next: DWT Initialization,  Up: Wavelet Transforms

Definitions
===========

The continuous wavelet transform and its inverse are defined by the
relations,

     w(s,\tau) = \int f(t) * \psi^*_{s,\tau}(t) dt

and,

     f(t) = \int \int_{-\infty}^\infty w(s, \tau) * \psi_{s,\tau}(t) d\tau ds

where the basis functions \psi_{s,\tau} are obtained by scaling and
translation from a single function, referred to as the "mother wavelet".

   The discrete version of the wavelet transform acts on evenly sampled
data, with fixed scaling and translation steps (s, \tau).  The
frequency and time axes are sampled "dyadically" on scales of 2^j
through a level parameter j.  The resulting family of functions
{\psi_{j,n}} constitutes an orthonormal basis for square-integrable
signals.

   The discrete wavelet transform is an O(N) algorithm, and is also
referred to as the "fast wavelet transform".


File: gsl-ref.info,  Node: DWT Initialization,  Next: DWT Transform Functions,  Prev: DWT Definitions,  Up: Wavelet Transforms

Initialization
==============

The `gsl_wavelet' structure contains the filter coefficients defining
the wavelet and associated offset parameters (for wavelets with
centered support).

 - Function: gsl_wavelet * gsl_wavelet_alloc (const gsl_wavelet_type *
          T, size_t K)
     This function allocates and initializes a wavelet object of type
     T.  The parameter K selects the specific member of the wavelet
     family.  A null pointer is returned if insufficient memory is
     available or if a unsupported member is selected.

   The following wavelet types are implemented:

 - Wavelet: gsl_wavelet_daubechies
 - Wavelet: gsl_wavelet_daubechies_centered
     The is the Daubechies wavelet family of maximum phase with k/2
     vanishing moments.  The implemented wavelets are k=4, 6, ..., 20,
     with K even.

 - Wavelet: gsl_wavelet_haar
 - Wavelet: gsl_wavelet_haar_centered
     This is the Haar wavelet.  The only valid choice of k for the Haar
     wavelet is k=2.

 - Wavelet: gsl_wavelet_bspline
 - Wavelet: gsl_wavelet_bspline_centered
     This is the biorthogonal B-spline wavelet family of order (i,j).
     The implemented values of k = 100*i + j are 103, 105, 202, 204,
     206, 208, 301, 303, 305 307, 309.

The centered forms of the wavelets align the coefficients of the various
sub-bands on edges.  Thus the resulting visualization of the
coefficients of the wavelet transform in the phase plane is easier to
understand.

 - Function: const char * gsl_wavelet_name (const gsl_wavelet * W)
     This function returns a pointer to the name of the wavelet family
     for W.

 - Function: void gsl_wavelet_free (gsl_wavelet * W)
     This function frees the wavelet object W.

   The `gsl_wavelet_workspace' structure contains scratch space of the
same size as the input data, for holding intermediate results during the
transform.

 - Function: gsl_wavelet_workspace * gsl_wavelet_workspace_alloc
          (size_t N)
     This function allocates a workspace for the discrete wavelet
     transform.  To perform a one-dimensional transform on N elements,
     a workspace of size N must be provided.  For two-dimensional
     transforms of N-by-N matrices it is sufficient to allocate a
     workspace of size N, since the transform operates on individual
     rows and columns.

 - Function: void gsl_wavelet_workspace_free (gsl_wavelet_workspace *
          WORKSPACE)
     This function frees the allocated workspace WORKSPACE.


File: gsl-ref.info,  Node: DWT Transform Functions,  Next: DWT Examples,  Prev: DWT Initialization,  Up: Wavelet Transforms

Transform Functions
===================

This sections describes the actual functions performing the discrete
wavelet transform.  Note that the transforms use periodic boundary
conditions.  If the signal is not periodic in the sample length then
spurious coefficients will appear at the beginning and end of each level
of the transform.

* Menu:

* DWT in one dimension::
* DWT in two dimension::


File: gsl-ref.info,  Node: DWT in one dimension,  Next: DWT in two dimension,  Up: DWT Transform Functions

Wavelet transforms in one dimension
-----------------------------------

 - Function: int gsl_dwt_transform (const gsl_wavelet * W, double *
          DATA, size_t STRIDE, size_t N, gsl_wavelet_direction DIR,
          gsl_wavelet_workspace * WORK)
 - Function: int gsl_dwt_transform_forward (const gsl_wavelet * W,
          double * DATA, size_t STRIDE, size_t N, gsl_wavelet_workspace
          * WORK)
 - Function: int gsl_dwt_transform_inverse (const gsl_wavelet * W,
          double * DATA, size_t STRIDE, size_t N, gsl_wavelet_workspace
          * WORK)
     These functions compute in-place forward and inverse discrete
     wavelet transforms of length N with stride STRIDE on the array
     DATA. The length of the transform N is restricted to powers of
     two.  For the `transform' version of the function the argument DIR
     can be either `forward' (+1) or `backward' (-1).  A workspace WORK
     of length N must be provided.

     For the forward transform, the elements of the original array are
     replaced by the discrete wavelet transform f_i -> w_{j,k} in a
     packed triangular storage layout, where J is the index of the level
     j = 0 ... J-1 and K is the index of the coefficient within each
     level, k = 0 ... (2^j)-1.  The total number of levels is J =
     \log_2(n).  The output data has the following form,

          (s_{-1,0}, d_{0,0}, d_{1,0}, d_{1,1}, d_{2,0}, ...,
            d_{j,k}, ..., d_{J-1,2^{J-1}-1})

     where the first element is the smoothing coefficient s_-1,0,
     followed by the detail coefficients d_j,k for each level j.  The
     backward transform inverts these coefficients to obtain the
     original data.

     These functions return a status of `GSL_SUCCESS' upon successful
     completion.  `GSL_EINVAL' is returned if N is not an integer power
     of 2 or if insufficient workspace is provided.


File: gsl-ref.info,  Node: DWT in two dimension,  Prev: DWT in one dimension,  Up: DWT Transform Functions

Wavelet transforms in two dimension
-----------------------------------

The library provides functions to perform two-dimensional discrete
wavelet transforms on square matrices.  The matrix dimensions must be an
integer power of two.  There are two possible orderings of the rows and
columns in the two-dimensional wavelet transform, referred to as the
"standard" and "non-standard" forms.

   The "standard" transform performs a discrete wavelet transform on all
rows of the matrix, followed by a separate discrete wavelet transform on
the columns of the resulting row-transformed matrix.  This procedure
uses the same ordering as a two-dimensional fourier transform.

   The "non-standard" transform is performed in interleaved passes of
the each level of the transform on the rows and columns of the matrix.
The first level of the transform is carried out on the matrix rows, and
then the columns of the partially row-transformed data.  This procedure
is then repeated across the rows and columns of the data for the
subsequent levels of the transform, until the full discrete wavelet
transform is complete.  The non-standard form of the discrete wavelet
transform is typically used in image analysis.

   The functions described in this section are declared in the header
file `gsl_wavelet2d.h'.

 - Function: int gsl_wavelet2d_transform (const gsl_wavelet * W, double
          * DATA, size_t TDA, size_t SIZE1, size_t SIZE2,
          gsl_wavelet_direction DIR, gsl_wavelet_workspace * WORK)
 - Function: int gsl_wavelet2d_transform_forward (const gsl_wavelet *
          W, double * DATA, size_t TDA, size_t SIZE1, size_t SIZE2,
          gsl_wavelet_workspace * WORK)
 - Function: int gsl_wavelet2d_transform_inverse (const gsl_wavelet *
          W, double * DATA, size_t TDA, size_t SIZE1, size_t SIZE2,
          gsl_wavelet_workspace * WORK)
     These functions compute two-dimensional in-place forward and
     inverse discrete wavelet transforms in standard and non-standard
     forms on the array DATA stored in row-major form with dimensions
     SIZE1 and SIZE2 and physical row length TDA.  The dimensions must
     be equal (square matrix) and are restricted to powers of two.  For
     the `transform' version of the function the argument DIR can be
     either `forward' (+1) or `backward' (-1).  A workspace WORK of the
     appropriate size must be provided.  On exit, the appropriate
     elements of the array DATA are replaced by their two-dimensional
     wavelet transform.

     The functions return a status of `GSL_SUCCESS' upon successful
     completion.  `GSL_EINVAL' is returned if SIZE1 and SIZE2 are not
     equal and integer powers of 2, or if insufficient workspace is
     provided.

 - Function: int gsl_wavelet2d_transform_matrix (const gsl_wavelet * W,
          gsl_matrix * M, gsl_wavelet_direction DIR,
          gsl_wavelet_workspace * WORK)
 - Function: int gsl_wavelet2d_transform_matrix_forward (const
          gsl_wavelet * W, gsl_matrix * M, gsl_wavelet_workspace * WORK)
 - Function: int gsl_wavelet2d_transform_matrix_inverse (const
          gsl_wavelet * W, gsl_matrix * M, gsl_wavelet_workspace * WORK)
     These functions compute the two-dimensional in-place wavelet
     transform on a matrix A.

 - Function: int gsl_wavelet2d_nstransform (const gsl_wavelet * W,
          double * DATA, size_t TDA, size_t SIZE1, size_t SIZE2,
          gsl_wavelet_direction DIR, gsl_wavelet_workspace * WORK)
 - Function: int gsl_wavelet2d_nstransform_forward (const gsl_wavelet *
          W, double * DATA, size_t TDA, size_t SIZE1, size_t SIZE2,
          gsl_wavelet_workspace * WORK)
 - Function: int gsl_wavelet2d_nstransform_inverse (const gsl_wavelet *
          W, double * DATA, size_t TDA, size_t SIZE1, size_t SIZE2,
          gsl_wavelet_workspace * WORK)
     These functions compute the two-dimensional wavelet transform in
     non-standard form.

 - Function: int gsl_wavelet2d_nstransform_matrix (const gsl_wavelet *
          W, gsl_matrix * M, gsl_wavelet_direction DIR,
          gsl_wavelet_workspace * WORK)
 - Function: int gsl_wavelet2d_nstransform_matrix_forward (const
          gsl_wavelet * W, gsl_matrix * M, gsl_wavelet_workspace * WORK)
 - Function: int gsl_wavelet2d_nstransform_matrix_inverse (const
          gsl_wavelet * W, gsl_matrix * M, gsl_wavelet_workspace * WORK)
     These functions compute the non-standard form of the
     two-dimensional in-place wavelet transform on a matrix A.


File: gsl-ref.info,  Node: DWT Examples,  Next: DWT References,  Prev: DWT Transform Functions,  Up: Wavelet Transforms

Example
=======

The following program demonstrates the use of the one-dimensional
wavelet transform functions.  It computes an approximation to an input
signal (of length 256) using the 20 largest components of the wavelet
transform, while setting the others to zero.

     #include <stdio.h>
     #include <math.h>
     #include <gsl/gsl_sort.h>
     #include <gsl/gsl_wavelet.h>
     
     int
     main (int argc, char **argv)
     {
       int i, n = 256, nc = 20;
       double *data = malloc (n * sizeof (double));
       double *abscoeff = malloc (n * sizeof (double));
       size_t *p = malloc (n * sizeof (size_t));
     
       FILE *f = fopen (argv[1], "r");
       for (i = 0; i < n; i++)
         {
           fscanf (f, "%lg", &data[i]);
         }
       fclose (f);
     
       {
         gsl_wavelet *w = gsl_wavelet_alloc (gsl_wavelet_daubechies, 4);
         gsl_wavelet_workspace *work = gsl_wavelet_workspace_alloc (n);
     
         gsl_wavelet_transform_forward (w, data, 1, n, work);
     
         for (i = 0; i < n; i++)
           {
             abscoeff[i] = fabs (data[i]);
           }
     
         gsl_sort_index (p, abscoeff, 1, n);
     
         for (i = 0; (i + nc) < n; i++)
           data[p[i]] = 0;
     
         gsl_wavelet_transform_inverse (w, data, 1, n, work);
       }
     
       for (i = 0; i < n; i++)
         {
           printf ("%g\n", data[i]);
         }
     }

The output can be used with the GNU plotutils `graph' program,

     $ ./a.out ecg.dat > dwt.dat
     $ graph -T ps -x 0 256 32 -h 0.3 -a dwt.dat > dwt.ps


File: gsl-ref.info,  Node: DWT References,  Prev: DWT Examples,  Up: Wavelet Transforms

References and Further Reading
==============================

The mathematical background to wavelet transforms is covered in the
original lectures by Daubechies,

     Ingrid Daubechies.  Ten Lectures on Wavelets.  `CBMS-NSF Regional
     Conference Series in Applied Mathematics' (1992), SIAM, ISBN
     0898712742.

An easy to read introduction to the subject with an emphasis on the
application of the wavelet transform in various branches of science is,

     Paul S. Addison. `The Illustrated Wavelet Transform Handbook'.
     Institute of Physics Publishing (2002), ISBN 0750306920.

For extensive coverage of signal analysis by wavelets, wavelet packets
and local cosine bases see,

     Ste'phane Mallat.  `A wavelet tour of signal processing' (Second
     edition). Academic Press (1999), ISBN 012466606X.

The concept of multiresolution analysis underlying the wavelet transform
is described in,

     S. G. Mallat.  Multiresolution Approximations and Wavelet
     Orthonormal Bases of L^2(R).  `Transactions of the American
     Mathematical Society', 315(1), 1989, 69-87.

     S. G. Mallat.  A Theory for Multiresolution Signal Decomposition
     -- The Wavelet Representation.  `IEEE Transactions on Pattern
     Analysis and Machine Intelligence', 11, 1989, 674-693.

The coefficients for the individual wavelet families implemented by the
library can be found in the following papers,

     I. Daubechies.  Orthonormal Bases of Compactly Supported Wavelets.
     `Communications on Pure and Applied Mathematics', 41 (1988)
     909-996.

     A. Cohen, I. Daubechies, and J.-C. Feauveau.  Biorthogonal Bases
     of Compactly Supported Wavelets.  `Communications on Pure and
     Applied Mathematics', 45 (1992) 485-560.

The PhysioNet archive of physiological datasets can be found online at
`http://www.physionet.org/' and is described in the following paper,

     Goldberger et al.  PhysioBank, PhysioToolkit, and PhysioNet:
     Components of a New Research Resource for Complex Physiologic
     Signals.  `Circulation' 101(23):e215-e220 2000.


File: gsl-ref.info,  Node: Discrete Hankel Transforms,  Next: One dimensional Root-Finding,  Prev: Wavelet Transforms,  Up: Top

Discrete Hankel Transforms
**************************

This chapter describes functions for performing Discrete Hankel
Transforms (DHTs).  The functions are declared in the header file
`gsl_dht.h'.

* Menu:

* Discrete Hankel Transform Definition::
* Discrete Hankel Transform Functions::
* Discrete Hankel Transform References::


File: gsl-ref.info,  Node: Discrete Hankel Transform Definition,  Next: Discrete Hankel Transform Functions,  Up: Discrete Hankel Transforms

Definitions
===========

The discrete Hankel transform acts on a vector of sampled data, where
the samples are assumed to have been taken at points related to the
zeroes of a Bessel function of fixed order; compare this to the case of
the discrete Fourier transform, where samples are taken at points
related to the zeroes of the sine or cosine function.

   Specifically, let f(t) be a function on the unit interval.  Then the
finite \nu-Hankel transform of f(t) is defined to be the set of numbers
g_m given by,

so that,

Suppose that f is band-limited in the sense that g_m=0 for m > M. Then
we have the following fundamental sampling theorem.

It is this discrete expression which defines the discrete Hankel
transform. The kernel in the summation above defines the matrix of the
\nu-Hankel transform of size M-1.  The coefficients of this matrix,
being dependent on \nu and M, must be precomputed and stored; the
`gsl_dht' object encapsulates this data.  The allocation function
`gsl_dht_alloc' returns a `gsl_dht' object which must be properly
initialized with `gsl_dht_init' before it can be used to perform
transforms on data sample vectors, for fixed \nu and M, using the
`gsl_dht_apply' function. The implementation allows a scaling of the
fundamental interval, for convenience, so that one can assume the
function is defined on the interval [0,X], rather than the unit
interval.

   Notice that by assumption f(t) vanishes at the endpoints of the
interval, consistent with the inversion formula and the sampling
formula given above. Therefore, this transform corresponds to an
orthogonal expansion in eigenfunctions of the Dirichlet problem for the
Bessel differential equation.


File: gsl-ref.info,  Node: Discrete Hankel Transform Functions,  Next: Discrete Hankel Transform References,  Prev: Discrete Hankel Transform Definition,  Up: Discrete Hankel Transforms

Functions
=========

 - Function: gsl_dht * gsl_dht_alloc (size_t SIZE)
     This function allocates a Discrete Hankel transform object of size
     SIZE.

 - Function: int gsl_dht_init (gsl_dht * T, double NU, double XMAX)
     This function initializes the transform T for the given values of
     NU and X.

 - Function: gsl_dht * gsl_dht_new (size_t SIZE, double NU, double XMAX)
     This function allocates a Discrete Hankel transform object of size
     SIZE and initializes it for the given values of NU and X.

 - Function: void gsl_dht_free (gsl_dht * T)
     This function frees the transform T.

 - Function: int gsl_dht_apply (const gsl_dht * T, double * F_IN,
          double * F_OUT)
     This function applies the transform T to the array F_IN whose size
     is equal to the size of the transform.  The result is stored in
     the array F_OUT which must be of the same length.

 - Function: double gsl_dht_x_sample (const gsl_dht * T, int N)
     This function returns the value of the n'th sample point in the
     unit interval, (j_{\nu,n+1}/j_{\nu,M}) X. These are the points
     where the function f(t) is assumed to be sampled.

 - Function: double gsl_dht_k_sample (const gsl_dht * T, int N)
     This function returns the value of the n'th sample point in
     "k-space", j_{\nu,n+1}/X.


File: gsl-ref.info,  Node: Discrete Hankel Transform References,  Prev: Discrete Hankel Transform Functions,  Up: Discrete Hankel Transforms

References and Further Reading
==============================

The algorithms used by these functions are described in the following
papers,

     H. Fisk Johnson, Comp. Phys. Comm. 43, 181 (1987).

     D. Lemoine, J. Chem. Phys. 101, 3936 (1994).


File: gsl-ref.info,  Node: One dimensional Root-Finding,  Next: One dimensional Minimization,  Prev: Discrete Hankel Transforms,  Up: Top

One dimensional Root-Finding
****************************

This chapter describes routines for finding roots of arbitrary
one-dimensional functions.  The library provides low level components
for a variety of iterative solvers and convergence tests.  These can be
combined by the user to achieve the desired solution, with full access
to the intermediate steps of the iteration.  Each class of methods uses
the same framework, so that you can switch between solvers at runtime
without needing to recompile your program.  Each instance of a solver
keeps track of its own state, allowing the solvers to be used in
multi-threaded programs.

   The header file `gsl_roots.h' contains prototypes for the root
finding functions and related declarations.

* Menu:

* Root Finding Overview::
* Root Finding Caveats::
* Initializing the Solver::
* Providing the function to solve::
* Search Bounds and Guesses::
* Root Finding Iteration::
* Search Stopping Parameters::
* Root Bracketing Algorithms::
* Root Finding Algorithms using Derivatives::
* Root Finding Examples::
* Root Finding References and Further Reading::


File: gsl-ref.info,  Node: Root Finding Overview,  Next: Root Finding Caveats,  Up: One dimensional Root-Finding

Overview
========

One-dimensional root finding algorithms can be divided into two classes,
"root bracketing" and "root polishing".  Algorithms which proceed by
bracketing a root are guaranteed to converge.  Bracketing algorithms
begin with a bounded region known to contain a root.  The size of this
bounded region is reduced, iteratively, until it encloses the root to a
desired tolerance.  This provides a rigorous error estimate for the
location of the root.

   The technique of "root polishing" attempts to improve an initial
guess to the root.  These algorithms converge only if started "close
enough" to a root, and sacrifice a rigorous error bound for speed.  By
approximating the behavior of a function in the vicinity of a root they
attempt to find a higher order improvement of an initial guess.  When
the behavior of the function is compatible with the algorithm and a good
initial guess is available a polishing algorithm can provide rapid
convergence.

   In GSL both types of algorithm are available in similar frameworks.
The user provides a high-level driver for the algorithms, and the
library provides the individual functions necessary for each of the
steps.  There are three main phases of the iteration.  The steps are,

   * initialize solver state, S, for algorithm T

   * update S using the iteration T

   * test S for convergence, and repeat iteration if necessary

The state for bracketing solvers is held in a `gsl_root_fsolver'
struct.  The updating procedure uses only function evaluations (not
derivatives).  The state for root polishing solvers is held in a
`gsl_root_fdfsolver' struct.  The updates require both the function and
its derivative (hence the name `fdf') to be supplied by the user.


File: gsl-ref.info,  Node: Root Finding Caveats,  Next: Initializing the Solver,  Prev: Root Finding Overview,  Up: One dimensional Root-Finding

Caveats
=======

Note that root finding functions can only search for one root at a time.
When there are several roots in the search area, the first root to be
found will be returned; however it is difficult to predict which of the
roots this will be. _In most cases, no error will be reported if you
try to find a root in an area where there is more than one._

   Care must be taken when a function may have a multiple root (such as
f(x) = (x-x_0)^2 or f(x) = (x-x_0)^3).  It is not possible to use
root-bracketing algorithms on even-multiplicity roots.  For these
algorithms the initial interval must contain a zero-crossing, where the
function is negative at one end of the interval and positive at the
other end.  Roots with even-multiplicity do not cross zero, but only
touch it instantaneously.  Algorithms based on root bracketing will
still work for odd-multiplicity roots (e.g. cubic, quintic, ...).  Root
polishing algorithms generally work with higher multiplicity roots, but
at reduced rate of convergence.  In these cases the "Steffenson
algorithm" can be used to accelerate the convergence of multiple roots.

   While it is not absolutely required that f have a root within the
search region, numerical root finding functions should not be used
haphazardly to check for the _existence_ of roots.  There are better
ways to do this.  Because it is easy to create situations where
numerical root finders go awry, it is a bad idea to throw a root finder
at a function you do not know much about.  In general it is best to
examine the function visually by plotting before searching for a root.


File: gsl-ref.info,  Node: Initializing the Solver,  Next: Providing the function to solve,  Prev: Root Finding Caveats,  Up: One dimensional Root-Finding

Initializing the Solver
=======================

 - Function: gsl_root_fsolver * gsl_root_fsolver_alloc (const
          gsl_root_fsolver_type * T)
     This function returns a pointer to a newly allocated instance of a
     solver of type T.  For example, the following code creates an
     instance of a bisection solver,

          const gsl_root_fsolver_type * T
            = gsl_root_fsolver_bisection;
          gsl_root_fsolver * s
            = gsl_root_fsolver_alloc (T);

     If there is insufficient memory to create the solver then the
     function returns a null pointer and the error handler is invoked
     with an error code of `GSL_ENOMEM'.

 - Function: gsl_root_fdfsolver * gsl_root_fdfsolver_alloc (const
          gsl_root_fdfsolver_type * T)
     This function returns a pointer to a newly allocated instance of a
     derivative-based solver of type T.  For example, the following
     code creates an instance of a Newton-Raphson solver,

          const gsl_root_fdfsolver_type * T
            = gsl_root_fdfsolver_newton;
          gsl_root_fdfsolver * s
            = gsl_root_fdfsolver_alloc (T);

     If there is insufficient memory to create the solver then the
     function returns a null pointer and the error handler is invoked
     with an error code of `GSL_ENOMEM'.

 - Function: int gsl_root_fsolver_set (gsl_root_fsolver * S,
          gsl_function * F, double X_LOWER, double X_UPPER)
     This function initializes, or reinitializes, an existing solver S
     to use the function F and the initial search interval [X_LOWER,
     X_UPPER].

 - Function: int gsl_root_fdfsolver_set (gsl_root_fdfsolver * S,
          gsl_function_fdf * FDF, double ROOT)
     This function initializes, or reinitializes, an existing solver S
     to use the function and derivative FDF and the initial guess ROOT.

 - Function: void gsl_root_fsolver_free (gsl_root_fsolver * S)
 - Function: void gsl_root_fdfsolver_free (gsl_root_fdfsolver * S)
     These functions free all the memory associated with the solver S.

 - Function: const char * gsl_root_fsolver_name (const gsl_root_fsolver
          * S)
 - Function: const char * gsl_root_fdfsolver_name (const
          gsl_root_fdfsolver * S)
     These functions return a pointer to the name of the solver.  For
     example,

          printf ("s is a '%s' solver\n",
                  gsl_root_fsolver_name (s));

     would print something like `s is a 'bisection' solver'.


File: gsl-ref.info,  Node: Providing the function to solve,  Next: Search Bounds and Guesses,  Prev: Initializing the Solver,  Up: One dimensional Root-Finding

Providing the function to solve
===============================

You must provide a continuous function of one variable for the root
finders to operate on, and, sometimes, its first derivative.  In order
to allow for general parameters the functions are defined by the
following data types:

 - Data Type: gsl_function
     This data type defines a general function with parameters.

    `double (* function) (double X, void * PARAMS)'
          this function should return the value f(x,params) for
          argument X and parameters PARAMS

    `void * params'
          a pointer to the parameters of the function

   Here is an example for the general quadratic function,

     f(x) = a x^2 + b x + c

with a = 3, b = 2, c = 1.  The following code defines a `gsl_function'
`F' which you could pass to a root finder:

     struct my_f_params { double a; double b; double c; };
     
     double
     my_f (double x, void * p) {
        struct my_f_params * params
          = (struct my_f_params *)p;
        double a = (params->a);
        double b = (params->b);
        double c = (params->c);
     
        return  (a * x + b) * x + c;
     }
     
     gsl_function F;
     struct my_f_params params = { 3.0, 2.0, 1.0 };
     
     F.function = &my_f;
     F.params = &params;

The function f(x) can be evaluated using the following macro,

     #define GSL_FN_EVAL(F,x)
         (*((F)->function))(x,(F)->params)

 - Data Type: gsl_function_fdf
     This data type defines a general function with parameters and its
     first derivative.

    `double (* f) (double X, void * PARAMS)'
          this function should return the value of f(x,params) for
          argument X and parameters PARAMS

    `double (* df) (double X, void * PARAMS)'
          this function should return the value of the derivative of F
          with respect to X, f'(x,params), for argument X and
          parameters PARAMS

    `void (* fdf) (double X, void * PARAMS, double * F, double * Df)'
          this function should set the values of the function F to
          f(x,params) and its derivative DF to f'(x,params) for
          argument X and parameters PARAMS.  This function provides an
          optimization of the separate functions for f(x) and f'(x)--it
          is always faster to compute the function and its derivative
          at the same time.

    `void * params'
          a pointer to the parameters of the function

   Here is an example where f(x) = 2\exp(2x):

     double
     my_f (double x, void * params)
     {
        return exp (2 * x);
     }
     
     double
     my_df (double x, void * params)
     {
        return 2 * exp (2 * x);
     }
     
     void
     my_fdf (double x, void * params,
             double * f, double * df)
     {
        double t = exp (2 * x);
     
        *f = t;
        *df = 2 * t;   /* uses existing value */
     }
     
     gsl_function_fdf FDF;
     
     FDF.f = &my_f;
     FDF.df = &my_df;
     FDF.fdf = &my_fdf;
     FDF.params = 0;

The function f(x) can be evaluated using the following macro,

     #define GSL_FN_FDF_EVAL_F(FDF,x)
          (*((FDF)->f))(x,(FDF)->params)

The derivative f'(x) can be evaluated using the following macro,

     #define GSL_FN_FDF_EVAL_DF(FDF,x)
          (*((FDF)->df))(x,(FDF)->params)

and both the function y = f(x) and its derivative dy = f'(x) can be
evaluated at the same time using the following macro,

     #define GSL_FN_FDF_EVAL_F_DF(FDF,x,y,dy)
          (*((FDF)->fdf))(x,(FDF)->params,(y),(dy))

The macro stores f(x) in its Y argument and f'(x) in its DY
argument--both of these should be pointers to `double'.


File: gsl-ref.info,  Node: Search Bounds and Guesses,  Next: Root Finding Iteration,  Prev: Providing the function to solve,  Up: One dimensional Root-Finding

Search Bounds and Guesses
=========================

You provide either search bounds or an initial guess; this section
explains how search bounds and guesses work and how function arguments
control them.

   A guess is simply an x value which is iterated until it is within
the desired precision of a root.  It takes the form of a `double'.

   Search bounds are the endpoints of a interval which is iterated until
the length of the interval is smaller than the requested precision.  The
interval is defined by two values, the lower limit and the upper limit.
Whether the endpoints are intended to be included in the interval or not
depends on the context in which the interval is used.


File: gsl-ref.info,  Node: Root Finding Iteration,  Next: Search Stopping Parameters,  Prev: Search Bounds and Guesses,  Up: One dimensional Root-Finding

Iteration
=========

The following functions drive the iteration of each algorithm.  Each
function performs one iteration to update the state of any solver of the
corresponding type.  The same functions work for all solvers so that
different methods can be substituted at runtime without modifications to
the code.

 - Function: int gsl_root_fsolver_iterate (gsl_root_fsolver * S)
 - Function: int gsl_root_fdfsolver_iterate (gsl_root_fdfsolver * S)
     These functions perform a single iteration of the solver S.  If the
     iteration encounters an unexpected problem then an error code will
     be returned,

    `GSL_EBADFUNC'
          the iteration encountered a singular point where the function
          or its derivative evaluated to `Inf' or `NaN'.

    `GSL_EZERODIV'
          the derivative of the function vanished at the iteration
          point, preventing the algorithm from continuing without a
          division by zero.

   The solver maintains a current best estimate of the root at all
times.  The bracketing solvers also keep track of the current best
interval bounding the root.  This information can be accessed with the
following auxiliary functions,

 - Function: double gsl_root_fsolver_root (const gsl_root_fsolver * S)
 - Function: double gsl_root_fdfsolver_root (const gsl_root_fdfsolver *
          S)
     These functions return the current estimate of the root for the
     solver S.

 - Function: double gsl_root_fsolver_x_lower (const gsl_root_fsolver *
          S)
 - Function: double gsl_root_fsolver_x_upper (const gsl_root_fsolver *
          S)
     These functions return the current bracketing interval for the
     solver S.


File: gsl-ref.info,  Node: Search Stopping Parameters,  Next: Root Bracketing Algorithms,  Prev: Root Finding Iteration,  Up: One dimensional Root-Finding

Search Stopping Parameters
==========================

A root finding procedure should stop when one of the following
conditions is true:

   * A root has been found to within the user-specified precision.

   * A user-specified maximum number of iterations has been reached.

   * An error has occurred.

The handling of these conditions is under user control.  The functions
below allow the user to test the precision of the current result in
several standard ways.

 - Function: int gsl_root_test_interval (double X_LOWER, double
          X_UPPER, double EPSABS, double EPSREL)
     This function tests for the convergence of the interval [X_LOWER,
     X_UPPER] with absolute error EPSABS and relative error EPSREL.
     The test returns `GSL_SUCCESS' if the following condition is
     achieved,

          |a - b| < epsabs + epsrel min(|a|,|b|)

     when the interval x = [a,b] does not include the origin.  If the
     interval includes the origin then \min(|a|,|b|) is replaced by
     zero (which is the minimum value of |x| over the interval).  This
     ensures that the relative error is accurately estimated for roots
     close to the origin.

     This condition on the interval also implies that any estimate of
     the root r in the interval satisfies the same condition with
     respect to the true root r^*,

          |r - r^*| < epsabs + epsrel r^*

     assuming that the true root r^* is contained within the interval.

 - Function: int gsl_root_test_delta (double X1, double X0, double
          EPSABS, double EPSREL)
     This function tests for the convergence of the sequence ..., X0,
     X1 with absolute error EPSABS and relative error EPSREL.  The test
     returns `GSL_SUCCESS' if the following condition is achieved,

          |x_1 - x_0| < epsabs + epsrel |x_1|

     and returns `GSL_CONTINUE' otherwise.

 - Function: int gsl_root_test_residual (double F, double EPSABS)
     This function tests the residual value F against the absolute
     error bound EPSABS.  The test returns `GSL_SUCCESS' if the
     following condition is achieved,

          |f| < epsabs

     and returns `GSL_CONTINUE' otherwise.  This criterion is suitable
     for situations where the precise location of the root, x, is
     unimportant provided a value can be found where the residual,
     |f(x)|, is small enough.


File: gsl-ref.info,  Node: Root Bracketing Algorithms,  Next: Root Finding Algorithms using Derivatives,  Prev: Search Stopping Parameters,  Up: One dimensional Root-Finding

Root Bracketing Algorithms
==========================

The root bracketing algorithms described in this section require an
initial interval which is guaranteed to contain a root--if a and b are
the endpoints of the interval then f(a) must differ in sign from f(b).
This ensures that the function crosses zero at least once in the
interval.  If a valid initial interval is used then these algorithm
cannot fail, provided the function is well-behaved.

   Note that a bracketing algorithm cannot find roots of even degree,
since these do not cross the x-axis.

 - Solver: gsl_root_fsolver_bisection
     The "bisection algorithm" is the simplest method of bracketing the
     roots of a function.   It is the slowest algorithm provided by the
     library, with linear convergence.

     On each iteration, the interval is bisected and the value of the
     function at the midpoint is calculated.  The sign of this value is
     used to determine which half of the interval does not contain a
     root.  That half is discarded to give a new, smaller interval
     containing the root.  This procedure can be continued indefinitely
     until the interval is sufficiently small.

     At any time the current estimate of the root is taken as the
     midpoint of the interval.


 - Solver: gsl_root_fsolver_falsepos
     The "false position algorithm" is a method of finding roots based
     on linear interpolation.  Its convergence is linear, but it is
     usually faster than bisection.

     On each iteration a line is drawn between the endpoints (a,f(a))
     and (b,f(b)) and the point where this line crosses the x-axis
     taken as a "midpoint".  The value of the function at this point is
     calculated and its sign is used to determine which side of the
     interval does not contain a root.  That side is discarded to give a
     new, smaller interval containing the root.  This procedure can be
     continued indefinitely until the interval is sufficiently small.

     The best estimate of the root is taken from the linear
     interpolation of the interval on the current iteration.


 - Solver: gsl_root_fsolver_brent
     The "Brent-Dekker method" (referred to here as "Brent's method")
     combines an interpolation strategy with the bisection algorithm.
     This produces a fast algorithm which is still robust.

     On each iteration Brent's method approximates the function using an
     interpolating curve.  On the first iteration this is a linear
     interpolation of the two endpoints.  For subsequent iterations the
     algorithm uses an inverse quadratic fit to the last three points,
     for higher accuracy.  The intercept of the interpolating curve
     with the x-axis is taken as a guess for the root.  If it lies
     within the bounds of the current interval then the interpolating
     point is accepted, and used to generate a smaller interval.  If
     the interpolating point is not accepted then the algorithm falls
     back to an ordinary bisection step.

     The best estimate of the root is taken from the most recent
     interpolation or bisection.


File: gsl-ref.info,  Node: Root Finding Algorithms using Derivatives,  Next: Root Finding Examples,  Prev: Root Bracketing Algorithms,  Up: One dimensional Root-Finding

Root Finding Algorithms using Derivatives
=========================================

The root polishing algorithms described in this section require an
initial guess for the location of the root.  There is no absolute
guarantee of convergence--the function must be suitable for this
technique and the initial guess must be sufficiently close to the root
for it to work.  When these conditions are satisfied then convergence is
quadratic.

   These algorithms make use of both the function and its derivative.

 - Derivative Solver: gsl_root_fdfsolver_newton
     Newton's Method is the standard root-polishing algorithm.  The
     algorithm begins with an initial guess for the location of the
     root.  On each iteration, a line tangent to the function f is
     drawn at that position.  The point where this line crosses the
     x-axis becomes the new guess.  The iteration is defined by the
     following sequence,

          x_{i+1} = x_i - f(x_i)/f'(x_i)

     Newton's method converges quadratically for single roots, and
     linearly for multiple roots.


 - Derivative Solver: gsl_root_fdfsolver_secant
     The "secant method" is a simplified version of Newton's method
     which does not require the computation of the derivative on every
     step.

     On its first iteration the algorithm begins with Newton's method,
     using the derivative to compute a first step,

          x_1 = x_0 - f(x_0)/f'(x_0)

     Subsequent iterations avoid the evaluation of the derivative by
     replacing it with a numerical estimate, the slope of the line
     through the previous two points,

          x_{i+1} = x_i f(x_i) / f'_{est} where
           f'_{est} = (f(x_i) - f(x_{i-1})/(x_i - x_{i-1})

     When the derivative does not change significantly in the vicinity
     of the root the secant method gives a useful saving.
     Asymptotically the secant method is faster than Newton's method
     whenever the cost of evaluating the derivative is more than 0.44
     times the cost of evaluating the function itself.  As with all
     methods of computing a numerical derivative the estimate can
     suffer from cancellation errors if the separation of the points
     becomes too small.

     On single roots, the method has a convergence of order (1 + \sqrt
     5)/2 (approximately 1.62).  It converges linearly for multiple
     roots.


 - Derivative Solver: gsl_root_fdfsolver_steffenson
     The "Steffenson Method" provides the fastest convergence of all the
     routines.  It combines the basic Newton algorithm with an Aitken
     "delta-squared" acceleration.  If the Newton iterates are x_i then
     the acceleration procedure generates a new sequence R_i,

          R_i = x_i - (x_{i+1} - x_i)^2 / (x_{i+2} - 2 x_{i+1} + x_{i})

     which converges faster than the original sequence under reasonable
     conditions.  The new sequence requires three terms before it can
     produce its first value so the method returns accelerated values
     on the second and subsequent iterations.  On the first iteration
     it returns the ordinary Newton estimate.  The Newton iterate is
     also returned if the denominator of the acceleration term ever
     becomes zero.

     As with all acceleration procedures this method can become
     unstable if the function is not well-behaved.


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

Examples
========

For any root finding algorithm we need to prepare the function to be
solved.  For this example we will use the general quadratic equation
described earlier.  We first need a header file (`demo_fn.h') to define
the function parameters,

     struct quadratic_params
       {
         double a, b, c;
       };
     
     double quadratic (double x, void *params);
     double quadratic_deriv (double x, void *params);
     void quadratic_fdf (double x, void *params, 
                         double *y, double *dy);

We place the function definitions in a separate file (`demo_fn.c'),

     double
     quadratic (double x, void *params)
     {
       struct quadratic_params *p 
         = (struct quadratic_params *) params;
     
       double a = p->a;
       double b = p->b;
       double c = p->c;
     
       return (a * x + b) * x + c;
     }
     
     double
     quadratic_deriv (double x, void *params)
     {
       struct quadratic_params *p 
         = (struct quadratic_params *) params;
     
       double a = p->a;
       double b = p->b;
       double c = p->c;
     
       return 2.0 * a * x + b;
     }
     
     void
     quadratic_fdf (double x, void *params, 
                    double *y, double *dy)
     {
       struct quadratic_params *p 
         = (struct quadratic_params *) params;
     
       double a = p->a;
       double b = p->b;
       double c = p->c;
     
       *y = (a * x + b) * x + c;
       *dy = 2.0 * a * x + b;
     }

The first program uses the function solver `gsl_root_fsolver_brent' for
Brent's method and the general quadratic defined above to solve the
following equation,

     x^2 - 5 = 0

with solution x = \sqrt 5 = 2.236068...

     #include <stdio.h>
     #include <gsl/gsl_errno.h>
     #include <gsl/gsl_math.h>
     #include <gsl/gsl_roots.h>
     
     #include "demo_fn.h"
     #include "demo_fn.c"
     
     int
     main (void)
     {
       int status;
       int iter = 0, max_iter = 100;
       const gsl_root_fsolver_type *T;
       gsl_root_fsolver *s;
       double r = 0, r_expected = sqrt (5.0);
       double x_lo = 0.0, x_hi = 5.0;
       gsl_function F;
       struct quadratic_params params = {1.0, 0.0, -5.0};
     
       F.function = &quadratic;
       F.params = &params;
     
       T = gsl_root_fsolver_brent;
       s = gsl_root_fsolver_alloc (T);
       gsl_root_fsolver_set (s, &F, x_lo, x_hi);
     
       printf ("using %s method\n", 
               gsl_root_fsolver_name (s));
     
       printf ("%5s [%9s, %9s] %9s %10s %9s\n",
               "iter", "lower", "upper", "root", 
               "err", "err(est)");
     
       do
         {
           iter++;
           status = gsl_root_fsolver_iterate (s);
           r = gsl_root_fsolver_root (s);
           x_lo = gsl_root_fsolver_x_lower (s);
           x_hi = gsl_root_fsolver_x_upper (s);
           status = gsl_root_test_interval (x_lo, x_hi,
                                            0, 0.001);
     
           if (status == GSL_SUCCESS)
             printf ("Converged:\n");
     
           printf ("%5d [%.7f, %.7f] %.7f %+.7f %.7f\n",
                   iter, x_lo, x_hi,
                   r, r - r_expected, 
                   x_hi - x_lo);
         }
       while (status == GSL_CONTINUE && iter < max_iter);
       return status;
     }

Here are the results of the iterations,

     bash$ ./a.out
     using brent method
      iter [    lower,     upper]      root        err  err(est)
         1 [1.0000000, 5.0000000] 1.0000000 -1.2360680 4.0000000
         2 [1.0000000, 3.0000000] 3.0000000 +0.7639320 2.0000000
         3 [2.0000000, 3.0000000] 2.0000000 -0.2360680 1.0000000
         4 [2.2000000, 3.0000000] 2.2000000 -0.0360680 0.8000000
         5 [2.2000000, 2.2366300] 2.2366300 +0.0005621 0.0366300
     Converged:
         6 [2.2360634, 2.2366300] 2.2360634 -0.0000046 0.0005666

If the program is modified to use the bisection solver instead of
Brent's method, by changing `gsl_root_fsolver_brent' to
`gsl_root_fsolver_bisection' the slower convergence of the Bisection
method can be observed,

     bash$ ./a.out
     using bisection method
      iter [    lower,     upper]      root        err  err(est)
         1 [0.0000000, 2.5000000] 1.2500000 -0.9860680 2.5000000
         2 [1.2500000, 2.5000000] 1.8750000 -0.3610680 1.2500000
         3 [1.8750000, 2.5000000] 2.1875000 -0.0485680 0.6250000
         4 [2.1875000, 2.5000000] 2.3437500 +0.1076820 0.3125000
         5 [2.1875000, 2.3437500] 2.2656250 +0.0295570 0.1562500
         6 [2.1875000, 2.2656250] 2.2265625 -0.0095055 0.0781250
         7 [2.2265625, 2.2656250] 2.2460938 +0.0100258 0.0390625
         8 [2.2265625, 2.2460938] 2.2363281 +0.0002601 0.0195312
         9 [2.2265625, 2.2363281] 2.2314453 -0.0046227 0.0097656
        10 [2.2314453, 2.2363281] 2.2338867 -0.0021813 0.0048828
        11 [2.2338867, 2.2363281] 2.2351074 -0.0009606 0.0024414
     Converged:
        12 [2.2351074, 2.2363281] 2.2357178 -0.0003502 0.0012207

   The next program solves the same function using a derivative solver
instead.

     #include <stdio.h>
     #include <gsl/gsl_errno.h>
     #include <gsl/gsl_math.h>
     #include <gsl/gsl_roots.h>
     
     #include "demo_fn.h"
     #include "demo_fn.c"
     
     int
     main (void)
     {
       int status;
       int iter = 0, max_iter = 100;
       const gsl_root_fdfsolver_type *T;
       gsl_root_fdfsolver *s;
       double x0, x = 5.0, r_expected = sqrt (5.0);
       gsl_function_fdf FDF;
       struct quadratic_params params = {1.0, 0.0, -5.0};
     
       FDF.f = &quadratic;
       FDF.df = &quadratic_deriv;
       FDF.fdf = &quadratic_fdf;
       FDF.params = &params;
     
       T = gsl_root_fdfsolver_newton;
       s = gsl_root_fdfsolver_alloc (T);
       gsl_root_fdfsolver_set (s, &FDF, x);
     
       printf ("using %s method\n", 
               gsl_root_fdfsolver_name (s));
     
       printf ("%-5s %10s %10s %10s\n",
               "iter", "root", "err", "err(est)");
       do
         {
           iter++;
           status = gsl_root_fdfsolver_iterate (s);
           x0 = x;
           x = gsl_root_fdfsolver_root (s);
           status = gsl_root_test_delta (x, x0, 0, 1e-3);
     
           if (status == GSL_SUCCESS)
             printf ("Converged:\n");
     
           printf ("%5d %10.7f %+10.7f %10.7f\n",
                   iter, x, x - r_expected, x - x0);
         }
       while (status == GSL_CONTINUE && iter < max_iter);
       return status;
     }

Here are the results for Newton's method,

     bash$ ./a.out
     using newton method
     iter        root        err   err(est)
         1  3.0000000 +0.7639320 -2.0000000
         2  2.3333333 +0.0972654 -0.6666667
         3  2.2380952 +0.0020273 -0.0952381
     Converged:
         4  2.2360689 +0.0000009 -0.0020263

Note that the error can be estimated more accurately by taking the
difference between the current iterate and next iterate rather than the
previous iterate.  The other derivative solvers can be investigated by
changing `gsl_root_fdfsolver_newton' to `gsl_root_fdfsolver_secant' or
`gsl_root_fdfsolver_steffenson'.


File: gsl-ref.info,  Node: Root Finding References and Further Reading,  Prev: Root Finding Examples,  Up: One dimensional Root-Finding

References and Further Reading
==============================

For information on the Brent-Dekker algorithm see the following two
papers,

     R. P. Brent, "An algorithm with guaranteed convergence for finding
     a zero of a function", `Computer Journal', 14 (1971) 422-425

     J. C. P. Bus and T. J. Dekker, "Two Efficient Algorithms with
     Guaranteed Convergence for Finding a Zero of a Function", `ACM
     Transactions of Mathematical Software', Vol. 1 No. 4 (1975) 330-345


File: gsl-ref.info,  Node: One dimensional Minimization,  Next: Multidimensional Root-Finding,  Prev: One dimensional Root-Finding,  Up: Top

One dimensional Minimization
****************************

This chapter describes routines for finding minima of arbitrary
one-dimensional functions.  The library provides low level components
for a variety of iterative minimizers and convergence tests.  These can
be combined by the user to achieve the desired solution, with full
access to the intermediate steps of the algorithms.  Each class of
methods uses the same framework, so that you can switch between
minimizers at runtime without needing to recompile your program.  Each
instance of a minimizer keeps track of its own state, allowing the
minimizers to be used in multi-threaded programs.

   The header file `gsl_min.h' contains prototypes for the minimization
functions and related declarations.  To use the minimization algorithms
to find the maximum of a function simply invert its sign.

* Menu:

* Minimization Overview::
* Minimization Caveats::
* Initializing the Minimizer::
* Providing the function to minimize::
* Minimization Iteration::
* Minimization Stopping Parameters::
* Minimization Algorithms::
* Minimization Examples::
* Minimization References and Further Reading::


File: gsl-ref.info,  Node: Minimization Overview,  Next: Minimization Caveats,  Up: One dimensional Minimization

Overview
========

The minimization algorithms begin with a bounded region known to contain
a minimum.  The region is described by a lower bound a and an upper
bound b, with an estimate of the location of the minimum x.

The value of the function at x must be less than the value of the
function at the ends of the interval,

     f(a) > f(x) < f(b)

This condition guarantees that a minimum is contained somewhere within
the interval.  On each iteration a new point x' is selected using one
of the available algorithms.  If the new point is a better estimate of
the minimum, f(x') < f(x), then the current estimate of the minimum x
is updated.  The new point also allows the size of the bounded interval
to be reduced, by choosing the most compact set of points which
satisfies the constraint f(a) > f(x) < f(b).  The interval is reduced
until it encloses the true minimum to a desired tolerance.  This
provides a best estimate of the location of the minimum and a rigorous
error estimate.

   Several bracketing algorithms are available within a single
framework.  The user provides a high-level driver for the algorithm,
and the library provides the individual functions necessary for each of
the steps.  There are three main phases of the iteration.  The steps
are,

   * initialize minimizer state, S, for algorithm T

   * update S using the iteration T

   * test S for convergence, and repeat iteration if necessary

The state for the minimizers is held in a `gsl_min_fminimizer' struct.
The updating procedure uses only function evaluations (not derivatives).


File: gsl-ref.info,  Node: Minimization Caveats,  Next: Initializing the Minimizer,  Prev: Minimization Overview,  Up: One dimensional Minimization

Caveats
=======

Note that minimization functions can only search for one minimum at a
time.  When there are several minima in the search area, the first
minimum to be found will be returned; however it is difficult to predict
which of the minima this will be. _In most cases, no error will be
reported if you try to find a minimum in an area where there is more
than one._

   With all minimization algorithms it can be difficult to determine the
location of the minimum to full numerical precision.  The behavior of
the function in the region of the minimum x^* can be approximated by a
Taylor expansion,

     y = f(x^*) + (1/2) f''(x^*) (x - x^*)^2

and the second term of this expansion can be lost when added to the
first term at finite precision.  This magnifies the error in locating
x^*, making it proportional to \sqrt \epsilon (where \epsilon is the
relative accuracy of the floating point numbers).  For functions with
higher order minima, such as x^4, the magnification of the error is
correspondingly worse.  The best that can be achieved is to converge to
the limit of numerical accuracy in the function values, rather than the
location of the minimum itself.


File: gsl-ref.info,  Node: Initializing the Minimizer,  Next: Providing the function to minimize,  Prev: Minimization Caveats,  Up: One dimensional Minimization

Initializing the Minimizer
==========================

 - Function: gsl_min_fminimizer * gsl_min_fminimizer_alloc (const
          gsl_min_fminimizer_type * T)
     This function returns a pointer to a newly allocated instance of a
     minimizer of type T.  For example, the following code creates an
     instance of a golden section minimizer,

          const gsl_min_fminimizer_type * T
            = gsl_min_fminimizer_goldensection;
          gsl_min_fminimizer * s
            = gsl_min_fminimizer_alloc (T);

     If there is insufficient memory to create the minimizer then the
     function returns a null pointer and the error handler is invoked
     with an error code of `GSL_ENOMEM'.

 - Function: int gsl_min_fminimizer_set (gsl_min_fminimizer * S,
          gsl_function * F, double X_MINIMUM, double X_LOWER, double
          X_UPPER)
     This function sets, or resets, an existing minimizer S to use the
     function F and the initial search interval [X_LOWER, X_UPPER],
     with a guess for the location of the minimum X_MINIMUM.

     If the interval given does not contain a minimum, then the function
     returns an error code of `GSL_FAILURE'.

 - Function: int gsl_min_fminimizer_set_with_values (gsl_min_fminimizer
          * S, gsl_function * F, double X_MINIMUM, double F_MINIMUM,
          double X_LOWER, double F_LOWER, double X_UPPER, double
          F_UPPER)
     This function is equivalent to `gsl_min_fminimizer_set' but uses
     the values F_MINIMUM, F_LOWER and F_UPPER instead of computing
     `f(x_minimum)', `f(x_lower)' and `f(x_upper)'.

 - Function: void gsl_min_fminimizer_free (gsl_min_fminimizer * S)
     This function frees all the memory associated with the minimizer S.

 - Function: const char * gsl_min_fminimizer_name (const
          gsl_min_fminimizer * S)
     This function returns a pointer to the name of the minimizer.  For
     example,

          printf ("s is a '%s' minimizer\n",
                  gsl_min_fminimizer_name (s));

     would print something like `s is a 'brent' minimizer'.


File: gsl-ref.info,  Node: Providing the function to minimize,  Next: Minimization Iteration,  Prev: Initializing the Minimizer,  Up: One dimensional Minimization

Providing the function to minimize
==================================

You must provide a continuous function of one variable for the
minimizers to operate on.  In order to allow for general parameters the
functions are defined by a `gsl_function' data type (*note Providing
the function to solve::).


File: gsl-ref.info,  Node: Minimization Iteration,  Next: Minimization Stopping Parameters,  Prev: Providing the function to minimize,  Up: One dimensional Minimization

Iteration
=========

The following functions drive the iteration of each algorithm.  Each
function performs one iteration to update the state of any minimizer of
the corresponding type.  The same functions work for all minimizers so
that different methods can be substituted at runtime without
modifications to the code.

 - Function: int gsl_min_fminimizer_iterate (gsl_min_fminimizer * S)
     This function performs a single iteration of the minimizer S.  If
     the iteration encounters an unexpected problem then an error code
     will be returned,

    `GSL_EBADFUNC'
          the iteration encountered a singular point where the function
          evaluated to `Inf' or `NaN'.

    `GSL_FAILURE'
          the algorithm could not improve the current best
          approximation or bounding interval.

   The minimizer maintains a current best estimate of the position of
the minimum at all times, and the current interval bounding the minimum.
This information can be accessed with the following auxiliary functions,

 - Function: double gsl_min_fminimizer_x_minimum (const
          gsl_min_fminimizer * S)
     This function returns the current estimate of the position of the
     minimum for the minimizer S.

 - Function: double gsl_min_fminimizer_x_upper (const
          gsl_min_fminimizer * S)
 - Function: double gsl_min_fminimizer_x_lower (const
          gsl_min_fminimizer * S)
     These functions return the current upper and lower bound of the
     interval for the minimizer S.

 - Function: double gsl_min_fminimizer_f_minimum (const
          gsl_min_fminimizer *S)
 - Function: double gsl_min_fminimizer_f_upper (const
          gsl_min_fminimizer *S)
 - Function: double gsl_min_fminimizer_f_lower (const
          gsl_min_fminimizer *S)
     These functions return the value of the function at the current
     estimate of the minimum and at the upper and lower bounds of
     interval for the minimizer S.


File: gsl-ref.info,  Node: Minimization Stopping Parameters,  Next: Minimization Algorithms,  Prev: Minimization Iteration,  Up: One dimensional Minimization

Stopping Parameters
===================

A minimization procedure should stop when one of the following
conditions is true:

   * A minimum has been found to within the user-specified precision.

   * A user-specified maximum number of iterations has been reached.

   * An error has occurred.

The handling of these conditions is under user control.  The function
below allows the user to test the precision of the current result.

 - Function: int gsl_min_test_interval (double X_LOWER, double X_UPPER,
          double EPSABS, double EPSREL)
     This function tests for the convergence of the interval [X_LOWER,
     X_UPPER] with absolute error EPSABS and relative error EPSREL.
     The test returns `GSL_SUCCESS' if the following condition is
     achieved,

          |a - b| < epsabs + epsrel min(|a|,|b|)

     when the interval x = [a,b] does not include the origin.  If the
     interval includes the origin then \min(|a|,|b|) is replaced by
     zero (which is the minimum value of |x| over the interval).  This
     ensures that the relative error is accurately estimated for minima
     close to the origin.

     This condition on the interval also implies that any estimate of
     the minimum x_m in the interval satisfies the same condition with
     respect to the true minimum x_m^*,

          |x_m - x_m^*| < epsabs + epsrel x_m^*

     assuming that the true minimum x_m^* is contained within the
     interval.


File: gsl-ref.info,  Node: Minimization Algorithms,  Next: Minimization Examples,  Prev: Minimization Stopping Parameters,  Up: One dimensional Minimization

Minimization Algorithms
=======================

The minimization algorithms described in this section require an initial
interval which is guaranteed to contain a minimum -- if a and b are the
endpoints of the interval and x is an estimate of the minimum then f(a)
> f(x) < f(b).  This ensures that the function has at least one minimum
somewhere in the interval.  If a valid initial interval is used then
these algorithm cannot fail, provided the function is well-behaved.

 - Minimizer: gsl_min_fminimizer_goldensection
     The "golden section algorithm" is the simplest method of bracketing
     the minimum of a function.  It is the slowest algorithm provided
     by the library, with linear convergence.

     On each iteration, the algorithm first compares the subintervals
     from the endpoints to the current minimum.  The larger subinterval
     is divided in a golden section (using the famous ratio (3-\sqrt
     5)/2 = 0.3189660...) and the value of the function at this new
     point is calculated.  The new value is used with the constraint
     f(a') > f(x') < f(b') to a select new interval containing the
     minimum, by discarding the least useful point.  This procedure can
     be continued indefinitely until the interval is sufficiently
     small.  Choosing the golden section as the bisection ratio can be
     shown to provide the fastest convergence for this type of
     algorithm.


 - Minimizer: gsl_min_fminimizer_brent
     The "Brent minimization algorithm" combines a parabolic
     interpolation with the golden section algorithm.  This produces a
     fast algorithm which is still robust.

     The outline of the algorithm can be summarized as follows: on each
     iteration Brent's method approximates the function using an
     interpolating parabola through three existing points.  The minimum
     of the parabola is taken as a guess for the minimum.  If it lies
     within the bounds of the current interval then the interpolating
     point is accepted, and used to generate a smaller interval.  If
     the interpolating point is not accepted then the algorithm falls
     back to an ordinary golden section step.  The full details of
     Brent's method include some additional checks to improve
     convergence.


File: gsl-ref.info,  Node: Minimization Examples,  Next: Minimization References and Further Reading,  Prev: Minimization Algorithms,  Up: One dimensional Minimization

Examples
========

The following program uses the Brent algorithm to find the minimum of
the function f(x) = \cos(x) + 1, which occurs at x = \pi.  The starting
interval is (0,6), with an initial guess for the minimum of 2.

     #include <stdio.h>
     #include <gsl/gsl_errno.h>
     #include <gsl/gsl_math.h>
     #include <gsl/gsl_min.h>
     
     double fn1 (double x, void * params)
     {
       return cos(x) + 1.0;
     }
     
     int
     main (void)
     {
       int status;
       int iter = 0, max_iter = 100;
       const gsl_min_fminimizer_type *T;
       gsl_min_fminimizer *s;
       double m = 2.0, m_expected = M_PI;
       double a = 0.0, b = 6.0;
       gsl_function F;
     
       F.function = &fn1;
       F.params = 0;
     
       T = gsl_min_fminimizer_brent;
       s = gsl_min_fminimizer_alloc (T);
       gsl_min_fminimizer_set (s, &F, m, a, b);
     
       printf ("using %s method\n",
               gsl_min_fminimizer_name (s));
     
       printf ("%5s [%9s, %9s] %9s %10s %9s\n",
               "iter", "lower", "upper", "min",
               "err", "err(est)");
     
       printf ("%5d [%.7f, %.7f] %.7f %+.7f %.7f\n",
               iter, a, b,
               m, m - m_expected, b - a);
     
       do
         {
           iter++;
           status = gsl_min_fminimizer_iterate (s);
     
           m = gsl_min_fminimizer_x_minimum (s);
           a = gsl_min_fminimizer_x_lower (s);
           b = gsl_min_fminimizer_x_upper (s);
     
           status 
             = gsl_min_test_interval (a, b, 0.001, 0.0);
     
           if (status == GSL_SUCCESS)
             printf ("Converged:\n");
     
           printf ("%5d [%.7f, %.7f] "
                   "%.7f %.7f %+.7f %.7f\n",
                   iter, a, b,
                   m, m_expected, m - m_expected, b - a);
         }
       while (status == GSL_CONTINUE && iter < max_iter);
     
       return status;
     }

Here are the results of the minimization procedure.

     bash$ ./a.out

         0 [0.0000000, 6.0000000] 2.0000000 -1.1415927 6.0000000
         1 [2.0000000, 6.0000000] 3.2758640 +0.1342713 4.0000000
         2 [2.0000000, 3.2831929] 3.2758640 +0.1342713 1.2831929
         3 [2.8689068, 3.2831929] 3.2758640 +0.1342713 0.4142862
         4 [2.8689068, 3.2831929] 3.2758640 +0.1342713 0.4142862
         5 [2.8689068, 3.2758640] 3.1460585 +0.0044658 0.4069572
         6 [3.1346075, 3.2758640] 3.1460585 +0.0044658 0.1412565
         7 [3.1346075, 3.1874620] 3.1460585 +0.0044658 0.0528545
         8 [3.1346075, 3.1460585] 3.1460585 +0.0044658 0.0114510
         9 [3.1346075, 3.1460585] 3.1424060 +0.0008133 0.0114510
        10 [3.1346075, 3.1424060] 3.1415885 -0.0000041 0.0077985
     Converged:                            
        11 [3.1415885, 3.1424060] 3.1415927 -0.0000000 0.0008175


File: gsl-ref.info,  Node: Minimization References and Further Reading,  Prev: Minimization Examples,  Up: One dimensional Minimization

References and Further Reading
==============================

Further information on Brent's algorithm is available in the following
book,

     Richard Brent, `Algorithms for minimization without derivatives',
     Prentice-Hall (1973), republished by Dover in paperback (2002),
     ISBN 0-486-41998-3.


File: gsl-ref.info,  Node: Multidimensional Root-Finding,  Next: Multidimensional Minimization,  Prev: One dimensional Minimization,  Up: Top

Multidimensional Root-Finding
*****************************

This chapter describes functions for multidimensional root-finding
(solving nonlinear systems with n equations in n unknowns).  The
library provides low level components for a variety of iterative
solvers and convergence tests.  These can be combined by the user to
achieve the desired solution, with full access to the intermediate
steps of the iteration.  Each class of methods uses the same framework,
so that you can switch between solvers at runtime without needing to
recompile your program.  Each instance of a solver keeps track of its
own state, allowing the solvers to be used in multi-threaded programs.
The solvers are based on the original Fortran library MINPACK.

   The header file `gsl_multiroots.h' contains prototypes for the
multidimensional root finding functions and related declarations.

* Menu:

* Overview of Multidimensional Root Finding::
* Initializing the Multidimensional Solver::
* Providing the multidimensional system of equations to solve::
* Iteration of the multidimensional solver::
* Search Stopping Parameters for the multidimensional solver::
* Algorithms using Derivatives::
* Algorithms without Derivatives::
* Example programs for Multidimensional Root finding::
* References and Further Reading for Multidimensional Root Finding::


File: gsl-ref.info,  Node: Overview of Multidimensional Root Finding,  Next: Initializing the Multidimensional Solver,  Up: Multidimensional Root-Finding

Overview
========

The problem of multidimensional root finding requires the simultaneous
solution of n equations, f_i, in n variables, x_i,

     f_i (x_1, ..., x_n) = 0    for i = 1 ... n.

In general there are no bracketing methods available for n dimensional
systems, and no way of knowing whether any solutions exist.  All
algorithms proceed from an initial guess using a variant of the Newton
iteration,

     x -> x' = x - J^{-1} f(x)

where x, f are vector quantities and J is the Jacobian matrix J_{ij} =
d f_i / d x_j.  Additional strategies can be used to enlarge the region
of convergence.  These include requiring a decrease in the norm |f| on
each step proposed by Newton's method, or taking steepest-descent steps
in the direction of the negative gradient of |f|.

   Several root-finding algorithms are available within a single
framework.  The user provides a high-level driver for the algorithms,
and the library provides the individual functions necessary for each of
the steps.  There are three main phases of the iteration.  The steps
are,

   * initialize solver state, S, for algorithm T

   * update S using the iteration T

   * test S for convergence, and repeat iteration if necessary

The evaluation of the Jacobian matrix can be problematic, either because
programming the derivatives is intractable or because computation of the
n^2 terms of the matrix becomes too expensive.  For these reasons the
algorithms provided by the library are divided into two classes
according to whether the derivatives are available or not.

   The state for solvers with an analytic Jacobian matrix is held in a
`gsl_multiroot_fdfsolver' struct.  The updating procedure requires both
the function and its derivatives to be supplied by the user.

   The state for solvers which do not use an analytic Jacobian matrix is
held in a `gsl_multiroot_fsolver' struct.  The updating procedure uses
only function evaluations (not derivatives).  The algorithms estimate
the matrix J or J^{-1} by approximate methods.


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

Initializing the Solver
=======================

The following functions initialize a multidimensional solver, either
with or without derivatives.  The solver itself depends only on the
dimension of the problem and the algorithm and can be reused for
different problems.

 - Function: gsl_multiroot_fsolver * gsl_multiroot_fsolver_alloc (const
          gsl_multiroot_fsolver_type * T, size_t N)
     This function returns a pointer to a newly allocated instance of a
     solver of type T for a system of N dimensions.  For example, the
     following code creates an instance of a hybrid solver, to solve a
     3-dimensional system of equations.

          const gsl_multiroot_fsolver_type * T
              = gsl_multiroot_fsolver_hybrid;
          gsl_multiroot_fsolver * s
              = gsl_multiroot_fsolver_alloc (T, 3);

     If there is insufficient memory to create the solver then the
     function returns a null pointer and the error handler is invoked
     with an error code of `GSL_ENOMEM'.

 - Function: gsl_multiroot_fdfsolver * gsl_multiroot_fdfsolver_alloc
          (const gsl_multiroot_fdfsolver_type * T, size_t N)
     This function returns a pointer to a newly allocated instance of a
     derivative solver of type T for a system of N dimensions.  For
     example, the following code creates an instance of a
     Newton-Raphson solver, for a 2-dimensional system of equations.

          const gsl_multiroot_fdfsolver_type * T
              = gsl_multiroot_fdfsolver_newton;
          gsl_multiroot_fdfsolver * s =
              gsl_multiroot_fdfsolver_alloc (T, 2);

     If there is insufficient memory to create the solver then the
     function returns a null pointer and the error handler is invoked
     with an error code of `GSL_ENOMEM'.

 - Function: int gsl_multiroot_fsolver_set (gsl_multiroot_fsolver * S,
          gsl_multiroot_function * F, gsl_vector * X)
     This function sets, or resets, an existing solver S to use the
     function F and the initial guess X.

 - Function: int gsl_multiroot_fdfsolver_set (gsl_multiroot_fdfsolver *
          S, gsl_multiroot_function_fdf * FDF, gsl_vector * X)
     This function sets, or resets, an existing solver S to use the
     function and derivative FDF and the initial guess X.

 - Function: void gsl_multiroot_fsolver_free (gsl_multiroot_fsolver * S)
 - Function: void gsl_multiroot_fdfsolver_free (gsl_multiroot_fdfsolver
          * S)
     These functions free all the memory associated with the solver S.

 - Function: const char * gsl_multiroot_fsolver_name (const
          gsl_multiroot_fsolver * S)
 - Function: const char * gsl_multiroot_fdfsolver_name (const
          gsl_multiroot_fdfsolver * S)
     These functions return a pointer to the name of the solver.  For
     example,

          printf ("s is a '%s' solver\n",
                  gsl_multiroot_fdfsolver_name (s));

     would print something like `s is a 'newton' solver'.


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

Providing the function to solve
===============================

You must provide n functions of n variables for the root finders to
operate on.  In order to allow for general parameters the functions are
defined by the following data types:

 - Data Type: gsl_multiroot_function
     This data type defines a general system of functions with
     parameters.

    `int (* f) (const gsl_vector * X, void * PARAMS, gsl_vector * F)'
          this function should store the vector result f(x,params) in F
          for argument X and parameters PARAMS, returning an
          appropriate error code if the function cannot be computed.

    `size_t n'
          the dimension of the system, i.e. the number of components of
          the vectors X and F.

    `void * params'
          a pointer to the parameters of the function.

Here is an example using Powell's test function,

     f_1(x) = A x_0 x_1 - 1,
     f_2(x) = exp(-x_0) + exp(-x_1) - (1 + 1/A)

with A = 10^4.  The following code defines a `gsl_multiroot_function'
system `F' which you could pass to a solver:

     struct powell_params { double A; };
     
     int
     powell (gsl_vector * x, void * p, gsl_vector * f) {
        struct powell_params * params
          = *(struct powell_params *)p;
        const double A = (params->A);
        const double x0 = gsl_vector_get(x,0);
        const double x1 = gsl_vector_get(x,1);
     
        gsl_vector_set (f, 0, A * x0 * x1 - 1);
        gsl_vector_set (f, 1, (exp(-x0) + exp(-x1)
                               - (1.0 + 1.0/A)));
        return GSL_SUCCESS
     }
     
     gsl_multiroot_function F;
     struct powell_params params = { 10000.0 };
     
     F.f = &powell;
     F.n = 2;
     F.params = &params;

 - Data Type: gsl_multiroot_function_fdf
     This data type defines a general system of functions with
     parameters and the corresponding Jacobian matrix of derivatives,

    `int (* f) (const gsl_vector * X, void * PARAMS, gsl_vector * F)'
          this function should store the vector result f(x,params) in F
          for argument X and parameters PARAMS, returning an
          appropriate error code if the function cannot be computed.

    `int (* df) (const gsl_vector * X, void * PARAMS, gsl_matrix * J)'
          this function should store the N-by-N matrix result J_ij = d
          f_i(x,params) / d x_j in J for argument X and parameters
          PARAMS, returning an appropriate error code if the function
          cannot be computed.

    `int (* fdf) (const gsl_vector * X, void * PARAMS, gsl_vector * F, gsl_matrix * J)'
          This function should set the values of the F and J as above,
          for arguments X and parameters PARAMS.  This function
          provides an optimization of the separate functions for f(x)
          and J(x)--it is always faster to compute the function and its
          derivative at the same time.

    `size_t n'
          the dimension of the system, i.e. the number of components of
          the vectors X and F.

    `void * params'
          a pointer to the parameters of the function.

The example of Powell's test function defined above can be extended to
include analytic derivatives using the following code,

     int
     powell_df (gsl_vector * x, void * p, gsl_matrix * J)
     {
        struct powell_params * params
          = *(struct powell_params *)p;
        const double A = (params->A);
        const double x0 = gsl_vector_get(x,0);
        const double x1 = gsl_vector_get(x,1);
        gsl_matrix_set (J, 0, 0, A * x1);
        gsl_matrix_set (J, 0, 1, A * x0);
        gsl_matrix_set (J, 1, 0, -exp(-x0));
        gsl_matrix_set (J, 1, 1, -exp(-x1));
        return GSL_SUCCESS
     }
     
     int
     powell_fdf (gsl_vector * x, void * p,
                 gsl_matrix * f, gsl_matrix * J) {
        struct powell_params * params
          = *(struct powell_params *)p;
        const double A = (params->A);
        const double x0 = gsl_vector_get(x,0);
        const double x1 = gsl_vector_get(x,1);
     
        const double u0 = exp(-x0);
        const double u1 = exp(-x1);
     
        gsl_vector_set (f, 0, A * x0 * x1 - 1);
        gsl_vector_set (f, 1, u0 + u1 - (1 + 1/A));
     
        gsl_matrix_set (J, 0, 0, A * x1);
        gsl_matrix_set (J, 0, 1, A * x0);
        gsl_matrix_set (J, 1, 0, -u0);
        gsl_matrix_set (J, 1, 1, -u1);
        return GSL_SUCCESS
     }
     
     gsl_multiroot_function_fdf FDF;
     
     FDF.f = &powell_f;
     FDF.df = &powell_df;
     FDF.fdf = &powell_fdf;
     FDF.n = 2;
     FDF.params = 0;

Note that the function `powell_fdf' is able to reuse existing terms
from the function when calculating the Jacobian, thus saving time.


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

Iteration
=========

The following functions drive the iteration of each algorithm.  Each
function performs one iteration to update the state of any solver of the
corresponding type.  The same functions work for all solvers so that
different methods can be substituted at runtime without modifications to
the code.

 - Function: int gsl_multiroot_fsolver_iterate (gsl_multiroot_fsolver *
          S)
 - Function: int gsl_multiroot_fdfsolver_iterate
          (gsl_multiroot_fdfsolver * S)
     These functions perform a single iteration of the solver S.  If the
     iteration encounters an unexpected problem then an error code will
     be returned,

    `GSL_EBADFUNC'
          the iteration encountered a singular point where the function
          or its derivative evaluated to `Inf' or `NaN'.

    `GSL_ENOPROG'
          the iteration is not making any progress, preventing the
          algorithm from continuing.

   The solver maintains a current best estimate of the root at all
times.  This information can be accessed with the following auxiliary
functions,

 - Function: gsl_vector * gsl_multiroot_fsolver_root (const
          gsl_multiroot_fsolver * S)
 - Function: gsl_vector * gsl_multiroot_fdfsolver_root (const
          gsl_multiroot_fdfsolver * S)
     These functions return the current estimate of the root for the
     solver S.

 - Function: gsl_vector * gsl_multiroot_fsolver_f (const
          gsl_multiroot_fsolver * S)
 - Function: gsl_vector * gsl_multiroot_fdfsolver_f (const
          gsl_multiroot_fdfsolver * S)
     These functions return the function value f(x) at the current
     estimate of the root for the solver S.

 - Function: gsl_vector * gsl_multiroot_fsolver_dx (const
          gsl_multiroot_fsolver * S)
 - Function: gsl_vector * gsl_multiroot_fdfsolver_dx (const
          gsl_multiroot_fdfsolver * S)
     These functions return the last step dx taken by the solver S.


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

Search Stopping Parameters
==========================

A root finding procedure should stop when one of the following
conditions is true:

   * A multidimensional root has been found to within the
     user-specified precision.

   * A user-specified maximum number of iterations has been reached.

   * An error has occurred.

The handling of these conditions is under user control.  The functions
below allow the user to test the precision of the current result in
several standard ways.

 - Function: int gsl_multiroot_test_delta (const gsl_vector * DX, const
          gsl_vector * X, double EPSABS, double EPSREL)
     This function tests for the convergence of the sequence by
     comparing the last step DX with the absolute error EPSABS and
     relative error EPSREL to the current position X.  The test returns
     `GSL_SUCCESS' if the following condition is achieved,

          |dx_i| < epsabs + epsrel |x_i|

     for each component of X and returns `GSL_CONTINUE' otherwise.

 - Function: int gsl_multiroot_test_residual (const gsl_vector * F,
          double EPSABS)
     This function tests the residual value F against the absolute
     error bound EPSABS.  The test returns `GSL_SUCCESS' if the
     following condition is achieved,

          \sum_i |f_i| < epsabs

     and returns `GSL_CONTINUE' otherwise.  This criterion is suitable
     for situations where the precise location of the root, x, is
     unimportant provided a value can be found where the residual is
     small enough.


File: gsl-ref.info,  Node: Algorithms using Derivatives,  Next: Algorithms without Derivatives,  Prev: Search Stopping Parameters for the multidimensional solver,  Up: Multidimensional Root-Finding

Algorithms using Derivatives
============================

The root finding algorithms described in this section make use of both
the function and its derivative.  They require an initial guess for the
location of the root, but there is no absolute guarantee of
convergence--the function must be suitable for this technique and the
initial guess must be sufficiently close to the root for it to work.
When the conditions are satisfied then convergence is quadratic.

 - Derivative Solver: gsl_multiroot_fdfsolver_hybridsj
     This is a modified version of Powell's Hybrid method as
     implemented in the HYBRJ algorithm in MINPACK.  Minpack was
     written by Jorge J. More', Burton S. Garbow and Kenneth E.
     Hillstrom.  The Hybrid algorithm retains the fast convergence of
     Newton's method but will also reduce the residual when Newton's
     method is unreliable.

     The algorithm uses a generalized trust region to keep each step
     under control.  In order to be accepted a proposed new position x'
     must satisfy the condition |D (x' - x)| < \delta, where D is a
     diagonal scaling matrix and \delta is the size of the trust
     region.  The components of D are computed internally, using the
     column norms of the Jacobian to estimate the sensitivity of the
     residual to each component of x.  This improves the behavior of the
     algorithm for badly scaled functions.

     On each iteration the algorithm first determines the standard
     Newton step by solving the system J dx = - f.  If this step falls
     inside the trust region it is used as a trial step in the next
     stage.  If not, the algorithm uses the linear combination of the
     Newton and gradient directions which is predicted to minimize the
     norm of the function while staying inside the trust region,

          dx = - \alpha J^{-1} f(x) - \beta \nabla |f(x)|^2.

     This combination of Newton and gradient directions is referred to
     as a "dogleg step".

     The proposed step is now tested by evaluating the function at the
     resulting point, x'.  If the step reduces the norm of the function
     sufficiently then it is accepted and size of the trust region is
     increased.  If the proposed step fails to improve the solution
     then the size of the trust region is decreased and another trial
     step is computed.

     The speed of the algorithm is increased by computing the changes
     to the Jacobian approximately, using a rank-1 update.  If two
     successive attempts fail to reduce the residual then the full
     Jacobian is recomputed.  The algorithm also monitors the progress
     of the solution and returns an error if several steps fail to make
     any improvement,

    `GSL_ENOPROG'
          the iteration is not making any progress, preventing the
          algorithm from continuing.

    `GSL_ENOPROGJ'
          re-evaluations of the Jacobian indicate that the iteration is
          not making any progress, preventing the algorithm from
          continuing.


 - Derivative Solver: gsl_multiroot_fdfsolver_hybridj
     This algorithm is an unscaled version of `hybridsj'.  The steps are
     controlled by a spherical trust region |x' - x| < \delta, instead
     of a generalized region.  This can be useful if the generalized
     region estimated by `hybridsj' is inappropriate.

 - Derivative Solver: gsl_multiroot_fdfsolver_newton
     Newton's Method is the standard root-polishing algorithm.  The
     algorithm begins with an initial guess for the location of the
     solution.  On each iteration a linear approximation to the
     function F is used to estimate the step which will zero all the
     components of the residual.  The iteration is defined by the
     following sequence,

          x -> x' = x - J^{-1} f(x)

     where the Jacobian matrix J is computed from the derivative
     functions provided by F.  The step dx is obtained by solving the
     linear system,

          J dx = - f(x)

     using LU decomposition.

 - Derivative Solver: gsl_multiroot_fdfsolver_gnewton
     This is a modified version of Newton's method which attempts to
     improve global convergence by requiring every step to reduce the
     Euclidean norm of the residual, |f(x)|.  If the Newton step leads
     to an increase in the norm then a reduced step of relative size,

          t = (\sqrt(1 + 6 r) - 1) / (3 r)

     is proposed, with r being the ratio of norms |f(x')|^2/|f(x)|^2.
     This procedure is repeated until a suitable step size is found.


File: gsl-ref.info,  Node: Algorithms without Derivatives,  Next: Example programs for Multidimensional Root finding,  Prev: Algorithms using Derivatives,  Up: Multidimensional Root-Finding

Algorithms without Derivatives
==============================

The algorithms described in this section do not require any derivative
information to be supplied by the user.  Any derivatives needed are
approximated from by finite difference.

 - Solver: gsl_multiroot_fsolver_hybrids
     This is a version of the Hybrid algorithm which replaces calls to
     the Jacobian function by its finite difference approximation.  The
     finite difference approximation is computed using
     `gsl_multiroots_fdjac' with a relative step size of
     `GSL_SQRT_DBL_EPSILON'.

 - Solver: gsl_multiroot_fsolver_hybrid
     This is a finite difference version of the Hybrid algorithm without
     internal scaling.

 - Solver: gsl_multiroot_fsolver_dnewton
     The "discrete Newton algorithm" is the simplest method of solving a
     multidimensional system.  It uses the Newton iteration

          x -> x - J^{-1} f(x)

     where the Jacobian matrix J is approximated by taking finite
     differences of the function F.  The approximation scheme used by
     this implementation is,

          J_{ij} = (f_i(x + \delta_j) - f_i(x)) /  \delta_j

     where \delta_j is a step of size \sqrt\epsilon |x_j| with \epsilon
     being the machine precision (\epsilon \approx 2.22 \times 10^-16).
     The order of convergence of Newton's algorithm is quadratic, but
     the finite differences require n^2 function evaluations on each
     iteration.  The algorithm may become unstable if the finite
     differences are not a good approximation to the true derivatives.

 - Solver: gsl_multiroot_fsolver_broyden
     The "Broyden algorithm" is a version of the discrete Newton
     algorithm which attempts to avoids the expensive update of the
     Jacobian matrix on each iteration.  The changes to the Jacobian
     are also approximated, using a rank-1 update,

          J^{-1} \to J^{-1} - (J^{-1} df - dx) dx^T J^{-1} / dx^T J^{-1} df

     where the vectors dx and df are the changes in x and f.  On the
     first iteration the inverse Jacobian is estimated using finite
     differences, as in the discrete Newton algorithm.

     This approximation gives a fast update but is unreliable if the
     changes are not small, and the estimate of the inverse Jacobian
     becomes worse as time passes.  The algorithm has a tendency to
     become unstable unless it starts close to the root.  The Jacobian
     is refreshed if this instability is detected (consult the source
     for details).

     This algorithm is not recommended and is included only for
     demonstration purposes.


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

Examples
========

The multidimensional solvers are used in a similar way to the
one-dimensional root finding algorithms.  This first example
demonstrates the `hybrids' scaled-hybrid algorithm, which does not
require derivatives. The program solves the Rosenbrock system of
equations,

     f_1 (x, y) = a (1 - x)
     f_2 (x, y) = b (y - x^2)

with a = 1, b = 10. The solution of this system lies at (x,y) = (1,1)
in a narrow valley.

   The first stage of the program is to define the system of equations,

     #include <stdlib.h>
     #include <stdio.h>
     #include <gsl/gsl_vector.h>
     #include <gsl/gsl_multiroots.h>
     
     struct rparams
       {
         double a;
         double b;
       };
     
     int
     rosenbrock_f (const gsl_vector * x, void *params,
                   gsl_vector * f)
     {
       double a = ((struct rparams *) params)->a;
       double b = ((struct rparams *) params)->b;
     
       const double x0 = gsl_vector_get (x, 0);
       const double x1 = gsl_vector_get (x, 1);
     
       const double y0 = a * (1 - x0);
       const double y1 = b * (x1 - x0 * x0);
     
       gsl_vector_set (f, 0, y0);
       gsl_vector_set (f, 1, y1);
     
       return GSL_SUCCESS;
     }

The main program begins by creating the function object `f', with the
arguments `(x,y)' and parameters `(a,b)'. The solver `s' is initialized
to use this function, with the `hybrids' method.

     int
     main (void)
     {
       const gsl_multiroot_fsolver_type *T;
       gsl_multiroot_fsolver *s;
     
       int status;
       size_t i, iter = 0;
     
       const size_t n = 2;
       struct rparams p = {1.0, 10.0};
       gsl_multiroot_function f = {&rosenbrock_f, n, &p};
     
       double x_init[2] = {-10.0, -5.0};
       gsl_vector *x = gsl_vector_alloc (n);
     
       gsl_vector_set (x, 0, x_init[0]);
       gsl_vector_set (x, 1, x_init[1]);
     
       T = gsl_multiroot_fsolver_hybrids;
       s = gsl_multiroot_fsolver_alloc (T, 2);
       gsl_multiroot_fsolver_set (s, &f, x);
     
       print_state (iter, s);
     
       do
         {
           iter++;
           status = gsl_multiroot_fsolver_iterate (s);
     
           print_state (iter, s);
     
           if (status)   /* check if solver is stuck */
             break;
     
           status =
             gsl_multiroot_test_residual (s->f, 1e-7);
         }
       while (status == GSL_CONTINUE && iter < 1000);
     
       printf ("status = %s\n", gsl_strerror (status));
     
       gsl_multiroot_fsolver_free (s);
       gsl_vector_free (x);
       return 0;
     }

Note that it is important to check the return status of each solver
step, in case the algorithm becomes stuck.  If an error condition is
detected, indicating that the algorithm cannot proceed, then the error
can be reported to the user, a new starting point chosen or a different
algorithm used.

   The intermediate state of the solution is displayed by the following
function.  The solver state contains the vector `s->x' which is the
current position, and the vector `s->f' with corresponding function
values.
     int
     print_state (size_t iter, gsl_multiroot_fsolver * s)
     {
       printf ("iter = %3u x = % .3f % .3f "
               "f(x) = % .3e % .3e\n",
               iter,
               gsl_vector_get (s->x, 0),
               gsl_vector_get (s->x, 1),
               gsl_vector_get (s->f, 0),
               gsl_vector_get (s->f, 1));
     }

Here are the results of running the program. The algorithm is started at
(-10,-5) far from the solution.  Since the solution is hidden in a
narrow valley the earliest steps follow the gradient of the function
downhill, in an attempt to reduce the large value of the residual. Once
the root has been approximately located, on iteration 8, the Newton
behavior takes over and convergence is very rapid.

     iter =  0 x = -10.000  -5.000  f(x) = 1.100e+01 -1.050e+03
     iter =  1 x = -10.000  -5.000  f(x) = 1.100e+01 -1.050e+03
     iter =  2 x =  -3.976  24.827  f(x) = 4.976e+00  9.020e+01
     iter =  3 x =  -3.976  24.827  f(x) = 4.976e+00  9.020e+01
     iter =  4 x =  -3.976  24.827  f(x) = 4.976e+00  9.020e+01
     iter =  5 x =  -1.274  -5.680  f(x) = 2.274e+00 -7.302e+01
     iter =  6 x =  -1.274  -5.680  f(x) = 2.274e+00 -7.302e+01
     iter =  7 x =   0.249   0.298  f(x) = 7.511e-01  2.359e+00
     iter =  8 x =   0.249   0.298  f(x) = 7.511e-01  2.359e+00
     iter =  9 x =   1.000   0.878  f(x) = 1.268e-10 -1.218e+00
     iter = 10 x =   1.000   0.989  f(x) = 1.124e-11 -1.080e-01
     iter = 11 x =   1.000   1.000  f(x) = 0.000e+00  0.000e+00
     status = success

Note that the algorithm does not update the location on every
iteration. Some iterations are used to adjust the trust-region
parameter, after trying a step which was found to be divergent, or to
recompute the Jacobian, when poor convergence behavior is detected.

   The next example program adds derivative information, in order to
accelerate the solution. There are two derivative functions
`rosenbrock_df' and `rosenbrock_fdf'. The latter computes both the
function and its derivative simultaneously. This allows the
optimization of any common terms.  For simplicity we substitute calls to
the separate `f' and `df' functions at this point in the code below.

     int
     rosenbrock_df (const gsl_vector * x, void *params,
                    gsl_matrix * J)
     {
       const double a = ((struct rparams *) params)->a;
       const double b = ((struct rparams *) params)->b;
     
       const double x0 = gsl_vector_get (x, 0);
     
       const double df00 = -a;
       const double df01 = 0;
       const double df10 = -2 * b  * x0;
       const double df11 = b;
     
       gsl_matrix_set (J, 0, 0, df00);
       gsl_matrix_set (J, 0, 1, df01);
       gsl_matrix_set (J, 1, 0, df10);
       gsl_matrix_set (J, 1, 1, df11);
     
       return GSL_SUCCESS;
     }
     
     int
     rosenbrock_fdf (const gsl_vector * x, void *params,
                     gsl_vector * f, gsl_matrix * J)
     {
       rosenbrock_f (x, params, f);
       rosenbrock_df (x, params, J);
     
       return GSL_SUCCESS;
     }

The main program now makes calls to the corresponding `fdfsolver'
versions of the functions,

     int
     main (void)
     {
       const gsl_multiroot_fdfsolver_type *T;
       gsl_multiroot_fdfsolver *s;
     
       int status;
       size_t i, iter = 0;
     
       const size_t n = 2;
       struct rparams p = {1.0, 10.0};
       gsl_multiroot_function_fdf f = {&rosenbrock_f,
                                       &rosenbrock_df,
                                       &rosenbrock_fdf,
                                       n, &p};
     
       double x_init[2] = {-10.0, -5.0};
       gsl_vector *x = gsl_vector_alloc (n);
     
       gsl_vector_set (x, 0, x_init[0]);
       gsl_vector_set (x, 1, x_init[1]);
     
       T = gsl_multiroot_fdfsolver_gnewton;
       s = gsl_multiroot_fdfsolver_alloc (T, n);
       gsl_multiroot_fdfsolver_set (s, &f, x);
     
       print_state (iter, s);
     
       do
         {
           iter++;
     
           status = gsl_multiroot_fdfsolver_iterate (s);
     
           print_state (iter, s);
     
           if (status)
             break;
     
           status = gsl_multiroot_test_residual (s->f, 1e-7);
         }
       while (status == GSL_CONTINUE && iter < 1000);
     
       printf ("status = %s\n", gsl_strerror (status));
     
       gsl_multiroot_fdfsolver_free (s);
       gsl_vector_free (x);
       return 0;
     }

The addition of derivative information to the `hybrids' solver does not
make any significant difference to its behavior, since it able to
approximate the Jacobian numerically with sufficient accuracy.  To
illustrate the behavior of a different derivative solver we switch to
`gnewton'. This is a traditional Newton solver with the constraint that
it scales back its step if the full step would lead "uphill". Here is
the output for the `gnewton' algorithm,

     iter = 0 x = -10.000  -5.000 f(x) =  1.100e+01 -1.050e+03
     iter = 1 x =  -4.231 -65.317 f(x) =  5.231e+00 -8.321e+02
     iter = 2 x =   1.000 -26.358 f(x) = -8.882e-16 -2.736e+02
     iter = 3 x =   1.000   1.000 f(x) = -2.220e-16 -4.441e-15
     status = success

The convergence is much more rapid, but takes a wide excursion out to
the point (-4.23,-65.3). This could cause the algorithm to go astray in
a realistic application.  The hybrid algorithm follows the downhill
path to the solution more reliably.


File: gsl-ref.info,  Node: References and Further Reading for Multidimensional Root Finding,  Prev: Example programs for Multidimensional Root finding,  Up: Multidimensional Root-Finding

References and Further Reading
==============================

The original version of the Hybrid method is described in the following
articles by Powell,

     M.J.D. Powell, "A Hybrid Method for Nonlinear Equations" (Chap 6, p
     87-114) and "A Fortran Subroutine for Solving systems of Nonlinear
     Algebraic Equations" (Chap 7, p 115-161), in `Numerical Methods for
     Nonlinear Algebraic Equations', P. Rabinowitz, editor.  Gordon and
     Breach, 1970.

The following papers are also relevant to the algorithms described in
this section,

     J.J. More', M.Y. Cosnard, "Numerical Solution of Nonlinear
     Equations", `ACM Transactions on Mathematical Software', Vol 5, No
     1, (1979), p 64-85

     C.G. Broyden, "A Class of Methods for Solving Nonlinear
     Simultaneous Equations", `Mathematics of Computation', Vol 19
     (1965), p 577-593

     J.J. More', B.S. Garbow, K.E. Hillstrom, "Testing Unconstrained
     Optimization Software", ACM Transactions on Mathematical Software,
     Vol 7, No 1 (1981), p 17-41


File: gsl-ref.info,  Node: Multidimensional Minimization,  Next: Least-Squares Fitting,  Prev: Multidimensional Root-Finding,  Up: Top

Multidimensional Minimization
*****************************

This chapter describes routines for finding minima of arbitrary
multidimensional functions.  The library provides low level components
for a variety of iterative minimizers and convergence tests.  These can
be combined by the user to achieve the desired solution, while providing
full access to the intermediate steps of the algorithms.  Each class of
methods uses the same framework, so that you can switch between
minimizers at runtime without needing to recompile your program.  Each
instance of a minimizer keeps track of its own state, allowing the
minimizers to be used in multi-threaded programs. The minimization
algorithms can be used to maximize a function by inverting its sign.

   The header file `gsl_multimin.h' contains prototypes for the
minimization functions and related declarations.

* Menu:

* Multimin Overview::
* Multimin Caveats::
* Initializing the Multidimensional Minimizer::
* Providing a function to minimize::
* Multimin Iteration::
* Multimin Stopping Criteria::
* Multimin Algorithms::
* Multimin Examples::
* Multimin References and Further Reading::


File: gsl-ref.info,  Node: Multimin Overview,  Next: Multimin Caveats,  Up: Multidimensional Minimization

Overview
========

The problem of multidimensional minimization requires finding a point x
such that the scalar function,

     f(x_1, ..., x_n)

takes a value which is lower than at any neighboring point. For smooth
functions the gradient g = \nabla f vanishes at the minimum. In general
there are no bracketing methods available for the minimization of
n-dimensional functions.  All algorithms proceed from an initial guess
using a search algorithm which attempts to move in a downhill direction.

   All algorithms making use of the gradient of the function perform a
one-dimensional line minimisation along this direction until the lowest
point is found to a suitable tolerance.  The search direction is then
updated with local information from the function and its derivatives,
and the whole process repeated until the true n-dimensional minimum is
found.

   The Nelder-Mead Simplex algorithm applies a different strategy.  It
maintains n+1 trial parameter vectors as the vertices of a
n-dimensional simplex.  In each iteration step it tries to improve the
worst vertex by a simple geometrical transformation until the size of
the simplex falls below a given tolerance.

   Several minimization algorithms are available within a single
framework. The user provides a high-level driver for the algorithms, and
the library provides the individual functions necessary for each of the
steps.  There are three main phases of the iteration.  The steps are,

   * initialize minimizer state, S, for algorithm T

   * update S using the iteration T

   * test S for convergence, and repeat iteration if necessary

Each iteration step consists either of an improvement to the
line-minimisation in the current direction or an update to the search
direction itself.  The state for the minimizers is held in a
`gsl_multimin_fdfminimizer' struct or a `gsl_multimin_fminimizer'
struct.


File: gsl-ref.info,  Node: Multimin Caveats,  Next: Initializing the Multidimensional Minimizer,  Prev: Multimin Overview,  Up: Multidimensional Minimization

Caveats
=======

Note that the minimization algorithms can only search for one local
minimum at a time.  When there are several local minima in the search
area, the first minimum to be found will be returned; however it is
difficult to predict which of the minima this will be.  In most cases,
no error will be reported if you try to find a local minimum in an area
where there is more than one.

   It is also important to note that the minimization algorithms find
local minima; there is no way to determine whether a minimum is a global
minimum of the function in question.


File: gsl-ref.info,  Node: Initializing the Multidimensional Minimizer,  Next: Providing a function to minimize,  Prev: Multimin Caveats,  Up: Multidimensional Minimization

Initializing the Multidimensional Minimizer
===========================================

The following function initializes a multidimensional minimizer.  The
minimizer itself depends only on the dimension of the problem and the
algorithm and can be reused for different problems.

 - Function: gsl_multimin_fdfminimizer *
          gsl_multimin_fdfminimizer_alloc (const
          gsl_multimin_fdfminimizer_type *T, size_t N)
 - Function: gsl_multimin_fminimizer * gsl_multimin_fminimizer_alloc
          (const gsl_multimin_fminimizer_type *T, size_t N)
     This function returns a pointer to a newly allocated instance of a
     minimizer of type T for an N-dimension function.  If there is
     insufficient memory to create the minimizer then the function
     returns a null pointer and the error handler is invoked with an
     error code of `GSL_ENOMEM'.

 - Function: int gsl_multimin_fdfminimizer_set
          (gsl_multimin_fdfminimizer * S, gsl_multimin_function_fdf
          *FDF, const gsl_vector * X, double STEP_SIZE, double TOL)
     This function initializes the minimizer S to minimize the function
     FDF starting from the initial point X.  The size of the first
     trial step is given by STEP_SIZE.  The accuracy of the line
     minimization is specified by TOL.  The precise meaning of this
     parameter depends on the method used.  Typically the line
     minimization is considered successful if the gradient of the
     function g is orthogonal to the current search direction p to a
     relative accuracy of TOL, where dot(p,g) < tol |p| |g|.

 - Function: int gsl_multimin_fminimizer_set (gsl_multimin_fminimizer *
          S, gsl_multimin_function *F, const gsl_vector * X, const
          gsl_vector * STEP_SIZE)
     This function initializes the minimizer S to minimize the function
     F, starting from the initial point X. The size of the initial
     trial steps is given in vector STEP_SIZE. The precise meaning of
     this parameter depends on the method used.

 - Function: void gsl_multimin_fdfminimizer_free
          (gsl_multimin_fdfminimizer *S)
 - Function: void gsl_multimin_fminimizer_free (gsl_multimin_fminimizer
          *S)
     This function frees all the memory associated with the minimizer S.

 - Function: const char * gsl_multimin_fdfminimizer_name (const
          gsl_multimin_fdfminimizer * S)
 - Function: const char * gsl_multimin_fminimizer_name (const
          gsl_multimin_fminimizer * S)
     This function returns a pointer to the name of the minimizer.  For
     example,

          printf ("s is a '%s' minimizer\n",
                  gsl_multimin_fdfminimizer_name (s));

     would print something like `s is a 'conjugate_pr' minimizer'.


File: gsl-ref.info,  Node: Providing a function to minimize,  Next: Multimin Iteration,  Prev: Initializing the Multidimensional Minimizer,  Up: Multidimensional Minimization

Providing a function to minimize
================================

You must provide a parametric function of n variables for the
minimizers to operate on.  You may also need to provide a routine which
calculates the gradient of the function and a third routine which
calculates both the function value and the gradient together.  In order
to allow for general parameters the functions are defined by the
following data type:

 - Data Type: gsl_multimin_function_fdf
     This data type defines a general function of n variables with
     parameters and the corresponding gradient vector of derivatives,

    `double (* f) (const gsl_vector * X, void * PARAMS)'
          this function should return the result f(x,params) for
          argument X and parameters PARAMS.

    `void (* df) (const gsl_vector * X, void * PARAMS, gsl_vector * G)'
          this function should store the N-dimensional gradient g_i = d
          f(x,params) / d x_i in the vector G for argument X and
          parameters PARAMS, returning an appropriate error code if the
          function cannot be computed.

    `void (* fdf) (const gsl_vector * X, void * PARAMS, double * f, gsl_vector * G)'
          This function should set the values of the F and G as above,
          for arguments X and parameters PARAMS.  This function
          provides an optimization of the separate functions for f(x)
          and g(x)--it is always faster to compute the function and its
          derivative at the same time.

    `size_t n'
          the dimension of the system, i.e. the number of components of
          the vectors X.

    `void * params'
          a pointer to the parameters of the function.

 - Data Type: gsl_multimin_function
     This data type defines a general function of n variables with
     parameters,

    `double (* f) (const gsl_vector * X, void * PARAMS)'
          this function should return the result f(x,params) for
          argument X and parameters PARAMS.

    `size_t n'
          the dimension of the system, i.e. the number of components of
          the vectors X.

    `void * params'
          a pointer to the parameters of the function.

The following example function defines a simple paraboloid with two
parameters,

     /* Paraboloid centered on (dp[0],dp[1]) */
     
     double
     my_f (const gsl_vector *v, void *params)
     {
       double x, y;
       double *dp = (double *)params;
     
       x = gsl_vector_get(v, 0);
       y = gsl_vector_get(v, 1);
     
       return 10.0 * (x - dp[0]) * (x - dp[0]) +
                20.0 * (y - dp[1]) * (y - dp[1]) + 30.0;
     }
     
     /* The gradient of f, df = (df/dx, df/dy). */
     void
     my_df (const gsl_vector *v, void *params,
            gsl_vector *df)
     {
       double x, y;
       double *dp = (double *)params;
     
       x = gsl_vector_get(v, 0);
       y = gsl_vector_get(v, 1);
     
       gsl_vector_set(df, 0, 20.0 * (x - dp[0]));
       gsl_vector_set(df, 1, 40.0 * (y - dp[1]));
     }
     
     /* Compute both f and df together. */
     void
     my_fdf (const gsl_vector *x, void *params,
             double *f, gsl_vector *df)
     {
       *f = my_f(x, params);
       my_df(x, params, df);
     }

The function can be initialized using the following code,

     gsl_multimin_function_fdf my_func;
     
     double p[2] = { 1.0, 2.0 }; /* center at (1,2) */
     
     my_func.f = &my_f;
     my_func.df = &my_df;
     my_func.fdf = &my_fdf;
     my_func.n = 2;
     my_func.params = (void *)p;


File: gsl-ref.info,  Node: Multimin Iteration,  Next: Multimin Stopping Criteria,  Prev: Providing a function to minimize,  Up: Multidimensional Minimization

Iteration
=========

The following function drives the iteration of each algorithm.  The
function performs one iteration to update the state of the minimizer.
The same function works for all minimizers so that different methods can
be substituted at runtime without modifications to the code.

 - Function: int gsl_multimin_fdfminimizer_iterate
          (gsl_multimin_fdfminimizer *S)
 - Function: int gsl_multimin_fminimizer_iterate
          (gsl_multimin_fminimizer *S)
     These functions perform a single iteration of the minimizer S.  If
     the iteration encounters an unexpected problem then an error code
     will be returned.

The minimizer maintains a current best estimate of the minimum at all
times.  This information can be accessed with the following auxiliary
functions,

 - Function: gsl_vector * gsl_multimin_fdfminimizer_x (const
          gsl_multimin_fdfminimizer * S)
 - Function: gsl_vector * gsl_multimin_fminimizer_x (const
          gsl_multimin_fminimizer * S)
 - Function: double gsl_multimin_fdfminimizer_minimum (const
          gsl_multimin_fdfminimizer * S)
 - Function: double gsl_multimin_fminimizer_minimum (const
          gsl_multimin_fminimizer * S)
 - Function: gsl_vector * gsl_multimin_fdfminimizer_gradient (const
          gsl_multimin_fdfminimizer * S)
 - Function: double gsl_multimin_fminimizer_size (const
          gsl_multimin_fminimizer * S)
     These functions return the current best estimate of the location
     of the minimum, the value of the function at that point, its
     gradient, and minimizer specific characteristic size for the
     minimizer S.

 - Function: int gsl_multimin_fdfminimizer_restart
          (gsl_multimin_fdfminimizer *S)
     This function resets the minimizer S to use the current point as a
     new starting point.


File: gsl-ref.info,  Node: Multimin Stopping Criteria,  Next: Multimin Algorithms,  Prev: Multimin Iteration,  Up: Multidimensional Minimization

Stopping Criteria
=================

A minimization procedure should stop when one of the following
conditions is true:

   * A minimum has been found to within the user-specified precision.

   * A user-specified maximum number of iterations has been reached.

   * An error has occurred.

The handling of these conditions is under user control.  The functions
below allow the user to test the precision of the current result.

 - Function: int gsl_multimin_test_gradient (const gsl_vector * G,
          double EPSABS)
     This function tests the norm of the gradient G against the
     absolute tolerance EPSABS. The gradient of a multidimensional
     function goes to zero at a minimum. The test returns `GSL_SUCCESS'
     if the following condition is achieved,

          |g| < epsabs

     and returns `GSL_CONTINUE' otherwise.  A suitable choice of EPSABS
     can be made from the desired accuracy in the function for small
     variations in x.  The relationship between these quantities is
     given by \delta f = g \delta x.

 - Function: int gsl_multimin_test_size (const double SIZE, double
          EPSABS)
     This function tests the minimizer specific characteristic size (if
     applicable to the used minimizer) against absolute tolerance
     EPSABS.  The test returns `GSL_SUCCESS' if the size is smaller
     than tolerance, otherwise `GSL_CONTINUE' is returned.


File: gsl-ref.info,  Node: Multimin Algorithms,  Next: Multimin Examples,  Prev: Multimin Stopping Criteria,  Up: Multidimensional Minimization

Algorithms
==========

There are several minimization methods available. The best choice of
algorithm depends on the problem.  All of the algorithms use the value
of the function and most of its gradient at each evaluation point, too.

 - Minimizer: gsl_multimin_fdfminimizer_conjugate_fr
     This is the Fletcher-Reeves conjugate gradient algorithm. The
     conjugate gradient algorithm proceeds as a succession of line
     minimizations. The sequence of search directions is used to build
     up an approximation to the curvature of the function in the
     neighborhood of the minimum.  An initial search direction P is
     chosen using the gradient and line minimization is carried out in
     that direction.  The accuracy of the line minimization is
     specified by the parameter TOL.  At the minimum along this line
     the function gradient G and the search direction P are orthogonal.
     The line minimization terminates when dot(p,g) < tol |p| |g|.  The
     search direction is updated  using the Fletcher-Reeves formula p'
     = g' - \beta g where \beta=-|g'|^2/|g|^2, and the line
     minimization is then repeated for the new search direction.

 - Minimizer: gsl_multimin_fdfminimizer_conjugate_pr
     This is the Polak-Ribiere conjugate gradient algorithm.  It is
     similar to the Fletcher-Reeves method, differing only in the
     choice of the coefficient \beta. Both methods work well when the
     evaluation point is close enough to the minimum of the objective
     function that it is well approximated by a quadratic hypersurface.

 - Minimizer: gsl_multimin_fdfminimizer_vector_bfgs
     This is the vector Broyden-Fletcher-Goldfarb-Shanno (BFGS)
     conjugate gradient algorithm.  It is a quasi-Newton method which
     builds up an approximation to the second derivatives of the
     function f using the difference between successive gradient
     vectors.  By combining the first and second derivatives the
     algorithm is able to take Newton-type steps towards the function
     minimum, assuming quadratic behavior in that region.

 - Minimizer: gsl_multimin_fdfminimizer_steepest_descent
     The steepest descent algorithm follows the downhill gradient of the
     function at each step. When a downhill step is successful the
     step-size is increased by a factor of two.  If the downhill step
     leads to a higher function value then the algorithm backtracks and
     the step size is decreased using the parameter TOL.  A suitable
     value of TOL for most applications is 0.1.  The steepest descent
     method is inefficient and is included only for demonstration
     purposes.

 - Minimizer: gsl_multimin_fminimizer_nmsimplex
     This is the Simplex algorithm of Nelder and Mead. It constructs n
     vectors p_i from the starting vector X and the vector STEP_SIZE as
     follows:
          p_0 = (x_0, x_1, ... , x_n)
          p_1 = (x_0 + step_size_0, x_1, ... , x_n)
          p_2 = (x_0, x_1 + step_size_1, ... , x_n)
          ... = ...
          p_n = (x_0, x_1, ... , x_n+step_size_n)

     These vectors form the n+1 vertices of a simplex in n dimensions.
     On each iteration step the algorithm tries to improve the
     parameter vector p_i corresponding to the highest function value
     by simple geometrical transformations.  These are reflection,
     reflection followed by expansion, contraction and multiple
     contraction. Using these transformations the simplex moves through
     the parameter space towards the minimum, where it contracts itself.

     After each iteration, the best vertex is returned.  Note, that due
     to the nature of the algorithm not every step improves the current
     best parameter vector.  Usually several iterations are required.

     The routine calculates the minimizer specific characteristic size
     as the average distance from the geometrical center of the simplex
     to all its vertices.  This size can be used as a stopping
     criteria, as the simplex contracts itself near the minimum. The
     size is returned by the function `gsl_multimin_fminimizer_size'.


File: gsl-ref.info,  Node: Multimin Examples,  Next: Multimin References and Further Reading,  Prev: Multimin Algorithms,  Up: Multidimensional Minimization

Examples
========

This example program finds the minimum of the paraboloid function
defined earlier.  The location of the minimum is offset from the origin
in x and y, and the function value at the minimum is non-zero. The main
program is given below, it requires the example function given earlier
in this chapter.

     int
     main (void)
     {
       size_t iter = 0;
       int status;
     
       const gsl_multimin_fdfminimizer_type *T;
       gsl_multimin_fdfminimizer *s;
     
       /* Position of the minimum (1,2). */
       double par[2] = { 1.0, 2.0 };
     
       gsl_vector *x;
       gsl_multimin_function_fdf my_func;
     
       my_func.f = &my_f;
       my_func.df = &my_df;
       my_func.fdf = &my_fdf;
       my_func.n = 2;
       my_func.params = &par;
     
       /* Starting point, x = (5,7) */
       x = gsl_vector_alloc (2);
       gsl_vector_set (x, 0, 5.0);
       gsl_vector_set (x, 1, 7.0);
     
       T = gsl_multimin_fdfminimizer_conjugate_fr;
       s = gsl_multimin_fdfminimizer_alloc (T, 2);
     
       gsl_multimin_fdfminimizer_set (s, &my_func, x, 0.01, 1e-4);
     
       do
         {
           iter++;
           status = gsl_multimin_fdfminimizer_iterate (s);
     
           if (status)
             break;
     
           status = gsl_multimin_test_gradient (s->gradient, 1e-3);
     
           if (status == GSL_SUCCESS)
             printf ("Minimum found at:\n");
     
           printf ("%5d %.5f %.5f %10.5f\n", iter,
                   gsl_vector_get (s->x, 0),
                   gsl_vector_get (s->x, 1),
                   s->f);
     
         }
       while (status == GSL_CONTINUE && iter < 100);
     
       gsl_multimin_fdfminimizer_free (s);
       gsl_vector_free (x);
     
       return 0;
     }

The initial step-size is chosen as 0.01, a conservative estimate in this
case, and the line minimization parameter is set at 0.0001.  The program
terminates when the norm of the gradient has been reduced below 0.001.
The output of the program is shown below,

              x       y         f
         1 4.99629 6.99072  687.84780
         2 4.98886 6.97215  683.55456
         3 4.97400 6.93501  675.01278
         4 4.94429 6.86073  658.10798
         5 4.88487 6.71217  625.01340
         6 4.76602 6.41506  561.68440
         7 4.52833 5.82083  446.46694
         8 4.05295 4.63238  261.79422
         9 3.10219 2.25548   75.49762
        10 2.85185 1.62963   67.03704
        11 2.19088 1.76182   45.31640
        12 0.86892 2.02622   30.18555
     Minimum found at:
        13 1.00000 2.00000   30.00000

Note that the algorithm gradually increases the step size as it
successfully moves downhill, as can be seen by plotting the successive
points.

The conjugate gradient algorithm finds the minimum on its second
direction because the function is purely quadratic. Additional
iterations would be needed for a more complicated function.

   Here is another example using the Nelder-Mead Simplex algorithm to
minimize the same example object function, as above.

     int
     main(void)
     {
       size_t np = 2;
       double par[2] = {1.0, 2.0};
     
       const gsl_multimin_fminimizer_type *T =
         gsl_multimin_fminimizer_nmsimplex;
       gsl_multimin_fminimizer *s = NULL;
       gsl_vector *ss, *x;
       gsl_multimin_function minex_func;
     
       size_t iter = 0, i;
       int status;
       double size;
     
       /* Initial vertex size vector */
       ss = gsl_vector_alloc (np);
     
       /* Set all step sizes to 1 */
       gsl_vector_set_all (ss, 1.0);
     
       /* Starting point */
       x = gsl_vector_alloc (np);
     
       gsl_vector_set (x, 0, 5.0);
       gsl_vector_set (x, 1, 7.0);
     
       /* Initialize method and iterate */
       minex_func.f = &my_f;
       minex_func.n = np;
       minex_func.params = (void *)&par;
     
       s = gsl_multimin_fminimizer_alloc (T, np);
       gsl_multimin_fminimizer_set (s, &minex_func, x, ss);
     
       do
         {
           iter++;
           status = gsl_multimin_fminimizer_iterate(s);
     
           if (status)
             break;
     
           size = gsl_multimin_fminimizer_size (s);
           status = gsl_multimin_test_size (size, 1e-2);
     
           if (status == GSL_SUCCESS)
             {
               printf ("converged to minimum at\n");
             }
     
           printf ("%5d ", iter);
           for (i = 0; i < np; i++)
             {
               printf ("%10.3e ", gsl_vector_get (s->x, i));
             }
           printf ("f() = %7.3f size = %.3f\n", s->fval, size);
         }
       while (status == GSL_CONTINUE && iter < 100);
     
       gsl_vector_free(x);
       gsl_vector_free(ss);
       gsl_multimin_fminimizer_free (s);
     
       return status;
     }

The minimum search stops when the Simplex size drops to 0.01. The
output is shown below.

         1  6.500e+00  5.000e+00 f() = 512.500 size = 1.082
         2  5.250e+00  4.000e+00 f() = 290.625 size = 1.372
         3  5.250e+00  4.000e+00 f() = 290.625 size = 1.372
         4  5.500e+00  1.000e+00 f() = 252.500 size = 1.372
         5  2.625e+00  3.500e+00 f() = 101.406 size = 1.823
         6  3.469e+00  1.375e+00 f() = 98.760  size = 1.526
         7  1.820e+00  3.156e+00 f() = 63.467  size = 1.105
         8  1.820e+00  3.156e+00 f() = 63.467  size = 1.105
         9  1.016e+00  2.812e+00 f() = 43.206  size = 1.105
        10  2.041e+00  2.008e+00 f() = 40.838  size = 0.645
        11  1.236e+00  1.664e+00 f() = 32.816  size = 0.645
        12  1.236e+00  1.664e+00 f() = 32.816  size = 0.447
        13  5.225e-01  1.980e+00 f() = 32.288  size = 0.447
        14  1.103e+00  2.073e+00 f() = 30.214  size = 0.345
        15  1.103e+00  2.073e+00 f() = 30.214  size = 0.264
        16  1.103e+00  2.073e+00 f() = 30.214  size = 0.160
        17  9.864e-01  1.934e+00 f() = 30.090  size = 0.132
        18  9.190e-01  1.987e+00 f() = 30.069  size = 0.092
        19  1.028e+00  2.017e+00 f() = 30.013  size = 0.056
        20  1.028e+00  2.017e+00 f() = 30.013  size = 0.046
        21  1.028e+00  2.017e+00 f() = 30.013  size = 0.033
        22  9.874e-01  1.985e+00 f() = 30.006  size = 0.028
        23  9.846e-01  1.995e+00 f() = 30.003  size = 0.023
        24  1.007e+00  2.003e+00 f() = 30.001  size = 0.012
     converged to minimum at
        25  1.007e+00  2.003e+00 f() = 30.001  size = 0.010

The simplex size first increases, while the simplex moves towards the
minimum. After a while the size begins to decrease as the simplex
contracts around the minimum.


File: gsl-ref.info,  Node: Multimin References and Further Reading,  Prev: Multimin Examples,  Up: Multidimensional Minimization

References and Further Reading
==============================

A brief description of multidimensional minimization algorithms and
further references can be found in the following book,

     C.W. Ueberhuber, `Numerical Computation (Volume 2)', Chapter 14,
     Section 4.4 "Minimization Methods", p. 325--335, Springer (1997),
     ISBN 3-540-62057-5.

     J.A. Nelder and R. Mead, `A simplex method for function
     minimization', Computer Journal vol. 7 (1965), 308--315.


File: gsl-ref.info,  Node: Least-Squares Fitting,  Next: Nonlinear Least-Squares Fitting,  Prev: Multidimensional Minimization,  Up: Top

Least-Squares Fitting
*********************

This chapter describes routines for performing least squares fits to
experimental data using linear combinations of functions. The data may
be weighted or unweighted.  For weighted data the functions compute the
best fit parameters and their associated covariance matrix.  For
unweighted data the covariance matrix is estimated from the scatter of
the points, giving a variance-covariance matrix. The functions are
divided into separate versions for simple one- or two-parameter
regression and multiple-parameter fits.  The functions are declared in
the header file `gsl_fit.h'.

* Menu:

* Linear regression::
* Linear fitting without a constant term::
* Multi-parameter fitting::
* Fitting Examples::
* Fitting References and Further Reading::


File: gsl-ref.info,  Node: Linear regression,  Next: Linear fitting without a constant term,  Up: Least-Squares Fitting

Linear regression
=================

The functions described in this section can be used to perform
least-squares fits to a straight line model, Y = c_0 + c_1 X.  For
weighted data the best-fit is found by minimizing the weighted sum of
squared residuals, \chi^2,

     \chi^2 = \sum_i w_i (y_i - (c_0 + c_1 x_i))^2

for the parameters c_0, c_1.  For unweighted data the sum is computed
with w_i = 1.

 - Function: int gsl_fit_linear (const double * X, const size_t
          XSTRIDE, const double * Y, const size_t YSTRIDE, size_t N,
          double * C0, double * C1, double * COV00, double * COV01,
          double * COV11, double * SUMSQ)
     This function computes the best-fit linear regression coefficients
     (C0,C1) of the model Y = c_0 + c_1 X for the datasets (X, Y), two
     vectors of length N with strides XSTRIDE and YSTRIDE.  The
     variance-covariance matrix for the parameters (C0, C1) is
     estimated from the scatter of the points around the best-fit line
     and returned via the parameters (COV00, COV01, COV11).  The sum of
     squares of the residuals from the best-fit line is returned in
     SUMSQ.

 - Function: int gsl_fit_wlinear (const double * X, const size_t
          XSTRIDE, const double * W, const size_t WSTRIDE, const double
          * Y, const size_t YSTRIDE, size_t N, double * C0, double *
          C1, double * COV00, double * COV01, double * COV11, double *
          CHISQ)
     This function computes the best-fit linear regression coefficients
     (C0,C1) of the model Y = c_0 + c_1 X for the weighted datasets (X,
     Y), two vectors of length N with strides XSTRIDE and YSTRIDE.  The
     vector W, of length N and stride WSTRIDE, specifies the weight of
     each datapoint. The weight is the reciprocal of the variance for
     each datapoint in Y.

     The covariance matrix for the parameters (C0, C1) is estimated
     from weighted data and returned via the parameters (COV00, COV01,
     COV11).  The weighted sum of squares of the residuals from the
     best-fit line, \chi^2, is returned in CHISQ.

 - Function: int gsl_fit_linear_est (double X, double C0, double C1,
          double C00, double C01, double C11, double *Y, double *Y_ERR)
     This function uses the best-fit linear regression coefficients
     C0,C1 and their estimated covariance COV00,COV01,COV11 to compute
     the fitted function Y and its standard deviation Y_ERR for the
     model Y = c_0 + c_1 X at the point X.


File: gsl-ref.info,  Node: Linear fitting without a constant term,  Next: Multi-parameter fitting,  Prev: Linear regression,  Up: Least-Squares Fitting

Linear fitting without a constant term
======================================

The functions described in this section can be used to perform
least-squares fits to a straight line model without a constant term, Y
= c_1 X.  For weighted data the best-fit is found by minimizing the
weighted sum of squared residuals, \chi^2,

     \chi^2 = \sum_i w_i (y_i - c_1 x_i)^2

for the parameter c_1.  For unweighted data the sum is computed with
w_i = 1.

 - Function: int gsl_fit_mul (const double * X, const size_t XSTRIDE,
          const double * Y, const size_t YSTRIDE, size_t N, double *
          C1, double * COV11, double * SUMSQ)
     This function computes the best-fit linear regression coefficient
     C1 of the model Y = c_1 X for the datasets (X, Y), two vectors of
     length N with strides XSTRIDE and YSTRIDE.  The variance of the
     parameter C1 is estimated from the scatter of the points around
     the best-fit line and returned via the parameter COV11.  The sum
     of squares of the residuals from the best-fit line is returned in
     SUMSQ.

 - Function: int gsl_fit_wmul (const double * X, const size_t XSTRIDE,
          const double * W, const size_t WSTRIDE, const double * Y,
          const size_t YSTRIDE, size_t N, double * C1, double * COV11,
          double * SUMSQ)
     This function computes the best-fit linear regression coefficient
     C1 of the model Y = c_1 X for the weighted datasets (X, Y), two
     vectors of length N with strides XSTRIDE and YSTRIDE.  The vector
     W, of length N and stride WSTRIDE, specifies the weight of each
     datapoint. The weight is the reciprocal of the variance for each
     datapoint in Y.

     The variance of the parameter C1 is estimated from the weighted
     data and returned via the parameters COV11.  The weighted sum of
     squares of the residuals from the best-fit line, \chi^2, is
     returned in CHISQ.

 - Function: int gsl_fit_mul_est (double X, double C1, double C11,
          double *Y, double *Y_ERR)
     This function uses the best-fit linear regression coefficient C1
     and its estimated covariance COV11 to compute the fitted function
     Y and its standard deviation Y_ERR for the model Y = c_1 X at the
     point X.


File: gsl-ref.info,  Node: Multi-parameter fitting,  Next: Fitting Examples,  Prev: Linear fitting without a constant term,  Up: Least-Squares Fitting

Multi-parameter fitting
=======================

The functions described in this section perform least-squares fits to a
general linear model, y = X c where y is a vector of n observations, X
is an n by p matrix of predictor variables, and c are the p unknown
best-fit parameters, which are to be estimated.

   The best-fit is found by minimizing the weighted sums of squared
residuals, \chi^2,

     \chi^2 = (y - X c)^T W (y - X c)

with respect to the parameters c. The weights are specified by the
diagonal elements of the n by n matrix W.  For unweighted data W is
replaced by the identity matrix.

   This formulation can be used for fits to any number of functions
and/or variables by preparing the n-by-p matrix X appropriately.  For
example, to fit to a p-th order polynomial in X, use the following
matrix,

     X_{ij} = x_i^j

where the index i runs over the observations and the index j runs from
0 to p-1.

   To fit to a set of p sinusoidal functions with fixed frequencies
\omega_1, \omega_2, ..., \omega_p, use,

     X_{ij} = sin(\omega_j x_i)

To fit to p independent variables x_1, x_2, ..., x_p, use,

     X_{ij} = x_j(i)

where x_j(i) is the i-th value of the predictor variable x_j.

   The functions described in this section are declared in the header
file `gsl_multifit.h'.

   The solution of the general linear least-squares system requires an
additional working space for intermediate results, such as the singular
value decomposition of the matrix X.

 - Function: gsl_multifit_linear_workspace * gsl_multifit_linear_alloc
          (size_t N, size_t P)
     This function allocates a workspace for fitting a model to N
     observations using P parameters.

 - Function: void gsl_multifit_linear_free
          (gsl_multifit_linear_workspace * WORK)
     This function frees the memory associated with the workspace W.

 - Function: int gsl_multifit_linear (const gsl_matrix * X, const
          gsl_vector * Y, gsl_vector * C, gsl_matrix * COV, double *
          CHISQ, gsl_multifit_linear_workspace * WORK)
     This function computes the best-fit parameters C of the model y =
     X c for the observations Y and the matrix of predictor variables
     X.  The variance-covariance matrix of the model parameters COV is
     estimated from the scatter of the observations about the best-fit.
     The sum of squares of the residuals from the best-fit, \chi^2, is
     returned in CHISQ.

     The best-fit is found by singular value decomposition of the matrix
     X using the preallocated workspace provided in WORK. The modified
     Golub-Reinsch SVD algorithm is used, with column scaling to
     improve the accuracy of the singular values. Any components which
     have zero singular value (to machine precision) are discarded from
     the fit.

 - Function: int gsl_multifit_wlinear (const gsl_matrix * X, const
          gsl_vector * W, const gsl_vector * Y, gsl_vector * C,
          gsl_matrix * COV, double * CHISQ,
          gsl_multifit_linear_workspace * WORK)
     This function computes the best-fit parameters C of the model y =
     X c for the observations Y and the matrix of predictor variables
     X.  The covariance matrix of the model parameters COV is estimated
     from the weighted data.  The weighted sum of squares of the
     residuals from the best-fit, \chi^2, is returned in CHISQ.

     The best-fit is found by singular value decomposition of the matrix
     X using the preallocated workspace provided in WORK. Any
     components which have zero singular value (to machine precision)
     are discarded from the fit.


File: gsl-ref.info,  Node: Fitting Examples,  Next: Fitting References and Further Reading,  Prev: Multi-parameter fitting,  Up: Least-Squares Fitting

Examples
========

The following program computes a least squares straight-line fit to a
simple (fictitious) dataset, and outputs the best-fit line and its
associated one standard-deviation error bars.

     #include <stdio.h>
     #include <gsl/gsl_fit.h>
     
     int
     main (void)
     {
       int i, n = 4;
       double x[4] = { 1970, 1980, 1990, 2000 };
       double y[4] = {   12,   11,   14,   13 };
       double w[4] = {  0.1,  0.2,  0.3,  0.4 };
     
       double c0, c1, cov00, cov01, cov11, chisq;
     
       gsl_fit_wlinear (x, 1, w, 1, y, 1, n, 
                        &c0, &c1, &cov00, &cov01, &cov11, 
                        &chisq);
     
       printf ("# best fit: Y = %g + %g X\n", c0, c1);
       printf ("# covariance matrix:\n");
       printf ("# [ %g, %g\n#   %g, %g]\n", 
               cov00, cov01, cov01, cov11);
       printf ("# chisq = %g\n", chisq);
     
       for (i = 0; i < n; i++)
         printf ("data: %g %g %g\n", 
                        x[i], y[i], 1/sqrt(w[i]));
     
       printf ("\n");
     
       for (i = -30; i < 130; i++)
         {
           double xf = x[0] + (i/100.0) * (x[n-1] - x[0]);
           double yf, yf_err;
     
           gsl_fit_linear_est (xf, 
                               c0, c1, 
                               cov00, cov01, cov11, 
                               &yf, &yf_err);
     
           printf ("fit: %g %g\n", xf, yf);
           printf ("hi : %g %g\n", xf, yf + yf_err);
           printf ("lo : %g %g\n", xf, yf - yf_err);
         }
       return 0;
     }

The following commands extract the data from the output of the program
and display it using the GNU plotutils `graph' utility,

     $ ./demo > tmp
     $ more tmp
     # best fit: Y = -106.6 + 0.06 X
     # covariance matrix:
     # [ 39602, -19.9
     #   -19.9, 0.01]
     # chisq = 0.8
     
     $ for n in data fit hi lo ;
        do
          grep "^$n" tmp | cut -d: -f2 > $n ;
        done
     $ graph -T X -X x -Y y -y 0 20 -m 0 -S 2 -Ie data
          -S 0 -I a -m 1 fit -m 2 hi -m 2 lo

   The next program performs a quadratic fit y = c_0 + c_1 x + c_2 x^2
to a weighted dataset using the generalised linear fitting function
`gsl_multifit_wlinear'.  The model matrix X for a quadratic fit is
given by,

     X = [ 1   , x_0  , x_0^2 ;
           1   , x_1  , x_1^2 ;
           1   , x_2  , x_2^2 ;
           ... , ...  , ...   ]

where the column of ones corresponds to the constant term c_0.  The two
remaining columns corresponds to the terms c_1 x and c_2 x^2.

   The program reads N lines of data in the format (X, Y, ERR) where
ERR is the error (standard deviation) in the value Y.

     #include <stdio.h>
     #include <gsl/gsl_multifit.h>
     
     int
     main (int argc, char **argv)
     {
       int i, n;
       double xi, yi, ei, chisq;
       gsl_matrix *X, *cov;
       gsl_vector *y, *w, *c;
     
       if (argc != 2)
         {
           fprintf (stderr,"usage: fit n < data\n");
           exit (-1);
         }
     
       n = atoi (argv[1]);
     
       X = gsl_matrix_alloc (n, 3);
       y = gsl_vector_alloc (n);
       w = gsl_vector_alloc (n);
     
       c = gsl_vector_alloc (3);
       cov = gsl_matrix_alloc (3, 3);
     
       for (i = 0; i < n; i++)
         {
           int count = fscanf (stdin, "%lg %lg %lg",
                               &xi, &yi, &ei);
     
           if (count != 3)
             {
               fprintf (stderr, "error reading file\n");
               exit (-1);
             }
     
           printf ("%g %g +/- %g\n", xi, yi, ei);
           
           gsl_matrix_set (X, i, 0, 1.0);
           gsl_matrix_set (X, i, 1, xi);
           gsl_matrix_set (X, i, 2, xi*xi);
           
           gsl_vector_set (y, i, yi);
           gsl_vector_set (w, i, 1.0/(ei*ei));
         }
     
       {
         gsl_multifit_linear_workspace * work 
           = gsl_multifit_linear_alloc (n, 3);
         gsl_multifit_wlinear (X, w, y, c, cov,
                               &chisq, work);
         gsl_multifit_linear_free (work);
       }
     
     #define C(i) (gsl_vector_get(c,(i)))
     #define COV(i,j) (gsl_matrix_get(cov,(i),(j)))
     
       {
         printf ("# best fit: Y = %g + %g X + %g X^2\n", 
                 C(0), C(1), C(2));
     
         printf ("# covariance matrix:\n");
         printf ("[ %+.5e, %+.5e, %+.5e  \n",
                    COV(0,0), COV(0,1), COV(0,2));
         printf ("  %+.5e, %+.5e, %+.5e  \n", 
                    COV(1,0), COV(1,1), COV(1,2));
         printf ("  %+.5e, %+.5e, %+.5e ]\n", 
                    COV(2,0), COV(2,1), COV(2,2));
         printf ("# chisq = %g\n", chisq);
       }
       return 0;
     }

A suitable set of data for fitting can be generated using the following
program.  It outputs a set of points with gaussian errors from the curve
y = e^x in the region 0 < x < 2.

     #include <stdio.h>
     #include <math.h>
     #include <gsl/gsl_randist.h>
     
     int
     main (void)
     {
       double x;
       const gsl_rng_type * T;
       gsl_rng * r;
       
       gsl_rng_env_setup ();
       
       T = gsl_rng_default;
       r = gsl_rng_alloc (T);
     
       for (x = 0.1; x < 2; x+= 0.1)
         {
           double y0 = exp (x);
           double sigma = 0.1 * y0;
           double dy = gsl_ran_gaussian (r, sigma);
     
           printf ("%g %g %g\n", x, y0 + dy, sigma);
         }
       return 0;
     }

The data can be prepared by running the resulting executable program,

     $ ./generate > exp.dat
     $ more exp.dat
     0.1 0.97935 0.110517
     0.2 1.3359 0.12214
     0.3 1.52573 0.134986
     0.4 1.60318 0.149182
     0.5 1.81731 0.164872
     0.6 1.92475 0.182212
     ....

To fit the data use the previous program, with the number of data points
given as the first argument.  In this case there are 19 data points.

     $ ./fit 19 < exp.dat
     0.1 0.97935 +/- 0.110517
     0.2 1.3359 +/- 0.12214
     ...
     # best fit: Y = 1.02318 + 0.956201 X + 0.876796 X^2
     # covariance matrix:
     [ +1.25612e-02, -3.64387e-02, +1.94389e-02
       -3.64387e-02, +1.42339e-01, -8.48761e-02
       +1.94389e-02, -8.48761e-02, +5.60243e-02 ]
     # chisq = 23.0987

The parameters of the quadratic fit match the coefficients of the
expansion of e^x, taking into account the errors on the parameters and
the O(x^3) difference between the exponential and quadratic functions
for the larger values of x.  The errors on the parameters are given by
the square-root of the corresponding diagonal elements of the
covariance matrix.  The chi-squared per degree of freedom is 1.4,
indicating a reasonable fit to the data.


File: gsl-ref.info,  Node: Fitting References and Further Reading,  Prev: Fitting Examples,  Up: Least-Squares Fitting

References and Further Reading
==============================

A summary of formulas and techniques for least squares fitting can be
found in the "Statistics" chapter of the Annual Review of Particle
Physics prepared by the Particle Data Group,

     `Review of Particle Properties', R.M. Barnett et al., Physical
     Review D54, 1 (1996) `http://pdg.lbl.gov/'

The Review of Particle Physics is available online at the website given
above.

The tests used to prepare these routines are based on the NIST
Statistical Reference Datasets. The datasets and their documentation are
available from NIST at the following website,

           `http://www.nist.gov/itl/div898/strd/index.html'.


File: gsl-ref.info,  Node: Nonlinear Least-Squares Fitting,  Next: Physical Constants,  Prev: Least-Squares Fitting,  Up: Top

Nonlinear Least-Squares Fitting
*******************************

This chapter describes functions for multidimensional nonlinear
least-squares fitting.  The library provides low level components for a
variety of iterative solvers and convergence tests.  These can be
combined by the user to achieve the desired solution, with full access
to the intermediate steps of the iteration.  Each class of methods uses
the same framework, so that you can switch between solvers at runtime
without needing to recompile your program.  Each instance of a solver
keeps track of its own state, allowing the solvers to be used in
multi-threaded programs.

   The header file `gsl_multifit_nlin.h' contains prototypes for the
multidimensional nonlinear fitting functions and related declarations.

* Menu:

* Overview of Nonlinear Least-Squares Fitting::
* Initializing the Nonlinear Least-Squares Solver::
* Providing the Function to be Minimized::
* Iteration of the Minimization Algorithm::
* Search Stopping Parameters for Minimization Algorithms::
* Minimization Algorithms using Derivatives::
* Minimization Algorithms without Derivatives::
* Computing the covariance matrix of best fit parameters::
* Example programs for Nonlinear Least-Squares Fitting::
* References and Further Reading for Nonlinear Least-Squares Fitting::


File: gsl-ref.info,  Node: Overview of Nonlinear Least-Squares Fitting,  Next: Initializing the Nonlinear Least-Squares Solver,  Up: Nonlinear Least-Squares Fitting

Overview
========

The problem of multidimensional nonlinear least-squares fitting requires
the minimization of the squared residuals of n functions, f_i, in p
parameters, x_i,

     \Phi(x) = (1/2) \sum_{i=1}^{n} f_i(x_1, ..., x_p)^2
             = (1/2) || F(x) ||^2

All algorithms proceed from an initial guess using the linearization,

     \psi(p) = || F(x+p) || ~=~ || F(x) + J p ||

where x is the initial point, p is the proposed step and J is the
Jacobian matrix J_{ij} = d f_i / d x_j.  Additional strategies are used
to enlarge the region of convergence.  These include requiring a
decrease in the norm ||F|| on each step or using a trust region to
avoid steps which fall outside the linear regime.


File: gsl-ref.info,  Node: Initializing the Nonlinear Least-Squares Solver,  Next: Providing the Function to be Minimized,  Prev: Overview of Nonlinear Least-Squares Fitting,  Up: Nonlinear Least-Squares Fitting

Initializing the Solver
=======================

 - Function: gsl_multifit_fsolver * gsl_multifit_fsolver_alloc (const
          gsl_multifit_fsolver_type * T, size_t N, size_t P)
     This function returns a pointer to a newly allocated instance of a
     solver of type T for N observations and P parameters.

     If there is insufficient memory to create the solver then the
     function returns a null pointer and the error handler is invoked
     with an error code of `GSL_ENOMEM'.

 - Function: gsl_multifit_fdfsolver * gsl_multifit_fdfsolver_alloc
          (const gsl_multifit_fdfsolver_type * T, size_t N, size_t P)
     This function returns a pointer to a newly allocated instance of a
     derivative solver of type T for N observations and P parameters.
     For example, the following code creates an instance of a
     Levenberg-Marquardt solver for 100 data points and 3 parameters,

          const gsl_multifit_fdfsolver_type * T
              = gsl_multifit_fdfsolver_lmder;
          gsl_multifit_fdfsolver * s
              = gsl_multifit_fdfsolver_alloc (T, 100, 3);

     If there is insufficient memory to create the solver then the
     function returns a null pointer and the error handler is invoked
     with an error code of `GSL_ENOMEM'.

 - Function: int gsl_multifit_fsolver_set (gsl_multifit_fsolver * S,
          gsl_multifit_function * F, gsl_vector * X)
     This function initializes, or reinitializes, an existing solver S
     to use the function F and the initial guess X.

 - Function: int gsl_multifit_fdfsolver_set (gsl_multifit_fdfsolver *
          S, gsl_function_fdf * FDF, gsl_vector * X)
     This function initializes, or reinitializes, an existing solver S
     to use the function and derivative FDF and the initial guess X.

 - Function: void gsl_multifit_fsolver_free (gsl_multifit_fsolver * S)
 - Function: void gsl_multifit_fdfsolver_free (gsl_multifit_fdfsolver *
          S)
     These functions free all the memory associated with the solver S.

 - Function: const char * gsl_multifit_fsolver_name (const
          gsl_multifit_fdfsolver * S)
 - Function: const char * gsl_multifit_fdfsolver_name (const
          gsl_multifit_fdfsolver * S)
     These functions return a pointer to the name of the solver.  For
     example,

          printf ("s is a '%s' solver\n",
                  gsl_multifit_fdfsolver_name (s));

     would print something like `s is a 'lmder' solver'.


File: gsl-ref.info,  Node: Providing the Function to be Minimized,  Next: Iteration of the Minimization Algorithm,  Prev: Initializing the Nonlinear Least-Squares Solver,  Up: Nonlinear Least-Squares Fitting

Providing the Function to be Minimized
======================================

You must provide n functions of p variables for the minimization
algorithms to operate on.  In order to allow for general parameters the
functions are defined by the following data types:

 - Data Type: gsl_multifit_function
     This data type defines a general system of functions with
     parameters.

    `int (* f) (const gsl_vector * X, void * PARAMS, gsl_vector * F)'
          this function should store the vector result f(x,params) in F
          for argument X and parameters PARAMS, returning an
          appropriate error code if the function cannot be computed.

    `size_t n'
          the number of functions, i.e. the number of components of the
          vector F.

    `size_t p'
          the number of independent variables, i.e. the number of
          components of the vectors X.

    `void * params'
          a pointer to the parameters of the function.

 - Data Type: gsl_multifit_function_fdf
     This data type defines a general system of functions with
     parameters and the corresponding Jacobian matrix of derivatives,

    `int (* f) (const gsl_vector * X, void * PARAMS, gsl_vector * F)'
          this function should store the vector result f(x,params) in F
          for argument X and parameters PARAMS, returning an
          appropriate error code if the function cannot be computed.

    `int (* df) (const gsl_vector * X, void * PARAMS, gsl_matrix * J)'
          this function should store the N-by-P matrix result J_ij = d
          f_i(x,params) / d x_j in J for argument X and parameters
          PARAMS, returning an appropriate error code if the function
          cannot be computed.

    `int (* fdf) (const gsl_vector * X, void * PARAMS, gsl_vector * F, gsl_matrix * J)'
          This function should set the values of the F and J as above,
          for arguments X and parameters PARAMS.  This function
          provides an optimization of the separate functions for f(x)
          and J(x)--it is always faster to compute the function and its
          derivative at the same time.

    `size_t n'
          the number of functions, i.e. the number of components of the
          vector F.

    `size_t p'
          the number of independent variables, i.e. the number of
          components of the vectors X.

    `void * params'
          a pointer to the parameters of the function.


File: gsl-ref.info,  Node: Iteration of the Minimization Algorithm,  Next: Search Stopping Parameters for Minimization Algorithms,  Prev: Providing the Function to be Minimized,  Up: Nonlinear Least-Squares Fitting

Iteration
=========

The following functions drive the iteration of each algorithm.  Each
function performs one iteration to update the state of any solver of the
corresponding type.  The same functions work for all solvers so that
different methods can be substituted at runtime without modifications to
the code.

 - Function: int gsl_multifit_fsolver_iterate (gsl_multifit_fsolver * S)
 - Function: int gsl_multifit_fdfsolver_iterate (gsl_multifit_fdfsolver
          * S)
     These functions perform a single iteration of the solver S.  If
     the iteration encounters an unexpected problem then an error code
     will be returned.  The solver maintains a current estimate of the
     best-fit parameters at all times. This information can be accessed
     with the following auxiliary functions,

   The solver struct S contains the following entries, which can be
used to track the progress of the solution:

`gsl_vector * x'
     The current position.

`gsl_vector * f'
     The function value at the current position.

`gsl_vector * dx'
     The difference between the current position and the previous
     position, i.e. the last step, taken as a vector.

`gsl_matrix * J'
     The Jacobian matrix at the current position (for the
     `gsl_multifit_fdfsolver' struct only)

 - Function: gsl_vector * gsl_multifit_fsolver_position (const
          gsl_multifit_fsolver * S)
 - Function: gsl_vector * gsl_multifit_fdfsolver_position (const
          gsl_multifit_fdfsolver * S)
     These functions return the current position (i.e. best-fit
     parameters) `s->x' of the solver S.


File: gsl-ref.info,  Node: Search Stopping Parameters for Minimization Algorithms,  Next: Minimization Algorithms using Derivatives,  Prev: Iteration of the Minimization Algorithm,  Up: Nonlinear Least-Squares Fitting

Search Stopping Parameters
==========================

A minimization procedure should stop when one of the following
conditions is true:

   * A minimum has been found to within the user-specified precision.

   * A user-specified maximum number of iterations has been reached.

   * An error has occurred.

The handling of these conditions is under user control.  The functions
below allow the user to test the current estimate of the best-fit
parameters in several standard ways.

 - Function: int gsl_multifit_test_delta (const gsl_vector * DX, const
          gsl_vector * X, double EPSABS, double EPSREL)
     This function tests for the convergence of the sequence by
     comparing the last step DX with the absolute error EPSABS and
     relative error EPSREL to the current position X.  The test returns
     `GSL_SUCCESS' if the following condition is achieved,

          |dx_i| < epsabs + epsrel |x_i|

     for each component of X and returns `GSL_CONTINUE' otherwise.

 - Function: int gsl_multifit_test_gradient (const gsl_vector * G,
          double EPSABS)
     This function tests the residual gradient G against the absolute
     error bound EPSABS.  Mathematically, the gradient should be
     exactly zero at the minimum. The test returns `GSL_SUCCESS' if the
     following condition is achieved,

          \sum_i |g_i| < epsabs

     and returns `GSL_CONTINUE' otherwise.  This criterion is suitable
     for situations where the precise location of the minimum, x, is
     unimportant provided a value can be found where the gradient is
     small enough.

 - Function: int gsl_multifit_gradient (const gsl_matrix * J, const
          gsl_vector * F, gsl_vector * G)
     This function computes the gradient G of \Phi(x) = (1/2)
     ||F(x)||^2 from the Jacobian matrix J and the function values F,
     using the formula g = J^T f.


File: gsl-ref.info,  Node: Minimization Algorithms using Derivatives,  Next: Minimization Algorithms without Derivatives,  Prev: Search Stopping Parameters for Minimization Algorithms,  Up: Nonlinear Least-Squares Fitting

Minimization Algorithms using Derivatives
=========================================

The minimization algorithms described in this section make use of both
the function and its derivative.  They require an initial guess for the
location of the minimum. There is no absolute guarantee of convergence
- the function must be suitable for this technique and the initial
guess must be sufficiently close to the minimum for it to work.

 - Derivative Solver: gsl_multifit_fdfsolver_lmsder
     This is a robust and efficient version of the Levenberg-Marquardt
     algorithm as implemented in the scaled LMDER routine in MINPACK.
     Minpack was written by Jorge J. More', Burton S. Garbow and
     Kenneth E. Hillstrom.

     The algorithm uses a generalized trust region to keep each step
     under control.  In order to be accepted a proposed new position x'
     must satisfy the condition |D (x' - x)| < \delta, where D is a
     diagonal scaling matrix and \delta is the size of the trust
     region.  The components of D are computed internally, using the
     column norms of the Jacobian to estimate the sensitivity of the
     residual to each component of x.  This improves the behavior of the
     algorithm for badly scaled functions.

     On each iteration the algorithm attempts to minimize the linear
     system |F + J p| subject to the constraint |D p| < \Delta.  The
     solution to this constrained linear system is found using the
     Levenberg-Marquardt method.

     The proposed step is now tested by evaluating the function at the
     resulting point, x'.  If the step reduces the norm of the function
     sufficiently, and follows the predicted behavior of the function
     within the trust region, then it is accepted and the size of the
     trust region is increased.  If the proposed step fails to improve
     the solution, or differs significantly from the expected behavior
     within the trust region, then the size of the trust region is
     decreased and another trial step is computed.

     The algorithm also monitors the progress of the solution and
     returns an error if the changes in the solution are smaller than
     the machine precision.  The possible error codes are,

    `GSL_ETOLF'
          the decrease in the function falls below machine precision

    `GSL_ETOLX'
          the change in the position vector falls below machine
          precision

    `GSL_ETOLG'
          the norm of the gradient, relative to the norm of the
          function, falls below machine precision

     These error codes indicate that further iterations will be
     unlikely to change the solution from its current value.


 - Derivative Solver: gsl_multifit_fdfsolver_lmder
     This is an unscaled version of the LMDER algorithm.  The elements
     of the diagonal scaling matrix D are set to 1.  This algorithm may
     be useful in circumstances where the scaled version of LMDER
     converges too slowly, or the function is already scaled
     appropriately.


File: gsl-ref.info,  Node: Minimization Algorithms without Derivatives,  Next: Computing the covariance matrix of best fit parameters,  Prev: Minimization Algorithms using Derivatives,  Up: Nonlinear Least-Squares Fitting

Minimization Algorithms without Derivatives
===========================================

There are no algorithms implemented in this section at the moment.


File: gsl-ref.info,  Node: Computing the covariance matrix of best fit parameters,  Next: Example programs for Nonlinear Least-Squares Fitting,  Prev: Minimization Algorithms without Derivatives,  Up: Nonlinear Least-Squares Fitting

Computing the covariance matrix of best fit parameters
======================================================

 - Function: int gsl_multifit_covar (const gsl_matrix * J, double
          EPSREL, gsl_matrix * COVAR)
     This function uses the Jacobian matrix J to compute the covariance
     matrix of the best-fit parameters, COVAR.  The parameter EPSREL is
     used to remove linear-dependent columns when J is rank deficient.

     The covariance matrix is given by,

          covar = (J^T J)^{-1}

     and is computed by QR decomposition of J with column-pivoting.  Any
     columns of R which satisfy

          |R_{kk}| <= epsrel |R_{11}|

     are considered linearly-dependent and are excluded from the
     covariance matrix (the corresponding rows and columns of the
     covariance matrix are set to zero).


File: gsl-ref.info,  Node: Example programs for Nonlinear Least-Squares Fitting,  Next: References and Further Reading for Nonlinear Least-Squares Fitting,  Prev: Computing the covariance matrix of best fit parameters,  Up: Nonlinear Least-Squares Fitting

Examples
========

The following example program fits a weighted exponential model with
background to experimental data, Y = A \exp(-\lambda t) + b. The first
part of the program sets up the functions `expb_f' and `expb_df' to
calculate the model and its Jacobian.  The appropriate fitting function
is given by,

     f_i = ((A \exp(-\lambda t_i) + b) - y_i)/\sigma_i

where we have chosen t_i = i.  The Jacobian matrix J is the derivative
of these functions with respect to the three parameters (A, \lambda,
b).  It is given by,

     J_{ij} = d f_i / d x_j

where x_0 = A, x_1 = \lambda and x_2 = b.

     #include <stdlib.h>
     #include <stdio.h>
     #include <gsl/gsl_rng.h>
     #include <gsl/gsl_randist.h>
     #include <gsl/gsl_vector.h>
     #include <gsl/gsl_blas.h>
     #include <gsl/gsl_multifit_nlin.h>
     
     struct data {
       size_t n;
       double * y;
       double * sigma;
     };
     
     int
     expb_f (const gsl_vector * x, void *params,
             gsl_vector * f)
     {
       size_t n = ((struct data *)params)->n;
       double *y = ((struct data *)params)->y;
       double *sigma = ((struct data *) params)->sigma;
     
       double A = gsl_vector_get (x, 0);
       double lambda = gsl_vector_get (x, 1);
       double b = gsl_vector_get (x, 2);
     
       size_t i;
     
       for (i = 0; i < n; i++)
         {
           /* Model Yi = A * exp(-lambda * i) + b */
           double t = i;
           double Yi = A * exp (-lambda * t) + b;
           gsl_vector_set (f, i, (Yi - y[i])/sigma[i]);
         }
     
       return GSL_SUCCESS;
     }
     
     int
     expb_df (const gsl_vector * x, void *params,
              gsl_matrix * J)
     {
       size_t n = ((struct data *)params)->n;
       double *sigma = ((struct data *) params)->sigma;
     
       double A = gsl_vector_get (x, 0);
       double lambda = gsl_vector_get (x, 1);
     
       size_t i;
     
       for (i = 0; i < n; i++)
         {
           /* Jacobian matrix J(i,j) = dfi / dxj, */
           /* where fi = (Yi - yi)/sigma[i],      */
           /*       Yi = A * exp(-lambda * i) + b  */
           /* and the xj are the parameters (A,lambda,b) */
           double t = i;
           double s = sigma[i];
           double e = exp(-lambda * t);
           gsl_matrix_set (J, i, 0, e/s);
           gsl_matrix_set (J, i, 1, -t * A * e/s);
           gsl_matrix_set (J, i, 2, 1/s);
     
     
         }
       return GSL_SUCCESS;
     }
     
     int
     expb_fdf (const gsl_vector * x, void *params,
               gsl_vector * f, gsl_matrix * J)
     {
       expb_f (x, params, f);
       expb_df (x, params, J);
     
       return GSL_SUCCESS;
     }

The main part of the program sets up a Levenberg-Marquardt solver and
some simulated random data. The data uses the known parameters
(1.0,5.0,0.1) combined with gaussian noise (standard deviation = 0.1)
over a range of 40 timesteps. The initial guess for the parameters is
chosen as (0.0, 1.0, 0.0).

     #define N 40
     
     int
     main (void)
     {
       const gsl_multifit_fdfsolver_type *T;
       gsl_multifit_fdfsolver *s;
     
       int status;
       size_t i, iter = 0;
     
       const size_t n = N;
       const size_t p = 3;
     
       gsl_matrix *covar = gsl_matrix_alloc (p, p);
     
       double y[N], sigma[N];
     
       struct data d = { n, y, sigma};
     
       gsl_multifit_function_fdf f;
     
       double x_init[3] = { 1.0, 0.0, 0.0 };
     
       gsl_vector_view x = gsl_vector_view_array (x_init, p);
     
       const gsl_rng_type * type;
       gsl_rng * r;
     
       gsl_rng_env_setup();
     
       type = gsl_rng_default;
       r = gsl_rng_alloc (type);
     
       f.f = &expb_f;
       f.df = &expb_df;
       f.fdf = &expb_fdf;
       f.n = n;
       f.p = p;
       f.params = &d;
     
       /* This is the data to be fitted */
     
       for (i = 0; i < n; i++)
         {
           double t = i;
           y[i] = 1.0 + 5 * exp (-0.1 * t)
                      + gsl_ran_gaussian (r, 0.1);
           sigma[i] = 0.1;
           printf ("data: %d %g %g\n", i, y[i], sigma[i]);
         };
     
     
       T = gsl_multifit_fdfsolver_lmsder;
       s = gsl_multifit_fdfsolver_alloc (T, n, p);
       gsl_multifit_fdfsolver_set (s, &f, &x.vector);
     
       print_state (iter, s);
     
       do
         {
           iter++;
           status = gsl_multifit_fdfsolver_iterate (s);
     
           printf ("status = %s\n", gsl_strerror (status));
     
           print_state (iter, s);
     
           if (status)
             break;
     
           status = gsl_multifit_test_delta (s->dx, s->x,
                                             1e-4, 1e-4);
         }
       while (status == GSL_CONTINUE && iter < 500);
     
       gsl_multifit_covar (s->J, 0.0, covar);
     
     #define FIT(i) gsl_vector_get(s->x, i)
     #define ERR(i) sqrt(gsl_matrix_get(covar,i,i))
     
       printf ("A      = %.5f +/- %.5f\n", FIT(0), ERR(0));
       printf ("lambda = %.5f +/- %.5f\n", FIT(1), ERR(1));
       printf ("b      = %.5f +/- %.5f\n", FIT(2), ERR(2));
     
       {
         double chi = gsl_blas_dnrm2(s->f);
         printf("chisq/dof = %g\n",  pow(chi, 2.0)/ (n - p));
       }
     
       printf ("status = %s\n", gsl_strerror (status));
     
       gsl_multifit_fdfsolver_free (s);
       return 0;
     }
     
     void
     print_state (size_t iter, gsl_multifit_fdfsolver * s)
     {
       printf ("iter: %3u x = % 15.8f % 15.8f % 15.8f "
               "|f(x)| = %g\n",
               iter,
               gsl_vector_get (s->x, 0),
               gsl_vector_get (s->x, 1),
               gsl_vector_get (s->x, 2),
               gsl_blas_dnrm2 (s->f));
     }

The iteration terminates when the change in x is smaller than 0.0001, as
both an absolute and relative change.  Here are the results of running
the program:

     iter: 0 x = 1.00000000 0.00000000 0.00000000 |f(x)| = 118.574
     iter: 1 x = 1.64919392 0.01780040 0.64919392 |f(x)| = 77.2068
     iter: 2 x = 2.86269020 0.08032198 1.45913464 |f(x)| = 38.0579
     iter: 3 x = 4.97908864 0.11510525 1.06649948 |f(x)| = 10.1548
     iter: 4 x = 5.03295496 0.09912462 1.00939075 |f(x)| = 6.4982
     iter: 5 x = 5.05811477 0.10055914 0.99819876 |f(x)| = 6.33121
     iter: 6 x = 5.05827645 0.10051697 0.99756444 |f(x)| = 6.33119
     iter: 7 x = 5.05828006 0.10051819 0.99757710 |f(x)| = 6.33119
     
     A      = 5.05828 +/- 0.05983
     lambda = 0.10052 +/- 0.00309
     b      = 0.99758 +/- 0.03944
     chisq/dof = 1.08335
     status = success

The approximate values of the parameters are found correctly, and the
chi-squared value indicates a good fit (the chi-squared per degree of
freedom is approximately 1).  In this case the errors on the parameters
can be estimated from the square roots of the diagonal elements of the
covariance matrix.  If the chi-squared value indicates a poor fit then
error estimates obtained from the covariance matrix are not valid, since
the Gaussian approximation would not apply.


File: gsl-ref.info,  Node: References and Further Reading for Nonlinear Least-Squares Fitting,  Prev: Example programs for Nonlinear Least-Squares Fitting,  Up: Nonlinear Least-Squares Fitting

References and Further Reading
==============================

The MINPACK algorithm is described in the following article,

     J.J. More', `The Levenberg-Marquardt Algorithm: Implementation and
     Theory', Lecture Notes in Mathematics, v630 (1978), ed G. Watson.

The following paper is also relevant to the algorithms described in this
section,

     J.J. More', B.S. Garbow, K.E. Hillstrom, "Testing Unconstrained
     Optimization Software", ACM Transactions on Mathematical Software,
     Vol 7, No 1 (1981), p 17-41.


File: gsl-ref.info,  Node: Physical Constants,  Next: IEEE floating-point arithmetic,  Prev: Nonlinear Least-Squares Fitting,  Up: Top

Physical Constants
******************

This chapter describes macros for the values of physical constants, such
as the speed of light, c, and gravitational constant, G.  The values
are available in different unit systems, including the standard MKSA
system (meters, kilograms, seconds, amperes) and the CGSM system
(centimeters, grams, seconds, gauss), which is commonly used in
Astronomy.

   The definitions of constants in the MKSA system are available in the
file `gsl_const_mksa.h'.  The constants in the CGSM system are defined
in `gsl_const_cgsm.h'.  Dimensionless constants, such as the fine
structure constant, which are pure numbers are defined in
`gsl_const_num.h'.

* Menu:

* Fundamental Constants::
* Astronomy and Astrophysics::
* Atomic and Nuclear Physics::
* Measurement of Time::
* Imperial Units ::
* Nautical Units::
* Printers Units::
* Volume::
* Mass and Weight ::
* Thermal Energy and Power::
* Pressure::
* Viscosity::
* Light and Illumination::
* Radioactivity::
* Force and Energy::
* Prefixes::
* Physical Constant Examples::
* Physical Constant References and Further Reading::

   The full list of constants is described briefly below.  Consult the
header files themselves for the values of the constants used in the
library.


File: gsl-ref.info,  Node: Fundamental Constants,  Next: Astronomy and Astrophysics,  Up: Physical Constants

Fundamental Constants
=====================

`GSL_CONST_MKSA_SPEED_OF_LIGHT'
     The speed of light in vacuum, c.

`GSL_CONST_MKSA_VACUUM_PERMEABILITY'
     The permeability of free space, \mu_0. This constant is defined in
     the MKSA system only.

`GSL_CONST_MKSA_VACUUM_PERMITTIVITY'
     The permittivity of free space, \epsilon_0.  This constant is
     defined in the MKSA system only.

`GSL_CONST_MKSA_PLANCKS_CONSTANT_H'
     Planck's constant, h.

`GSL_CONST_MKSA_PLANCKS_CONSTANT_HBAR'
     Planck's constant divided by 2\pi, \hbar.

`GSL_CONST_NUM_AVOGADRO'
     Avogadro's number, N_a.

`GSL_CONST_MKSA_FARADAY'
     The molar charge of 1 Faraday.

`GSL_CONST_MKSA_BOLTZMANN'
     The Boltzmann constant, k.

`GSL_CONST_MKSA_MOLAR_GAS'
     The molar gas constant, R_0.

`GSL_CONST_MKSA_STANDARD_GAS_VOLUME'
     The standard gas volume, V_0.

`GSL_CONST_MKSA_STEFAN_BOLTZMANN_CONSTANT'
     The Stefan-Boltzmann radiation constant, \sigma.

`GSL_CONST_MKSA_GAUSS'
     The magnetic field of 1 Gauss.

`GSL_CONST_MKSA_MICRON'
     The length of 1 micron.

`GSL_CONST_MKSA_HECTARE'
     The area of 1 hectare.

`GSL_CONST_MKSA_MILES_PER_HOUR'
     The speed of 1 mile per hour.

`GSL_CONST_MKSA_KILOMETERS_PER_HOUR'
     The speed of 1 kilometer per hour.


File: gsl-ref.info,  Node: Astronomy and Astrophysics,  Next: Atomic and Nuclear Physics,  Prev: Fundamental Constants,  Up: Physical Constants

Astronomy and Astrophysics
==========================

`GSL_CONST_MKSA_ASTRONOMICAL_UNIT'
     The length of 1 astronomical unit (mean earth-sun distance), au.

`GSL_CONST_MKSA_GRAVITATIONAL_CONSTANT'
     The gravitational constant, G.

`GSL_CONST_MKSA_LIGHT_YEAR'
     The distance of 1 light-year, ly.

`GSL_CONST_MKSA_PARSEC'
     The distance of 1 parsec, pc.

`GSL_CONST_MKSA_GRAV_ACCEL'
     The standard gravitational acceleration on Earth, g.

`GSL_CONST_MKSA_SOLAR_MASS'
     The mass of the Sun.


File: gsl-ref.info,  Node: Atomic and Nuclear Physics,  Next: Measurement of Time,  Prev: Astronomy and Astrophysics,  Up: Physical Constants

Atomic and Nuclear Physics
==========================

`GSL_CONST_MKSA_ELECTRON_CHARGE'
     The charge of the electron, e.

`GSL_CONST_MKSA_ELECTRON_VOLT'
     The energy of 1 electron volt, eV.

`GSL_CONST_MKSA_UNIFIED_ATOMIC_MASS'
     The unified atomic mass, amu.

`GSL_CONST_MKSA_MASS_ELECTRON'
     The mass of the electron, m_e.

`GSL_CONST_MKSA_MASS_MUON'
     The mass of the muon, m_\mu.

`GSL_CONST_MKSA_MASS_PROTON'
     The mass of the proton, m_p.

`GSL_CONST_MKSA_MASS_NEUTRON'
     The mass of the neutron, m_n.

`GSL_CONST_NUM_FINE_STRUCTURE'
     The electromagnetic fine structure constant \alpha.

`GSL_CONST_MKSA_RYDBERG'
     The Rydberg constant, Ry, in units of energy.  This is related to
     the Rydberg inverse wavelength R by Ry = h c R.

`GSL_CONST_MKSA_BOHR_RADIUS'
     The Bohr radius, a_0.

`GSL_CONST_MKSA_ANGSTROM'
     The length of 1 angstrom.

`GSL_CONST_MKSA_BARN'
     The area of 1 barn.

`GSL_CONST_MKSA_BOHR_MAGNETON'
     The Bohr Magneton, \mu_B.

`GSL_CONST_MKSA_NUCLEAR_MAGNETON'
     The Nuclear Magneton, \mu_N.

`GSL_CONST_MKSA_ELECTRON_MAGNETIC_MOMENT'
     The absolute value of the magnetic moment of the electron, \mu_e.
     The physical magnetic moment of the electron is negative.

`GSL_CONST_MKSA_PROTON_MAGNETIC_MOMENT'
     The magnetic moment of the proton, \mu_p.

`GSL_CONST_MKSA_THOMSON_CROSS_SECTION'
     The Thomson cross section, \sigma_T.


File: gsl-ref.info,  Node: Measurement of Time,  Next: Imperial Units,  Prev: Atomic and Nuclear Physics,  Up: Physical Constants

Measurement of Time
===================

`GSL_CONST_MKSA_MINUTE'
     The number of seconds in 1 minute.

`GSL_CONST_MKSA_HOUR'
     The number of seconds in 1 hour.

`GSL_CONST_MKSA_DAY'
     The number of seconds in 1 day.

`GSL_CONST_MKSA_WEEK'
     The number of seconds in 1 week.


File: gsl-ref.info,  Node: Imperial Units,  Next: Nautical Units,  Prev: Measurement of Time,  Up: Physical Constants

Imperial Units
==============

`GSL_CONST_MKSA_INCH'
     The length of 1 inch.

`GSL_CONST_MKSA_FOOT'
     The length of 1 foot.

`GSL_CONST_MKSA_YARD'
     The length of 1 yard.

`GSL_CONST_MKSA_MILE'
     The length of 1 mile.

`GSL_CONST_MKSA_MIL'
     The length of 1 mil (1/1000th of an inch).


File: gsl-ref.info,  Node: Nautical Units,  Next: Printers Units,  Prev: Imperial Units,  Up: Physical Constants

Nautical Units
==============

`GSL_CONST_MKSA_NAUTICAL_MILE'
     The length of 1 nautical mile.

`GSL_CONST_MKSA_FATHOM'
     The length of 1 fathom.

`GSL_CONST_MKSA_KNOT'
     The speed of 1 knot.


File: gsl-ref.info,  Node: Printers Units,  Next: Volume,  Prev: Nautical Units,  Up: Physical Constants

Printers Units
==============

`GSL_CONST_MKSA_POINT'
     The length of 1 printer's point (1/72 inch).

`GSL_CONST_MKSA_TEXPOINT'
     The length of 1 TeX point (1/72.27 inch).


File: gsl-ref.info,  Node: Volume,  Next: Mass and Weight,  Prev: Printers Units,  Up: Physical Constants

Volume
======

`GSL_CONST_MKSA_ACRE'
     The area of 1 acre.

`GSL_CONST_MKSA_LITER'
     The volume of 1 liter.

`GSL_CONST_MKSA_US_GALLON'
     The volume of 1 US gallon.

`GSL_CONST_MKSA_CANADIAN_GALLON'
     The volume of 1 Canadian gallon.

`GSL_CONST_MKSA_UK_GALLON'
     The volume of 1 UK gallon.

`GSL_CONST_MKSA_QUART'
     The volume of 1 quart.

`GSL_CONST_MKSA_PINT'
     The volume of 1 pint.


File: gsl-ref.info,  Node: Mass and Weight,  Next: Thermal Energy and Power,  Prev: Volume,  Up: Physical Constants

Mass and Weight
===============

`GSL_CONST_MKSA_POUND_MASS'
     The mass of 1 pound.

`GSL_CONST_MKSA_OUNCE_MASS'
     The mass of 1 ounce.

`GSL_CONST_MKSA_TON'
     The mass of 1 ton.

`GSL_CONST_MKSA_METRIC_TON'
     The mass of 1 metric ton (1000 kg).

`GSL_CONST_MKSA_UK_TON'
     The mass of 1 UK ton.

`GSL_CONST_MKSA_TROY_OUNCE'
     The mass of 1 troy ounce.

`GSL_CONST_MKSA_CARAT'
     The mass of 1 carat.

`GSL_CONST_MKSA_GRAM_FORCE'
     The force of 1 gram weight.

`GSL_CONST_MKSA_POUND_FORCE'
     The force of 1 pound weight.

`GSL_CONST_MKSA_KILOPOUND_FORCE'
     The force of 1 kilopound weight.

`GSL_CONST_MKSA_POUNDAL'
     The force of 1 poundal.


File: gsl-ref.info,  Node: Thermal Energy and Power,  Next: Pressure,  Prev: Mass and Weight,  Up: Physical Constants

Thermal Energy and Power
========================

`GSL_CONST_MKSA_CALORIE'
     The energy of 1 calorie.

`GSL_CONST_MKSA_BTU'
     The energy of 1 British Thermal Unit, btu.

`GSL_CONST_MKSA_THERM'
     The energy of 1 Therm.

`GSL_CONST_MKSA_HORSEPOWER'
     The power of 1 horsepower.


File: gsl-ref.info,  Node: Pressure,  Next: Viscosity,  Prev: Thermal Energy and Power,  Up: Physical Constants

Pressure
========

`GSL_CONST_MKSA_BAR'
     The pressure of 1 bar.

`GSL_CONST_MKSA_STD_ATMOSPHERE'
     The pressure of 1 standard atmosphere.

`GSL_CONST_MKSA_TORR'
     The pressure of 1 torr.

`GSL_CONST_MKSA_METER_OF_MERCURY'
     The pressure of 1 meter of mercury.

`GSL_CONST_MKSA_INCH_OF_MERCURY'
     The pressure of 1 inch of mercury.

`GSL_CONST_MKSA_INCH_OF_WATER'
     The pressure of 1 inch of water.

`GSL_CONST_MKSA_PSI'
     The pressure of 1 pound per square inch.


File: gsl-ref.info,  Node: Viscosity,  Next: Light and Illumination,  Prev: Pressure,  Up: Physical Constants

Viscosity
=========

`GSL_CONST_MKSA_POISE'
     The dynamic viscosity of 1 poise.

`GSL_CONST_MKSA_STOKES'
     The kinematic viscosity of 1 stokes.


File: gsl-ref.info,  Node: Light and Illumination,  Next: Radioactivity,  Prev: Viscosity,  Up: Physical Constants

Light and Illumination
======================

`GSL_CONST_MKSA_STILB'
     The luminance of 1 stilb.

`GSL_CONST_MKSA_LUMEN'
     The luminous flux of 1 lumen.

`GSL_CONST_MKSA_LUX'
     The illuminance of 1 lux.

`GSL_CONST_MKSA_PHOT'
     The illuminance of 1 phot.

`GSL_CONST_MKSA_FOOTCANDLE'
     The illuminance of 1 footcandle.

`GSL_CONST_MKSA_LAMBERT'
     The luminance of 1 lambert.

`GSL_CONST_MKSA_FOOTLAMBERT'
     The luminance of 1 footlambert.


File: gsl-ref.info,  Node: Radioactivity,  Next: Force and Energy,  Prev: Light and Illumination,  Up: Physical Constants

Radioactivity
=============

`GSL_CONST_MKSA_CURIE'
     The activity of 1 curie.

`GSL_CONST_MKSA_ROENTGEN'
     The exposure of 1 roentgen.

`GSL_CONST_MKSA_RAD'
     The absorbed dose of 1 rad.

