[拉格朗日插值]C实现求拉格朗日插值多项式系数
作者:互联网
题目
实现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);
}
参考
标签:拉格朗,...,插值,double,coeff,int,多项式 来源: https://www.cnblogs.com/umiu/p/15498395.html