具体思路现在,没空补,看代码之前,建议推一下newton插值公式。
#include <iostream>
#include <cstdio>
#include <queue>
#include <stdlib.h>
#include <ctime>
/// Newton 插值法的实现方法
#define NAN 0x3f3f3f3f
class Newdon{
// (x-x0)(x-x1).....(x-xn) *a
private:
double fac[1000]={0}; // o ~ n 的u多项式系数, 系数从0到n
double fac_ploy[1000]={0}; // 新增一项的(x-x0)(x-x1).....(x-xn)的各个多项式的系数,用于构造新的多项式系数
int fac_len = 0; //有效系数的长度 从0开始
double get_fac_ploy(double);// (x-x0)(x-x1).....(x-xn)的值
public:
Newdon(double *,double *,int);
void add_node(double,double);
double get_value(double x);
void test_print();
};
void Newdon::test_print()
{
for(int i=0;i<fac_len;i++)
printf(" %lf ",fac[i]);
printf("\n");
for(int i=0;i<fac_len;i++)
printf(" %lf ",fac_ploy[i]);
printf("\n");
}
Newdon::Newdon(double *x_value,double *y_value, int len)
{
if(len<1) printf("len is too small to creat Newton\n");
// 系数构造方式 比较巧妙
//动态构造 新建新newton项式的系数, 当然newton项式之间的系数有递推关系
for(int i=0;i<len;i++)
{
//test_print();
add_node(x_value[i],y_value[i]);
}
}
// 已有一个节点的基础上 添加一个新节点
void Newdon::add_node(double x,double y)
{
//第一个节点
if(fac_len == 0)
{
fac[0] = y;
fac_ploy[0] = -x; fac_ploy[1] = 1;//下一个多项式的(x-x0)
fac_len ++;
return;
}
// 之后的节点
double tmp1 = get_value(x);
fac_len++;// 有效的fac_ploy 比当前fac多一个
double tmp2 = get_fac_ploy(x);
double ploy_a = (y - tmp1) / tmp2; // 求出系数a,之后将系数相乘相加即可
for(int i=0;i<fac_len;i++)
{
fac[i] += fac_ploy[i] * ploy_a;
}
//随之,增加有效的ploy系数
for(int i=fac_len-1;i>=0;i--)
{
fac_ploy[i+1] += fac_ploy[i];
fac_ploy[i] = fac_ploy[i] * (-x);
}
}
double Newdon::get_fac_ploy(double x)
{
if(fac_len <= 0){
printf("fac_len is <= 0\n");
return 0;
}
double res = fac_ploy[fac_len-1];
for(int i=fac_len-2;i>=0;i--)
res = res * x + fac_ploy[i];
return res;
}
double Newdon::get_value(double x)
{
if(fac_len <= 0){
printf("fac_len is <= 0\n");
return 0;
}
double res = fac[fac_len-1];
for(int i=fac_len-2;i>=0;i--)
res = res * x + fac[i];
return res;
}
int main(){
srand(time(0));
double y[10] = {10,20,49,50,40,0,-40,-30,-10,20};
double x[10] = {1,2,3,4,5,6,7,8,9,10};
int len = 10;
/* for(int i=0;i<len;i++)
{
x[i] = rand()%10000;
y[i] = rand()%20000;
}
*/
Newdon first = Newdon(x,y,len);
first.test_print();
int a = 11, b =15;
for(int i=a;i<b;i++)
printf("c(%d) = %.2lf\n",i,first.get_value(i));
}
可用MATLAB的Newton的插值法源代码进行验证!
%newton.m
%求牛顿插值多项式、差商、插值及其误差估计的MATLAB主程序
%输入的量:X是n+1个节点(x_i,y_i)(i = 1,2, ... , n+1)横坐标向量,Y是纵坐标向量,
%x是以向量形式输入的m个插值点,M在[a,b]上满足|f~(n+1)(x)|≤M
%注:f~(n+1)(x)表示f(x)的n+1阶导数
%输出的量:向量y是向量x处的插值,误差限R,n次牛顿插值多项式L及其系数向量C,
%差商的矩阵A
function[y,R,A,C,L] = newton(X,Y,x,M)
n = length(X);
m = length(x);
for t = 1 : m
z = x(t);
A = zeros(n,n);
A(:,1) = Y';
s = 0.0; p = 1.0; q1 = 1.0; c1 = 1.0;
for j = 2 : n
for i = j : n
A(i,j) = (A(i,j-1) - A(i-1,j-1))/(X(i)-X(i-j+1));
end
q1 = abs(q1*(z-X(j-1)));
c1 = c1 * j;
end
C = A(n, n); q1 = abs(q1*(z-X(n)));
for k = (n-1):-1:1
C = conv(C, poly(X(k)));
d = length(C);
C(d) = C(d) + A(k,k);%在最后一维,也就是常数项加上新的差商
end
y(t) = polyval(C,z);
R(t) = M * q1 / c1;
end
L = poly2sym(C);
运行
M = 10000000000000
X = [1 2 3 4 5 6 7 8 9 10 ]
Y = [ 10 20 49 50 40 0 -40 -30 -10 20 ]
x = [11,12,13,14]
newton(X,Y,x,M)