36.3 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.