GSL 中的错误 - 根查找

Question

大家好，我在使我的程序与 GSL 一起工作时遇到问题 - 根查找。我正在尝试找到我的方程式的解决方案。我正在寻找 64 行数据的解决方案，但在某些特定行中程序无法继续，可能是因为不存在好的解决方案。但我希望程序在他找不到解决方案时跳过行。但是我的程序有时会停止并出现以下消息： gsl: brent.c:74: 错误：端点不跨越 y=0 调用了 passou1passou2passou3Default GSL 错误处理程序。 中止陷阱：6

所以，我做了一些打印来检查我的程序到底在哪里停止，我发现它在 gsl_root_fsolver_set(s,&F, x_lo, x_hi) 中，但我没有找到如何打印这个值或者这个功能给了我什么。

我的程序到了，谢谢大家！

#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

加布里埃拉。输出已经告诉你为什么你的程序是错误的。 端点不跨越 y=0。

本节中描述的根包围算法需要一个保证包含根的初始区间——如果 a 和 b 是区间的端点，则 f(a) 的符号必须与 f(b) 不同。

以上来自manual of gsl，所以如果端点有相同的符号，程序会停止并告诉你这个错误。

你有没有尝试过 gsl 中的错误处理程序。在手册的第 3 章中，他们提供了一个名为 gsl_set_error_handler_off() 的函数，如果您将此函数放在 gsl_root_fsolver_set(s,&F, x_lo, x_hi) 之前，您可以将此函数分配给 int 类型变量，让我们比如status，从手册中我们可以看出它（gsl_root_fsolver_set(s,&F, x_lo, x_hi)）是一个int function，那么如果你打印出status的值，并用gsl_errno.h文件检查它，你就会知道是什么那个值的意思。

这个 gsl_set_error_handler_off() 可以取消执行中止。

至于您的代码，您应该执行以下操作：

1. 添加#include <gsl/gsl_errno.h>
2. 在status=gsl_root_fsolver_set(s,&F, x_lo, x_hi)之前添加gsl_set_error_handler_off()
3. 使用status的值，就像你程序的下面的循环一样，对你的端点做一些小的改变，扩大或缩小或翻译你的区间，当它们满足初始条件时，程序将再次运行