回调函数有一个奇怪的行为,我想用二分法计算一个函数的根,要解决的系统涉及到一个函数的积分,其中上限是另一个函数的函数,这不是什么大不了的事情,当我把数值放进去的时候,程序结束得很好,但是如果我做函数的回调,终端杀死了进程,程序运行得更慢。顺便说一下,这些是我用c++写的第一段代码。
代码块是下一个;
#include <gsl/gsl_integration.h>
#include <functional>
#include <iostream>
#include <cmath>
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 e = 8.542 * pow(10,-2);
double fm = 1/197.3269602 * 1/MeV;
double cm = pow(10,13) * fm;
double Gauss = 5.788381749 * pow(10,-15) * MeV * 2 * me/e;
double meter = pow(10,2) * cm;
double sec = 299792458*meter;
// Function to find the root using the bisection method
double bisection(const std::function<double(double)> &func, double a, double b, double tolerance){
if (func(a) * func(b) >= 0){
printf("Bisection method cannot guarantee convergence in this interval.\n");
return NAN; // Not a Number to indicate failure
}
double c = a;
while ((b - a) >= tolerance) {
// Find middle point
c = (a + b) / 2;
// Check if the root is found
if (func(c) == 0.0)
break;
// Update the interval based on the sign of the function at the middle point
if (func(c) * func(a) < 0)
b = c;
else
a = c;
}
return c;
}
double beta(double Gamma){
return sqrt(1. - pow(Gamma, -2.));
}
double R(double t, double Rs, double Gamma){
return (Rs + beta(Gamma)*c*Gamma*t)*1./cm;
}
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;
}
double Q(double gammae, void * args){
Ninj_Params& params = *(Ninj_Params*)args;
if(params.gamma_m <= gammae){
return params.Q0*pow(gammae/params.gamma_m, -params.p)*sec;
}else{
return 0.0;
}
}
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;
}
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\n", 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\n",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)。When设置为:
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;
型
程序以期望的根值结束:
% ./主要 * 7.957e+06 1.50212e+42
2条答案
按热度按时间svujldwt1#
你不能指望用
double
数学找到精确的值.即使是简单的方程3x=1
也不能用浮点数学精确求解(数字1/3不能用double
精确表示)。您的退出条件需要基于某个值。
IEEE 754双精度浮点数(通常是C++
double
的精度)只有53位尾数;因此应该考虑所涉及的值的大小来计算尾数值。如果说你的函数输出的值可以是“大”的(例如1 E50),那么一个单位顺序的整数是不够的(1 E50 +1 == 1 E50)。
还要注意的是,即使在中间计算中,浮点数的不准确也会造成问题,特别是当执行大的取消时(例如,当你减去两个非常大的数字时)。
正如Ted Lyngmo指出的那样,原始代码中也存在内存泄漏。
o4hqfura2#
你分配了太多的内存,它很可能会变得越来越慢,以寻找新的内存。如果你释放它时,完成它,它运行良好,并结束:
字符串
这是我对你的程序所做的唯一更改:
型