Providing the function to solve

You must provide @math{n} functions of @math{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 @math{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,

with @math{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;
   double A = (params->A);
   double x0 = gsl_vector_get(x,0);
   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.function = &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 @math{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 @math{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 @math{f(x)} and @math{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;
   double A = (params->A);
   double x0 = gsl_vector_get(x,0);
   double x1 = gsl_vector_get(x,1);
   gsl_vector_set (J, 0, 0, A * x1)
   gsl_vector_set (J, 0, 1, A * x0)
   gsl_vector_set (J, 1, 0, -exp(-x0))
   gsl_vector_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;
   double A = (params->A);
   double x0 = gsl_vector_get(x,0);
   double x1 = gsl_vector_get(x,1);

   double u0 = exp(-x0);
   double u1 = exp(-x1);

   gsl_vector_set (f, 0, A * x0 * x1 - 1)
   gsl_vector_set (f, 1, u0 + u1 - (1 + 1/A))

   gsl_vector_set (J, 0, 0, A * x1)
   gsl_vector_set (J, 0, 1, A * x0)
   gsl_vector_set (J, 1, 0, -u0)
   gsl_vector_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.