在C++故障中集成Goempertz功能

Integrating Goempertz Function in C++ glitch

本文关键字：Goempertz 功能 集成 C++ 故障      更新时间：2023-10-16

我试图找到Goempertz函数的梯形规则估计，并使用它来测量50岁吸烟者和50岁非吸烟者的预期寿命之间的差异，但我的代码一直给了我废话答案。

50岁人的Goempertz函数可以编码为：

exp((-b/log(c))*pow(c,50)*(pow(c,t)-1))

其中b和c是常数，我们需要将其从 0 积分到无穷大(一个非常大的数字)以获得预期寿命。

对于非吸烟者，预期寿命可以通过以下方式计算： 常量 b = 0.0005，c = 1.07。 对于吸烟者，预期寿命可以用 常量 b = 0.0010，c = 1.07。

const double A = 0; // lower limit of integration
const double B = 1000000000000; // Upper limit to represent infinity
const int N = 10000; //# number of steps of the approximation

double g(double b, double c, double t)  // 
{//b and c are constants, t is the variable of integration.
return exp((-b/log(c))*pow(c,50)*(pow(c,t)-1));
}
double trapezoidal(double Bconst, double Cconst)
{
double deltaX = (B-A)/N; //The "horizontal height" of each tiny trapezoid
double innerTrap = 0; //innerTrap is summation of terms inside Trapezoidal rule
for (int i = 0; i <= N; i++)
{
double xvalue;
if (i == 0) // at the beginning, evaluate function of innerTrap at x0=A
{
xvalue = A;
}
else if (i == N) //at the end, evaluate function at xN=B
{
xvalue = B;
}
else //in the middle terms, evaluate function at xi=x0+i(dX)
{
xvalue = A + i * deltaX;
}
if ((i == 0) || (i == N)) //coefficient is 1 at beginning and end
{
innerTrap = innerTrap + 1*g(Bconst, Cconst, xvalue);
}
else // for all other terms in the middle, has coefficient 2
{
innerTrap = innerTrap + 2*g(Bconst, Cconst, xvalue);
}
}
return (deltaX/2)*innerTrap;
}
int main()
{
cout << "years 50 year old nonsmoker lives: " << trapezoidal(0.0005,1.07) << endl;
cout << "years 50 year old smoker lives: " << trapezoidal(0.0010,1.07) << endl;
cout << "difference between life expectancies: " << trapezoidal(0.0005,1.07)-trapezoidal(0.0010,1.07) << endl;
return 0;
}

问题在于你选择的末端x坐标以及你对面积求和的切片数量：

const double A = 0;
const double B = 1000000000000;
const int N = 10000;
double deltaX = (B-A) / N;  //100 million!

当你做这样的离散集成时，你希望你的deltaX与函数的变化相比很小。我猜Goempertz函数在0到1亿之间变化很大。

要修复它，只需进行两项更改：

const double B = 100;
const int N = 10000000;

这使得deltaX == 0.00001并且似乎给出了良好的结果(21.2和14.8)。放大B不会改变最终答案(如果有的话)，因为此范围内的函数值基本上为 0。

如果您希望能够弄清楚如何选择B的良好值并N该过程大致如下：

  1. 对于B查找函数结果足够小(或函数更改足够小)以忽略的x的值。这对于周期性或复杂函数来说可能很棘手。
  2. 从一个小N值开始，然后计算结果。将N增加 2 倍(或更多)，直到结果收敛到所需的精度。
  3. 您可以通过增加B来检查您选择的是否有效，并查看结果的变化是否小于您想要的准确性。

例如，我对B和N的选择非常保守。这些可以减少到低至B = 50和N = 10，并且仍然为 3 个有效数字提供相同的结果。

据我了解，您在常量B和N上犯了一个错误。B- 一个人可以以一定的概率和N生活多少年是积分步。因此B应该相对较小(<100，因为一个人活50+100岁或更长时间的概率非常小)，N应该尽可能大。您可以使用以下代码来解决您的任务

const double A = 0; // lower limit of integration
const double B = 100; // Upper limit to represent infinity
const int N = 1000000; //# number of steps of the approximation
double g(double b, double c, double t)  // 
{//b and c are constants, t is the variable of integration.
return exp((-b/log(c))*pow(c,50)*(pow(c,t)-1));
}
double trapezoidal(double Bconst, double Cconst)
{
double deltaX = (B-A)/double(N); //The "horizontal height" of each tiny trapezoid
double innerTrap = 0; //innerTrap is summation of terms inside Trapezoidal rule
double xvalue = A + deltaX/2;
for (int i = 0; i < N; i++)
{
xvalue += deltaX;
innerTrap += g(Bconst, Cconst, xvalue);
}
return deltaX*innerTrap;
}
int main()
{
double smk = trapezoidal(0.0010,1.07);
double nonsmk = trapezoidal(0.0005,1.07);
cout << "years 50 year old nonsmoker lives: " << nonsmk  << endl;
cout << "years 50 year old smoker lives: " << smk << endl;
cout << "difference between life expectancies: " << nonsmk-smk << endl;
return 0;
}
相关文章：