C++程序的奇怪行为

Question

这段话的翻译为：

在使用二分法计算一个函数根的过程中，我发现了一个奇怪的行为。该问题需要求解一个涉及积分的系统，其中积分的上限是另一个函数的函数，这本身并不是什么大问题。当我直接将这个数值写入程序时，程序运行正常结束；但如果我尝试回调这个函数，则终端会终止进程，并且程序运行变慢。顺便提一下，这是我用C++编写的首批代码之一。

相关的代码块如下：

#include 
#include 
#include 
#include 

// 结构体定义，用于存储参数
struct Ninj_Params{
    double p;
    double Q0;
    double gamma_m;
};

// 定义物理常量
double c = 1; // 光速
double eV = 1.0; // 电子伏特
double MeV = pow(10, 6) * eV; // 毫电子伏特
double me = 0.510998902 * MeV; // 电子质量
// ... 其他相关物理常量定义

// 二分法求根函数
double bisection(const std::function &func, double a, double b, double tolerance){
    if (func(a) * func(b) >= 0){
        printf("二分法在该区间无法保证收敛。
");
        return NAN; // 返回非数字表示失败
    }

    double c = a;
    while ((b - a) >= tolerance) {
        // 找到中间点
        c = (a + b) / 2;

        // 检查是否找到根
        if (func(c) == 0.0)
            break;

        // 根据中间点处函数值的符号更新区间
        if (func(c) * func(a) < 0)
            b = c;
        else
            a = c;
    }

    return c;
}

// 计算β函数
double beta(double Gamma){
    return sqrt(1. - pow(Gamma, -2.));
}

// 计算径向距离R
double R(double t, double Rs, double Gamma){
    return (Rs + beta(Gamma)*c*Gamma*t) * 1./cm;
}

// 计算磁场强度B
double B(double t, double Rs, double R0, double Gamma, double B0, double a){
    return B0 * pow((R(t, Rs, Gamma)*cm) / R0, -a) * 1./Gauss;
}

// 定义发射率Q函数
double Q(double gammae, void * args){
    Ninj_Params ¶ms = *(Ninj_Params*)args;
    
    if (params.gamma_m <= gammae){
        return params.Q0 * pow(gammae / params.gamma_m, -params.p) * sec;
    } else {
        return 0.0;
    }
}

// 计算总注入粒子数Ninj
double Ninj(double Q0, double p, double gamma_min, double gamma_max){
    struct Ninj_Params params = {
        .p = 2.5,
        .Q0 = Q0,
        .gamma_m = gamma_min,
    };

    gsl_integration_workspace *w = gsl_integration_workspace_alloc (1000);
    double result, error;
    gsl_function F;

    F.function = &Q;
    F.params = ¶ms;

    gsl_integration_qags(&F, gamma_min, gamma_max, 0, 1e-7, 1000, w, &result, &error);

    return result;
}

// 计算Q0随时间的变化率Q0rate
double Q0rate(double Ninj_Gest, double t, double Rs, double R0, double Gamma, double B0, double a){
    //double gamma_max = 1e8*pow(B(t, Rs, R0, Gamma, B0, a), -0.5);
    double gamma_max = 7.957e+06;
    printf("**** %.3e
", gamma_max);
    double p = 2.5;

    auto root2solve = [Ninj_Gest, p, &gamma_max](double Q0) -> double{
        return Ninj_Gest - Ninj(Q0*1/sec, p, 1e5, gamma_max);
    };

    double root = bisection(root2solve, 1e40*1./sec, 1e53*1./sec, 1e-7);
    return root*sec;
}

// 测试用例函数
void TEST_ROOT(double t, double Rs, double R0, double Gamma, double gamma_m, double B0, double a){
    double Nj = 1e47*1/sec;
    printf("%.5e
", Q0rate(Nj, t, Rs, R0, Gamma, B0, a));
}

// 主函数入口
int main(){
    double t = 10*sec;
    double Rs = 1e14*cm;
    double R0 = 1e15*cm;
    double B0 = 0.3*1e2*Gauss;
    double Gamma = 300;
    double gamma_m = 1e5;
    double a = 1.0;
    TEST_ROOT(t, Rs, R0, Gamma, gamma_m, B0, a);
}

问题行出现在第105行（以及注释掉的106行）。当设置为：

double gamma_max = 1e8*pow(B(t, Rs, R0, Gamma, B0, a), -0.5);

程序输出表现为（运行较慢）：

**** 7.957e+06
zsh: killed     ./Main

而如果将其改为：

double gamma_max = 7.957e+06;

程序会正常结束，并显示预期的根值：

% ./Main
**** 7.957e+06
1.50212e+42

miku · Answer

你不能期望通过double类型的浮点数运算得到精确的结果，即使是像 3x=1 这样简单的方程也无法在浮点数数学中精确解决（因为1/3这个数不能被double类型准确表示）。

你的退出条件应该基于某个ε（epsilon，即容许误差）值。

按照IEEE 754双精度浮点数标准（通常C++中的double类型所采用的精度），只有53位尾数部分，因此需要根据涉及数值的大小来计算ε值。例如，如果你的函数可能输出“很大”（如1E50）的值，那么单位级别的ε值就不足够了（对于这样的大数，1E50+1实际上等于1E50）。

另外要注意的是，浮点数不精确性甚至可能在中间计算过程中产生问题，特别是在进行大的数值抵消操作时（比如当你减去两个非常接近的大数时）。

正如Ted Lyngmo指出的那样，原始代码中还存在内存泄漏问题。

Mark Setchell · Answer

您分配了大量的内存，这很可能导致寻找新的内存时变得越来越慢。如果在使用完毕后释放这部分内存，程序就能正常运行，并且在我这里得到了如下输出：

**** 7.957e+06
1.50212e+42

我唯一对您的程序做出的修改是在以下位置添加了一行代码：

double Ninj(double Q0, double p, double gamma_min, double gamma_max) {
    struct Ninj_Params params = {
        .p = 2.5,
        .Q0 = Q0,
        .gamma_m = gamma_min,
    };

    gsl_integration_workspace* w = gsl_integration_workspace_alloc(1000);
    double result, error;
    gsl_function F;

    F.function = &Q;
    F.params = ¶ms;

    gsl_integration_qags(&F, gamma_min, gamma_max, 0, 1e-7, 1000, w, &result, &error);

    gsl_integration_workspace_free(w);  // <- 在这里添加了释放内存的语句

    return result;
}