gsl multiroot iteration trying nan?

Question

I am trying to find the right parameters to input to my code to produced the desired results. Instead of guessing and checking I am using a root find to find the parameters that give the desired results. There are two variables that are free to vary, but I was having difficulty running the root finder. I changed the code to solve for each variable individually and found out that I was having trouble optimizing one variable.

It seems that the problem is that gsl_multiroot_iterate is guessing nan for x1 after the first iteration. At least that is what the value of x1 is returning in the function() call after that point, when I put in a printf statement for x1.

The simulation I am running only allows values of x1 between 0 and 1. It could be possible that this is causing the issue, though I check in the simulation to make sure x1 is between 0 and 1, and never throws an issue besides when x1 is nan. Is there anyway to set a range for what values the iteration tries for x1? And would anyone know what the iteration tries using nan for x1?

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_multiroots.h>

struct rparams{
    double target1;
};

int function(const gsl_vector * x, void *params, gsl_vector * f);

int main(int argc, char* argv[]) {

    double target1;
    sscanf(argv[1],"%lf",&target1);

    const gsl_multiroot_fsolver_type *T;
    gsl_multiroot_fsolver *s;

    int status;
    unsigned int iter = 0;

    const size_t n = 1;
    struct rparams p;
    p.target1 = target1;

    gsl_multiroot_function f = {&function, n, &p};

    double x_init[1] = {.1};
    gsl_vector * x = gsl_vector_alloc(n);

    gsl_vector_set(x, 0, x_init[0]);

    T = gsl_multiroot_fsolver_hybrid;
    s = gsl_multiroot_fsolver_alloc(T, 1);
    gsl_multiroot_fsolver_set(s, &f, x);

    print_state(iter, s);

    do
    {
        iter++;
        status = gsl_multiroot_fsolver_iterate (s);

        print_state(iter, s);

        /* check if solver is stuck */
        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_fsolver_free (s);
    gsl_vector_free (x);
    return 0;
}

int function(const gsl_vector * x, void *params, gsl_vector * f){

    double target1 = ((struct rparams *) params)->target1;

    double x1 = gsl_vector_get(x, 0);

    /* Run simulation here using x1 parameter */
    /* Assign output to temp1, which I am trying to match to target1 */

    const double y1 = temp1 - target1;

    gsl_vector_set (f, 0, y1);

    return GSL_SUCCESS;
}

Answer 1

Be careful in designing the function you want to obtain the root from. In fact, for a test, I tried a function that had a constant output. This caused the algorithm to throw out the NaNs.

Answer 2

If you only need to find the root of a single equation, you can use the gsl_roots library instead of gsl_multiroots. The gsl_roots library has several bisection algorithms for which you specify a range instead of an initial guess. If you know your root is inside the interval (0, 1), you would set that as the target interval and the algorithm would never go outside that range. A minimal, complete example in C++ demonstrating the bisection method is below. If you can't use C++11 lambda functions, then you'd have to define the objective function like you did in your original question.

#include <iostream>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_roots.h>

using namespace std;

int
main (void)
{
    //Set the solver type (bisection method)
    gsl_root_fsolver* s = gsl_root_fsolver_alloc(gsl_root_fsolver_bisection);

    //Use a lambda to define the objective function.
    //This is a parabola with the equation: y = (x-1)^2 - 1
    //It has roots at x = 0 and x = 2.
    gsl_function F;
    F.function = [](double x, void*){return ((x-1) * (x-1)) - 1;};

    //Initialize the solver; make a guess that the root is between x = 0.5 and x = 10
    gsl_root_fsolver_set(s, &F, 0.5, 10.0);

    //Run the solver until the root is found to within 0.001
    int status;
    do {
        gsl_root_fsolver_iterate(s);
        double r = gsl_root_fsolver_root(s);
        double x_low = gsl_root_fsolver_x_lower(s);
        double x_high = gsl_root_fsolver_x_upper(s);
        status = gsl_root_test_interval(x_low, x_high, 0, 0.001);

        if (status == GSL_SUCCESS)
            cout << "Converged" << endl;

        cout << "x_low = " << x_low;
        cout << "; x_high = " << x_high;
        cout << "; root = " << r << endl;
    }
    while (status == GSL_CONTINUE);

    return status;
}