Блог

Я хочу использовать процедуру деления GSL на две части для поиска корней функции f(x)=1/(2sin(x)^2) log(sin(x))-k, константа k, на интервале [0;Pi/2].

      int main(void *params){
  
     struct func_params *part= (struct func_params*)params;
    
     int status;
     int iter = 0, max_iter = 10;
     const gsl_root_fsolver_type *R;
     gsl_root_fsolver *s;
     double x_lo = 1e-4, x_hi = M_PI/2.;

     gsl_function F;
     F.function = amp;my_func;
     F.params = amp;params;
     R = gsl_root_fsolver_bisection;
     s = gsl_root_fsolver_alloc(R);
     gsl_root_fsolver_set(s,amp;F,x_lo,x_hi);

     printf ("using %s methodn",gsl_root_fsolver_name(s));
     printf ("%5s [%9s, %9s] %9s s %9sn","iter", "lower", "upper", "root","err", "err(est)");

    return(0);
 

}

Этот код возвращает gsl: bisection.c:55: ERROR: function value is not finite Default GSL error handler invoked. Проблема в том, что, когда я смотрю на файл bisection.c, строка 55 не соответствует предыдущей ошибке…Есть идеи, что происходит?

Вы должны указать интервал таким образом, чтобы функция пересекалась y = 0 только один раз.

Например, [0.1, 2.0] здесь:

 #include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_roots.h>

double my_func(double x,  void* unused)
{
    double k = 1;
    return 1 / (2 * pow(sin(x), 2))   log(sin(x)) - k;
}

int main()
{
    gsl_function F;
    F.function = amp;my_func;
    F.params = 0;

    gsl_root_fsolver *s = gsl_root_fsolver_alloc(gsl_root_fsolver_bisection);

    double x_lo =  0.1, x_hi =  2.0;
    gsl_root_fsolver_set(s, amp;F, x_lo, x_hi);

    for(int i = 0; i < 1000; i  ) {

        int status = gsl_root_fsolver_iterate(s);

        double left_int  = gsl_root_fsolver_x_lower(s);
        double right_int = gsl_root_fsolver_x_upper(s);

        printf("iteration d: [%.04lf, %.04lf]n", i, left_int, right_int);

        status = gsl_root_test_interval(left_int, right_int, 1.0e-5, 1.0e-15);
        if(status != GSL_CONTINUE) {

            printf("status: %sn", gsl_strerror(status));
            printf("nRoot interval = [%.010lf, %.010lf]n", left_int, right_int);
            break;
        }
    }
}
 

Выход:

 $ gcc main.c -lgsl -lm amp;amp; ./a.out
iteration 000: [0.1000, 1.0500]
iteration 001: [0.5750, 1.0500]
iteration 002: [0.5750, 0.8125]
iteration 003: [0.5750, 0.6938]
iteration 004: [0.5750, 0.6344]
iteration 005: [0.5750, 0.6047]
iteration 006: [0.5898, 0.6047]
iteration 007: [0.5973, 0.6047]
iteration 008: [0.5973, 0.6010]
iteration 009: [0.5973, 0.5991]
iteration 010: [0.5982, 0.5991]
iteration 011: [0.5987, 0.5991]
iteration 012: [0.5989, 0.5991]
iteration 013: [0.5989, 0.5990]
iteration 014: [0.5989, 0.5990]
iteration 015: [0.5989, 0.5990]
iteration 016: [0.5989, 0.5990]
iteration 017: [0.5989, 0.5990]
status: success

Root interval = [0.5989471436, 0.5989543915]