在 C++ 故障中集成 Goempertz 函数

Question

我试图找到 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;
}

Answer 1

据我了解，您在常量 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;
}

Answer 2

问题在于您选择的结束 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 的最佳值，过程大致如下：

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

比如我选择的B和N就非常保守了。这些可以减少到 B = 50 和 N = 10，并且仍然对 3 个有效数字给出相同的结果。

在 C++ 故障中集成 Goempertz 函数

Integrating Goempertz Function in C++ glitch

c++

loops

function

numerical-integration