您的位置:首页 > 其它

牛顿插值多项式

2015-10-18 00:14 381 查看
//题目要求:用牛顿差商公式进行插值
#include <iostream>
#include <cmath>
#include <cstring>
using namespace std;
#define numofx 20                //定义数据量的大小
struct data{
double x;
double fx;                   //定义进行运算的点(x,fx)
}pointdata[numofx];

//求差商f(x0,x1,x2...xn)的公式
long double chashang(int n){     //计算f[x0,x1,x2...xn]
if(n<=0){
return 0;
}

long double result=0;
for(int i=0;i<=n;i++){
//计算这0~n的计算值的和
long double resproduct=1;
for(int j=0;j<=n;j++){
if(j!=i){
resproduct*=(pointdata[i].x-pointdata[j].x);
}
}
long double temp=pointdata[i].fx/resproduct;
result+=temp;
}
return result;
}

int main(){
int numofpoint;
cout<<"将要输入的插值点的个数为:";
cin>>numofpoint;
cout<<"请依次输入x,fx的值:\n";
for(int i=0;i<numofpoint;i++){
cin>>pointdata[i].x>>pointdata[i].fx;
}

//打印输出牛顿差商公式 ,未简约
cout<<"由以上插值点计算出的牛顿插值多项式为:\n";
cout<<"f(x)="<<pointdata[0].fx;
for(int i=1;i<numofpoint;i++){
long double temp4=chashang(i);
if(temp4>=0){
cout<<"+"<<temp4<<"*";
}
else{
cout<<temp4<<"*";
}

for(int j=0;j<i;j++){
cout<<"(x-"<<pointdata[j].x<<")";
}
}

//计算某个x的牛顿插值多项式近似值
double x;
cout<<"\n\n请输入所求变量x的值:";
cin>>x;
cout<<"牛顿插值多项式近似结果为:";

long double res=pointdata[0].fx;        //牛顿差商公式的f(x0)项
for(int i=0;i<numofpoint-1;i++){        //牛顿插值多项式项f[x0,x1...xn](x-x0)(x-x1)...(x-xn) n>=1

long double temp=1;
for(int j=0;j<=i;j++){
temp=temp*(x-pointdata[j].x);
}

res=res+chashang(i+1)*temp;
}
cout<<res<<endl;
}


拉格朗日插值的程序:(这个程序中运用了文件输入的方式读入数据,当然上面的牛顿插值程序也可以运用这个方法,减轻输入数据时的繁琐;同时也可以将生成的牛顿插值公式通过文件输出出来,便于复制粘贴)

#include <iostream>
#include <fstream>
using namespace std;

class lagrange{

private:
int i,j,n;
double mult,sum,z;
double *x,*y;

public:
void interpolation(){
ifstream fin("a.txt");
fin>>n;
x=new double
;
y=new double
;
for(i=0;i<n;i++){
fin>>x[i]>>y[i];
}
fin.close();
cout<<"\n输入需要插值的点:";
cin>>z;
sum=0.0;

for(i=0;i<n;i++){
mult=1.0;
for(j=0;j<n;j++){
if(j!=i){
mult*=(z-x[j])/(x[i]-x[j]);
}
}
sum+=mult*y[i];
}
cout<<"\n插补的值="<<sum<<endl;

}

~lagrange(){
delete[] x,y;
}

};

int main(){
lagrange interp;
interp.interpolation();
}
内容来自用户分享和网络整理,不保证内容的准确性,如有侵权内容,可联系管理员处理 点击这里给我发消息
标签: