[Help-gsl] Help using gsl_multiroot_fsolver_hybrid

[Yuu Suzumi]

Dear all, 


I am not skilled in GSL. I am using gsl_multiroot_fsolver_hybrid to solve a 
system of 3 equations with three unknowns. I have tried very hard to find my 
error and repeatedly consulted the GSL manual, but without success. I hope 
there is someone who can offer some insight. I will first post the part of the 
code where the error occurs. Please tell me if it is necessary to send a 
simpler but standalone code displaying the problem. Sofar I post only the two 
functions to keep it short.


The problem occurs in the second (iter =1)  iteration, inside int quasinormal_f 
,  where the vector x (independent variable of the problem, which the solver 
should vary) takes values (nan, nan, nan). I cannot see why.

The functions I defined are:
/* ---  --- */
int quasinormal_f (const gsl_vector * x, void *params,
                   gsl_vector * f)
     {

       const double x0 = gsl_vector_get (x, 0);
       const double x1 = gsl_vector_get (x, 1);
       const double x2 = gsl_vector_get (x, 2);


       qparams->N =x0;
       qparams->A =x1;
       qparams->S =x2;

                     printf("INNER POINT 1 -- Attention!\n"); 
                      printf("N,A,S = %e %e %e \n",
                        x0, x1, x2);
       const double y0 = norm_quasinormal(qparams) - 1.0;
       //const double y0 = 0;

       const double y1 = mean_quasinormal(qparams) - 0.0;
       //const double y1 =0;

       const double y2 = var_quasinormal(qparams) - kappa_var(gcosmo-> theta);


       gsl_vector_set (f, 0, y0);
       gsl_vector_set (f, 1, y1);
       gsl_vector_set (f, 2, y2);

       return GSL_SUCCESS;
     }

int COMPUTE_QUASINORMAL_PARAMETERS()
     {

     const gsl_multiroot_fsolver_type *T;
     gsl_multiroot_fsolver *s;

       int status;
       size_t i, iter = 0;

       const size_t ndim = 3;

       gsl_multiroot_function F;

    F.f = &quasinormal_f;
    F.params = NULL;
    F.n = 3;

       double x_init[3] = {1.5, 0.1,0.5*gcosmo->ausgelagertes_lognormal_sigma};
        gsl_vector *x = gsl_vector_calloc (ndim);

     printf("Check point 2: x = %e %e %e \n",
     gsl_vector_get(x,0), gsl_vector_get(x,1), gsl_vector_get(x,2));


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

     printf("Check point 3: x = %e %e %e \n",
     gsl_vector_get(x,0), gsl_vector_get(x,1), gsl_vector_get(x,2));


       T = gsl_multiroot_fsolver_hybrid;
       s = gsl_multiroot_fsolver_alloc (T, 3);

       gsl_multiroot_fsolver_set (s, &F, x);
       printf("THE SOLVER WAS SET!\n");
       do
         {
           iter++;
                 printf("ABOUT TO ITERATE SOLVER\n");
                             printf("iter = %i\n",
                        iter);


                       printf("x_init = %e %e %e \n",
                        x_init[0], x_init[1], x_init[2]);

                       status = gsl_multiroot_fsolver_iterate (s);

                       printf("THE SOLVER WAS ITERATED\n");
                             printf("iter = %i\n",
                        iter);

                       printf("x_init = %e %e %e \n",
                        x_init[0], x_init[1], x_init[2]);


           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_vector_free (x);

              return 0;
     }     
/* ---  --- */

The output of the programme is: 


Check point 2: x = 0.000000e+00 0.000000e+00 0.000000e+00 
Check point 3: x = 1.500000e+00 1.000000e-01 1.035659e-01 
INNER POINT 1 -- Attention!
N,A,S = 1.500000e+00 1.000000e-01 1.035659e-01 
INNER POINT 1 -- Attention!
N,A,S = 1.500000e+00 1.000000e-01 1.035659e-01 
INNER POINT 1 -- Attention!
N,A,S = 1.500000e+00 1.000000e-01 1.035659e-01 
INNER POINT 1 -- Attention!
N,A,S = 1.500000e+00 1.000000e-01 1.035659e-01 
THE SOLVER WAS SET!
ABOUT TO ITERATE SOLVER
iter = 1
x_init = 1.500000e+00 1.000000e-01 1.035659e-01 
INNER POINT 1 -- Attention!
N,A,S = nan nan nan 

Thank you for any time you invest in this!
YSZ