基于C++的拉格朗日插值实现问题排查（f(x)=(x²+1)^-1）

问题说明 #

• 目标函数：$f(x) = (x^2 + 1)^{-1}$，插值区间为$[-5,5]$
• 核心任务：
  1. 将区间按步长划分成201个点（h=200），计算每个点的$f(x)$值（此部分已验证正确）
  2. 根据用户输入的节点数n（可选5、11、21、41），计算每个划分点的拉格朗日插值多项式$P_n(x)$
• 异常表现：插值结果严重偏离原函数，例如n=5时x=4.95处，插值结果为46.1，而原函数值仅为0.039；多项式呈现指数增长后趋于平稳，未出现预期的Runge振荡，但误差极大，推测问题出在嵌套循环逻辑中。

问题代码 #

#include <stdio.h>
#include <math.h>
#include <iostream>
#include <iomanip>

using namespace std;
//THIS PROGRAM USES LAGRANGE INTERPOLATION
int main ()
{
    float userN;
    float deltaX, deltaXP, xp;
    float x[500], y[500],xP[50], yF[50], yP[500];
    float p, yp=0;
    float h=200;
    
    deltaX=10/h;
    x[0]=-5;
    y[0]=1/((x[0]*x[0])+1);
    for (int i=0; i<=h; i++)
    {
        x[i+1] = x[i]+deltaX;
        y[i+1] = 1/((x[i+1]*x[i+1])+1);
    }

    cout << "Select your desired nodes (degree n - you have 5, 11,21, and 41 to choose from):"<<endl;
    cin >> userN;
    if (userN==5 || userN==11 || userN==21 || userN==41)
    {
        deltaXP = 10/userN;
        xP[0]=-5;
        yF[0] = 1/((xP[0]*xP[0])+1);
        for (int i=0; i<=userN; i++)
        {
            xP[i+1] = xP[i]+deltaXP;
            yF[i+1] = 1/((xP[i+1]*xP[i+1])+1);
        }//TEST
            cout <<i<<" th x-value for p is: "<<xP[i]<<endl;
            cout <<i<<" th y-value for p is: "<<yF[i]<<endl;

        for (int i = 0; i<=h; i++)
        {
            xp=x[i];
            for (int j=0; j<=userN; j++)
            {
                p = 1;
                for (int k=0; k<=userN; k++)
                {
                    if (xP[j] != xP[k])
                    {
                        {
                            p = p* (xp - xP[k])/(xP[j] - xP[k]);
                        }
                    }
                }
                yp = yp+p*yF[j];
            }
            yP[i]=yp;
        }

        cout <<"*****The P(x_i) values for P are below*****"<<endl;
        for (int i=0; i<=h; i++)
        {
            cout <<yP[i]<<endl;
        }//use the similar loop to cout x and f(x) values
    }else {cout <<"You entered the wrong number of n. Exiting the program"<<endl;}

    return 0;
}

错误分析与修正方案 #

核心错误点 #

1. yp未重置：计算每个插值点的$P_n(x_i)$时，yp变量未在循环开头重置为0，导致每次计算都累加了上一个点的插值结果，最终数值爆炸，完全偏离正确值。
2. 浮点比较隐患：用if (xP[j] != xP[k])直接比较浮点数是否相等，可能因精度问题导致错误判断，改用索引判断k != j更可靠（等间距节点下，索引不同则节点值必然不同）。
3. 循环边界越界：原循环i<=h会生成超出区间的点，且节点生成循环也存在类似越界风险，导致数组元素访问异常。

修正后的代码 #

#include <stdio.h>
#include <math.h>
#include <iostream>
#include <iomanip>

using namespace std;
// 拉格朗日插值程序
int main ()
{
    int userN; // 用int存储节点数，避免浮点比较精度问题
    float deltaX, deltaXP, xp;
    float x[201], y[201], xP[50], yF[50], yP[201]; // 匹配实际需要的数组长度
    float p, yp;
    const int h = 200; // 用常量定义划分次数

    // 生成原函数的采样点
    deltaX = 10.0f / h;
    x[0] = -5.0f;
    y[0] = 1.0f / ((x[0] * x[0]) + 1.0f);
    for (int i = 0; i < h; i++)
    {
        x[i+1] = x[i] + deltaX;
        y[i+1] = 1.0f / ((x[i+1] * x[i+1]) + 1.0f);
    }

    // 输入节点数
    cout << "选择节点数（可选5、11、21、41）：" << endl;
    cin >> userN;
    if (userN == 5 || userN == 11 || userN == 21 || userN == 41)
    {
        // 生成插值节点
        deltaXP = 10.0f / userN;
        xP[0] = -5.0f;
        yF[0] = 1.0f / ((xP[0] * xP[0]) + 1.0f);
        for (int i = 0; i < userN; i++)
        {
            xP[i+1] = xP[i] + deltaXP;
            yF[i+1] = 1.0f / ((xP[i+1] * xP[i+1]) + 1.0f);
        }

        // 输出节点信息
        cout << "生成的插值节点：" << endl;
        for (int i = 0; i <= userN; i++)
        {
            cout << i << " 号节点：x=" << xP[i] << "，y=" << yF[i] << endl;
        }

        // 计算每个采样点的插值结果
        for (int i = 0; i <= h; i++)
        {
            xp = x[i];
            yp = 0.0f; // 每次计算新点前重置累加器
            for (int j = 0; j <= userN; j++)
            {
                p = 1.0f;
                for (int k = 0; k <= userN; k++)
                {
                    if (k != j) // 用索引判断代替浮点比较
                    {
                        p = p * (xp - xP[k]) / (xP[j] - xP[k]);
                    }
                }
                yp += p * yF[j];
            }
            yP[i] = yp;
        }

        // 输出插值结果与原函数值对比
        cout << "\n*****插值结果与原函数值对比*****" << endl;
        for (int i = 0; i <= h; i++)
        {
            cout << "x=" << fixed << setprecision(2) << x[i] 
                 << " | P(x)=" << fixed << setprecision(6) << yP[i] 
                 << " | f(x)=" << fixed << setprecision(6) << y[i] << endl;
        }
    }
    else 
    {
        cout << "输入的节点数无效，程序退出。" << endl;
    }

    return 0;
}

修正说明 #

• 重置累加器：在每个插值点的计算循环开头，将yp设为0，确保每个点的插值结果独立计算，不会累加之前的数值。
• 索引判断替代浮点比较：避免因浮点数精度问题导致的基函数计算错误。
• 修正循环边界：调整循环条件，避免数组越界，生成准确的区间内采样点和插值节点。
• 输出优化：增加原函数值的输出，方便对比插值结果与真实值的差异，同时格式化输出精度，提升可读性。

内容的提问来源于stack exchange，提问作者raspberry_flapper