牛顿插值法的原理,在维基百科上不太全面,具体可以参考这篇文章。同样贴出,楼主作为初学者认为好理解的代码。
function p=Newton1(x1,y,x2) %p为多项式估计出的插值 syms x n = length(x1); %差商的求法 for i=2:n f1(i,1)=(y(i)-y(i-1))/(x1(i)-x1(i-1)); end for i=2:n for j=i+1:n f1(j,i)=(f1(j,i-1)-f1(j-1,i-1))/(x1(j)-x1(j-i)); end end f1=[y',f1]% 输出带0阶差商的差商表格 %Newton插值函数 Newton=f1(1,1); for i=2:n tt=1; for j=1:i-1 tt=tt*(x-x1(j)); end Newton=Newton+f1(i,i)*tt; end fprintf('Newton插值函数为 ') expand(Newton) % 将连乘多项式合并展开 x = x2; p = eval(Newton); % 代入值计算 fprintf('Newton插值函数在所求点x2的函数值为 ') p
运行:
输出:
CPP实现代码如下:
注意此处求差商运用的是另外一种方法
#include<iostream> #include<string> #include<vector> using namespace std; // 函数声明,使得其在被完整定义之前可以被引用 double ChaShang(int n, vector<double>&X, vector<double>&Y); double Newton(double x, vector<double>&X, vector<double>&Y); int main(){ int n; cout<<"输入插值点的个数:"<<endl; cin>>n; // 先将X,Y填充为n个0 vector<double>X(n,0); vector<double>Y(n,0); cout<<"请输入X[i],Y[i]:"<<endl; for(int i=0;i<n;i++){ cin>>X[i]>>Y[i]; } double x; cout<<"请输入要进行插值的点的x值:"<<endl; cin>>x; cout<<Newton(x,X,Y)<<endl; return 0; } // 此处为差商的另一种求法,可有差商定义根据数学归纳法求出 double ChaShang(int n,vector<double>&X,vector<double>&Y){ double f=0; double temp=0; for(int i=0;i<n+1;i++){ temp=Y[i]; for(int j=0;j<n+1;j++){ if(i!=j){temp /= (X[i]-X[j]);} } f += temp; }return f; } double Newton(double x,vector<double>&X,vector<double>&Y){ double result=0; for(int i=0;i<X.size();i++){ //此处的temp用于生成牛顿插值多项式里面的多项式乘积因子,(x-1)(x+3)这些 double temp=1; double f=ChaShang(i,X,Y); for(int j=0;j<i;j++){ temp=temp*(x-X[j]); } // 差商乘以因子得到最终的牛顿插值多项式 result += f*temp; }return result; }
运行结果: