Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Help using gsl_multiroot_fsolver in C

  1. Nov 14, 2012 #1
    Dear all,

    Unfortunately I am not skilled in GSL. I am using gsl_multiroot_fsolver_hybrid in C 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 your time!
    Yuu
     
  2. jcsd
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Can you offer guidance or do you also need help?
Draft saved Draft deleted



Similar Discussions: Help using gsl_multiroot_fsolver in C
  1. Help with C (Replies: 9)

  2. Help with C++ (Replies: 5)

Loading...