最小二乘法拟合多项式

　本科实验报告

　课程名称：

　 计算机数值方法 B

　 实验项目：

　最小二乘法拟合多项式

　 实验地点：

　 逸夫楼 302

　 专业班级：

　软件 1127 班

　 学号：

　 学生姓名：

　 指导教师：

　2013 年 4 月

　21

　 日

　 一、

　实验目的和实验要求：

　1．熟练运用已学计算方法求解方程组 2．加深对计算方法技巧，选择正确的计算方法来求解各种方程组 3．培养使用电子计算机进行科学计算和解决问题的能力 二、实验内容和原理 1、实验内容：

　给定数据点（x i

　,y i ），用最小二乘法拟合数据的多项式，并求平方误差。

　x i

　0 0.5 0.6 0.7 0.8 0.9 1.0 y i

　1 1.75 1.96 2.19 2.44 2.71 3.00 2、实验原理：

　建立正规方程组：

　∑（∑x ij+k ）ak =∑x ij yi

　,j=0,1,…,n 平方误差：

　I=∑（∑a k x ik -yi ）

　三、主要仪器设备 笔记本电脑，

　四、实验源程序、 实验结果及分析 实验源程序：

　#include<iostream.h> #include<fstream.h> #define N 15

　double power(double &a,int n) {

　 double b=1;

　for(int i=0;i<n;i++)

　 b*=a;

　return b; }

　void Gauss(); double X[N],Y[N],sumX[N],sumY[N],a[N][N],b[N],l[N][N],x[N];

　void main() {

　ofstream outdata;

　ifstream indata;

　double s;

　int i,j,k,n,index;

　cout<<"请输入已知点的个数 n=";

　cin>>n;

　cout<<endl;

　cout<<"请输入 X 和 Y:"<<endl;

　 //输入给定数据

　 for(i=0;i<n;i++) {

　 cout<<"X["<<i<<"]=";

　 cin>>X[i];

　 sumX[1]+=X[i];

　cout<<"Y["<<i<<"]=";

　 cin>>Y[i];

　 sumY[1]+=Y[i];

　 cout<<endl;

　}

　cout<<"sumX[1]="<<sumX[1]<<"\t"<<"sumY[1]="<<sumY[1]<<endl;

　cout<<"请输入拟合次数 index=";

　cin>>index;

　cout<<endl;

　i=n;

　sumX[0]=i;

　 for(i=2;i<=2*index;i++) {

　 sumX[i]=0;

　 for(j=0;j<n;j++)

　sumX[i]+=power(X[j],i);

　 cout<<"sumX["<<i<<"]="<<sumX[i]<<endl;

　 }

　 for(i=2;i<=index+1;i++) {

　 sumY[i]=0;

　 for(j=0;j<n;j++)

　sumY[i]+=power(X[j],i-1)*Y[j];

　 cout<<"sumY["<<i<<"]="<<sumY[i]<<endl;

　}

　for(i=1;i<=index+1;i++) {

　 //建立正规方程组

　 for(j=1;j<=index+1;j++)

　a[i][j]=sumX[i+j-2];

　 b[i]=sumY[i];

　}

　k=1;

　//用高斯消元法解方程组

　do{

　 for(j=k+1;j<=index+1;j++) l[j][k]=a[j][k]/a[k][k];

　 for(i=k+1;i<=index+1;i++){

　for(j=k+1;j<=index+1;j++)

　 a[i][j]=a[i][j]-l[i][k]*a[k][j];

　b[i]=b[i]-l[i][k]*b[k];

　 }

　 if(k==index+1) break;

　 k++;

　}while(1);

　x[index+1]=b[index+1]/a[index+1][index+1];

　 for(i=index;i>=1;i--) {

　 s=0;

　 for(j=i+1;j<=index+1;j++)

　s=s+a[i][j]*x[j];

　 x[i]=(b[i]-s)/a[i][i];

　 }

　cout<<"拟合系数为:";

　 //输出拟合系数

　for(i=1;i<=index+1;i++)

　 cout<<x[i]<<"\t";

　double m=0;

　cout<<endl<<"平方误差为:";

　for(i=0;i<n;i++) {

　double t=x[1]+x[2]*X[i]-Y[i];

　 m=m+power(t,2);

　}

　cout<<m<<endl;

　} 运行结果：