题目

实现n次拉格朗日插值多项式L(x)的计算，插值函数原型为int lagPolynomial(int n, double* X, double* Y, double* a),其中X为插值节点数组$x_0$

实现

拉格朗日插值多项式可表示为，\(L(x)=\sum_{j=0}^{n} y_jl_j(x)\)

[1]中给出拉格朗日基函数的形式，考虑\(l_0(x)\)

\[l_0(x)=(\frac{x-x_1}{x_0-x_1})...(\frac{x-x_{n-1}}{x_0-x_{n-1}})(\frac{x-x_{n}}{x_0-x_{n}}) \]

\(l_0(x)\)共n项，n项分母的乘积是常数，将n项分母的乘积和\(y_0\)一起考虑为常数\(C_0\)。n项分子的乘积为\(P_j\)。则\(L(x)=\sum_{j=0}^nC_jP_j\)

总的框架

外层循环n+1次，内层循环n次， multiPolynomial函数时间复杂度O(N)，因此总的时间复杂度是O(N^3)。

void lagPolynomial(int n, double* X, double* Y, double* a){
    /* 
        ref: https://stackoverflow.com/questions/9860937/how-to-calculate-coefficients-of-polynomial-using-lagrange-interpolation

        拉格朗日插值多项式 n+1个数据点
        y_0 * (x-x_1)(x-x_2)...(x-x_n)/(x_0-x_1)(x_0-x_2)...(x_0-x_n) + 
        y_1 * (x-x_0)(x-x_2)...(x-x_n)/(x_1-x_0)(x_1-x_2)...(x_1-x_n) + 
        ...
    */

    double coeff[n+1] = {0};
    for (int i = 0; i <= n; i++)
    {
        /* 外层循环遍历基函数y_i*l_i */
        double l_coeff[n+1] = {0};
        l_coeff[0] = 1;
        
        int j;
        for (j = 0; j <= n; j++)
        {
            /* (x-x_j) j!=i */
            if (j == i)
                continue;
             
            /*  
                l_coeff初始化为1
                本次循环乘多项式x-x[j]，
            */
           multiPolynomial(n, l_coeff, -X[j]);
        }
        double denominator = calDenominator(n, X, Y, i);
        for (j = 0; j <= n; j++)
        {
            /* 多项式系数乘以 分母*yi */
            l_coeff[j] *= denominator;
            coeff[j] += l_coeff[j];
        }
    }

    for (int i = 0; i <= n; i++){
        a[i] = coeff[i];
    }
};

多项式的积

计算coeff表示的多项式乘以 x-constant。时间复杂度O(N)。

x使得coeff的所有系数右移一位，constant直接乘上每个系数。

void multiPolynomial(int n, double* coeff, int constant){
    /* coeff表示的多项式乘以 x-constant  x使得coeff的所有系数右移一位，constant直接乘上每个系数 */
    int pre = 0, tmp;
    for (int i = 0; i <= n; i++)
    { 
        tmp = coeff[i];
        coeff[i] = pre + coeff[i]*constant;
        pre = tmp;
    }
}

求常数

计算常数\(C_{idx}\)。时间复杂度O(N)。

double calDenominator(int n, double* X, double* Y, int idx){
    /* 
        计算分母，
        y_idx / (x_idx-x_1)(x_idx-x_2)...(x_idx-x_n) 
     */
    double sum = 1.0;
    for (int i = 0; i <= n; i++)
    {
        if (idx == i)
            continue;
        
        sum *= (X[idx]-X[i]);
    }
    return Y[idx] / sum;
}

测试

进行简单的测试

void show(double* nums, int n){
    for (int i = 0; i < n; i++)
    {
        std::cout << nums[i] << " ";
    }
    std::cout << std::endl;
}

void testCalDenominator(){
    /* 
        三个样本点(1,4), (2,7), (3,9)
        期望的结果
        C_0: 2
        C_1: -7
        C_2: 4.5
    */
    int n = 2;
    double X[n+1] = {1,2,3};
    double Y[n+1] = {4,7,9};
    double a[n+1] = {0};
    std::cout << " " << calDenominator(n, X, Y, 0)<<std::endl;
    std::cout << " " << calDenominator(n, X, Y, 1)<<std::endl;
    std::cout << " " << calDenominator(n, X, Y, 2)<<std::endl;
}
void testMultiPolynomial(){
    /* 
        (x-3)(x-5)=x^2-8x+15 
        [1,0,0]->[-3,1, 0]->[15,-8,1]
    */
   double coeff[3] = {1,0,0};
   show(coeff, 3);
   multiPolynomial(2, coeff, -3);
   show(coeff, 3);
   multiPolynomial(2, coeff, -5);
   show(coeff, 3);
}

void testLagPolynomial(){
    /* 
        (1,4), (2,7), (3,9)
        拉格朗日多项式：-0.5*x^2+4.5*x
        期望a输出: [0, 4.5, -0.5]
     */
    int n = 2;
    double X[n+1] = {1,2,3};
    double Y[n+1] = {4,7,9};
    double a[n+1] = {0};
    lagPolynomial(n, X, Y, a);
    show(a, n+1);
}

参考

1. javascript - How to calculate coefficients of polynomial using Lagrange interpolation - Stack Overflow

标签：拉格朗,...,插值,double,coeff,int,多项式
来源： https://www.cnblogs.com/umiu/p/15498395.html