error in GSL - root finding

Question

Hello everybody I'm with problems to make my program work with GSL - root finding. I'm trying find a solution to my equation. I'm looking solutions for a data with 64 lines, but in some specific lines the program can't continue, maybe because a good solution is not existent. But I want that the program just skip line when he don't find a solution. But my program just stop sometimes and appear this message: gsl: brent.c:74: ERROR: endpoints do not straddle y=0 passou1passou2passou3Default GSL error handler invoked. Abort trap: 6

So, I did some prints to check where exactly my program stop and I found that is in gsl_root_fsolver_set(s,&F, x_lo, x_hi), but I didn't find how print this value or what this function give to me.

My program is here, thanks everybody!

#include <stdio.h>
#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_roots.h>
#include "demo_fn.h"
#include "demo_fn.c"

int
main (void)
{

double NT, c;
c = 64.0/2.0;

char url[]="charm.txt";
double corr, core, a[64][3];
int nt, i = 0;
FILE *arq;
arq = fopen(url, "r");
if(arq == NULL)
    printf("Erro, nao foi possivel abrir o arquivo\n");
else
    while( (fscanf(arq,"%lf %lf %i\n", &corr, &core, &nt))!=EOF ) {
        a[i][0]=corr;
        a[i][1]=core;
        a[i][2]=nt;
        i++;
        //printf("%lf, %lf, %i\n",corr, core, nt);
    }

fclose(arq);




for (i= 0; i < 64; i++)
{
    int status;
    int iter = 0, max_iter = 200;
    const gsl_root_fsolver_type *T;
    gsl_root_fsolver *s;
    double r = 0, r_expected = 4.0;
    double x_lo = 0.0001, x_hi = 4.0;
    double ratio1, ratio2;

    ratio1 = a[i][0]/a[i+1][0];
    ratio2 = a[i+1][0]/a[i+2][0];
    printf ("ratio1: %lf, ratio2: %lf", ratio1, ratio2);
    printf ("\n");

 //   if (ratio1*ratio2 > 0)
 //      {

    printf("C(n_t) : %.15lf -- loop index : %i ----- ratio: %lf \n", a[i][0],i, ratio1);

    gsl_function F;
    struct quadratic_params params = {a[i][0], i, c, i+1, a[i+1][0]};
    F.function = &quadratic;
           printf ("passou1");
    F.params = &params;
    T = gsl_root_fsolver_brent;
           printf ("passou2");
    //T = gsl_root_fsolver_bisection;
    s = gsl_root_fsolver_alloc (T);
        printf ("passou3");
    gsl_root_fsolver_set (s, &F, x_lo, x_hi);
        printf ("passou4");
    printf ("using %s method\n", gsl_root_fsolver_name (s));
    printf ("%5s [%9s, %9s] %9s %10s %9s\n", "iter", "lower", "upper", "root", "err", "err(est)");


    do
    {
        iter++;
        status = gsl_root_fsolver_iterate (s);
        r = gsl_root_fsolver_root (s);
        x_lo = gsl_root_fsolver_x_lower (s);
        x_hi = gsl_root_fsolver_x_upper (s);
        status = gsl_root_test_interval (x_lo, x_hi,0, 0.001);
        if (status == GSL_SUCCESS)
        {
            printf ("Converged:\n");
        }
        printf ("%5d [%.7lf, %.7lf] %.7lf %+.7lf %.7lf\n", iter, x_lo, x_hi, r, r - r_expected, x_hi - x_lo);


    }
    while (status == GSL_CONTINUE && iter < max_iter);

    gsl_root_fsolver_free (s);

 //     }

    printf("\n");
}



return 0;

}

Answer 1

Gabriela. The output has already told you why you program is wrong. Endpoints do not straddle y=0.

The root bracketing algorithms described in this section require an initial interval which is guaranteed to contain a root—if a and b are the endpoints of the interval then f(a) must differ in sign from f(b).

The above is from the manual of gsl, so if the end points have the same sign, the program will stop and tell you this error.

Have you ever tried the error handler in gsl. In the chapter 3 of the manual, they offer a function called gsl_set_error_handler_off(), if you place this function just before gsl_root_fsolver_set(s,&F, x_lo, x_hi), and you can assign this function to an int type variable, let's say status, as we can see from the manual that it(gsl_root_fsolver_set(s,&F, x_lo, x_hi)) is an int function, then if you print out the value of status, and check it with gsl_errno.h file, you will know what that value mean.

This gsl_set_error_handler_off() can dismiss the abortion from execution.

as for you code, you should do the following:

1. add #include <gsl/gsl_errno.h>
2. add gsl_set_error_handler_off() before status=gsl_root_fsolver_set(s,&F, x_lo, x_hi)
3. use the value of status like the following loop of you program to make little changes of you end points,to expand or narrow or translate your interval, when they satisfied the initial condition, the program will run again