c++ - 在 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

问题在于您选择的结束 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 位有效数字相同的结果。