c++ - 为什么我得到 NaN 作为计算结果？ (拉格朗日标准形式的分段线性插值)

Question

我正在编写一个程序，它将以分段线性方式利用拉格朗日标准形式来插值 n 次多项式。我的代码在第一个子区间以及第三个和第四个子区间正常工作，但是，出于某种我无法理解的原因，我收到 NaN 作为第二个子区间的输出。第二个子区间在注释行//P2 的下方计算。我绞尽脑汁并尝试了所有我能想到的改变来解决这个问题，但没有运气。如果有人能提供一些见解，我将不胜感激。请注意，我只包含了直到第二个插值多项式的代码，因为第三个和第四个以类似的方式跟随。我提前为我的代码的残酷性道歉。我是 C++ 的新手，还没有时间获得对于如此大的问题可能期望的优雅。再次感谢。

ofstream Outfile;
Outfile.open ("PiecewiseLagrange_D.dat");
double *P1 = new double[201]; //Polynomial 1
double *P2 = new double[201]; //Polynomial 2
double *P3 = new double[201]; //Polynomial 3
double *P4 = new double[201]; //Polynomial 4
double *x = new double[201]; //Interpolating points/x's
double *x1 = new double[(int)n+1]; //First subinterval/mesh/xi's
double *x2 = new double[(int)n+1]; //Second subinterval
double *x3 = new double[(int)n+1]; //Third subinterval
double *x4 = new double[(int)n+1]; //Fourth subinterval
double a, b; //interval end points
char func; //function selection
double xDifference1;
double xDifference2;

cout << "Enter an interval with integer end points (lesser value first)";
cin >> a >> b;

for (int i=0; i<=n; i++) //Initialize
{
    P1[i] = 0;
    P2[i] = 0;
    P3[i] = 0;
    P4[i] = 0;
    x[i] = 0;
    x1[i] = 0;
    x2[i] = 0;
    x3[i] = 0;
    x4[i] = 0;
}


x1[0] = a;
for (int i=0; i<=n; i++)
{
    x1[i] = x1[0] + i*(((b-a)/4)/n);
    cout << x1[i] << endl;
}
cout << endl;
x2[0] = x1[(int) n];
for (int i=0; i<=n; i++)
{
    x2[i] = x2[0] + i*(((b-a)/4)/n);
    cout << x2[i] << endl;
}
cout << endl;
x3[0] = x2[(int) n];
for (int i=0; i<=n; i++)
{
    x3[i] = x3[0] + i*(((b-a)/4)/n);
    cout << x3[i] << endl;
}
cout << endl;
x4[0] = x3[(int) n];
for (int i=0; i<=n; i++)
{
    x4[i] = x4[0] + i*(((b-a)/4)/n);
    cout << x4[i] << endl;
}


cout << "Enter a function to evaluate (1,2, or 3):";
cin >> func;
//cout << "Polynomial is g1(x) on [" << a << "," << b << "]" << endl;

if (func == '1')
{
    //P1
    x[0] = a;
    for (int i=0; i<=200; i++)
    {
        x[i] = x[0] + i*((x1[(int) n] - x1[0])/200);
    }

    for (int j=0; j<=200; j++)
    {
        xDifference1 = 0;
        xDifference2 = 0;
        for (int i=0; i<=n; i++)
        {
            xDifference1 = (x1[i] - x1[i+1]);
            xDifference2 = (x1[i+1] - x1[i]);

            P1[j] = F1(x1[i])*((x[j] - x1[i+1])/xDifference1) + F1(x1[i+1])*((x[j] - x1[i])/xDifference2);
        }

        Outfile << x[j] << "      " << P1[j] << "     " << F1(x[j]) << endl;
        cout << setw(8) << x[j] << setw(12) << P1[j] << endl;
    }

    cout << endl;


    //P2
    x[0] = x1[(int) n];
    for (int i=0; i<=200; i++)
    {
        x[i] = x[0] + i*((x2[(int) n] - x2[0])/200);
    }

    for (int j=0; j<=200; j++)
    {
        xDifference1 = 0;
        xDifference2 = 0;
        for (int i=0; i<=n; i++)
        {
            xDifference1 = (x2[i] - x2[i+1]);
            xDifference2 = (x2[i+1] - x2[i]);

            P2[j] = F1(x2[i])*((x[j] - x2[i+1])/xDifference1) + F1(x2[i+1])*((x[j] - x2[i])/xDifference2);
        }

        Outfile << x[j] << "      " << P2[j] << "     " << F1(x[j]) << endl;
        cout << setw(8) << x[j] << setw(12) << P2[j] << "      " << F1(x[j]) << endl;
    }

    cout << endl;

Answer 1

由于您没有调用任何其他函数，因此当您将零除以零 (0.0/0.0) 时，您将得到一个 NaN。在某些时候，您的 xDifference1 和/或 xDifference2 为零。

将非零值除以零得到无穷大。

编辑 然而，显然情况并非如此，进一步的调查显示各种x 数组，包括x2，有n +1 个元素，索引为 0 到 n。在循环期间，您访问 x2[i+1]。由于 i 在最后一次迭代中将等于 n，因此您访问超出数组边界的元素 x2[n+1] 和结果在未定义的行为中。在这种情况下，数组后面的随机内存为 x2 而不是其他数组生成 NaN。

在一个不相关的注释中，你为每次迭代分配给 P2[j] 的内部 i 循环，因此你从循环中获得的唯一值是最后一次迭代。您是想使用 P2[j] += ... 吗？