最小二乘法 多项式拟合 实验报告

连续系统仿真实验报告实验数据拟合建模姓名：专业：学号：时间：2013年5月1日实验单元二 实验数据拟合建模一、实验目的1、用C 语言实现最小二乘的多项式拟合和LU 分解法；2、熟练掌握最小二乘拟合和LU 分解法的基本原理。

3、体会用计算机编程解决计算问题的方法。

二、需求说明(一)、需求阐述本次实验是要求根据已知的自变量和函数值，通过多项式拟合来分别计算2、3、4阶拟合多项式，并根据拟合结果分别计算出待求点的函数值。

其中解拟合系数方程组时采用LU 分解的方法计算拟合多项式的系数。

（二）、实验公式m 次拟合函数公式为：01x =++m m a a x a x ϕ()⋯计算系数i a 的方程组为：{0011++...+=m m s a s a s a t1021+11++...+=m m s a s a s a t ⋮0+112++...+=m m m m s a s a s a t○1其中 ++=0=nk jj k ii s x ∑, =0=nj i i i t y x ∑所以，在编程计算时，先计算出方程组○1，再用LU 分解法计算求出i a 的值，即可得到拟合多项式。

LU 分解法的公式为：[ 1l 21l 31 ⋮l n101l 32⋮l n2001⋮⋯⋯⋯⋱⋱⋯0 0⋮01] [u 110000 u 12 u 22000 u 13 u 23u 33⋱⋯⋯⋯⋱⋱0 u 1nu 2n ⋮ u n−1n u nn ] =[ a 11a 21a 31⋮a n1a 12a 22a 32⋮a n2a 13a 23a 33⋮a n3⋯⋯⋯⋱⋯a 1n a 2n a 3n ⋮a nn ]其中L 矩阵和U 矩阵的计算公式如下： 第一步，当k=1,有：{11= (=2,3,,n)j j u a j 1111111== (=2,3,,n)i i i a al i u a第二步，当k=2,3,⋯,n-1时，有：{-1=1=- (j=+1,+2,)k kjkj kr rj r u a l u k k n ∑-1(k)=1(k)-=(j=,+1,)k ikir rkr ik kka l u l k k n u∑最后求 u nn ：-1=1=-n nn nn nr rn r u a l u ∑三、设计说明（一）、数据结构程序采用一维数组的形式来读取文件中给出的已知点处的值和要计算的未知点处的自变量值，最终的拟合计算结果也是采用一维数组的形式输出到文件中。

拟合多项式的系数a 和拟合系数方程组的参数t 都是采用一维数组来存储的，而拟合系数方程组中的参数s 和L 、U 矩阵都是用二维数组来表示的。

由于要分别计算2、3、4阶拟合结果，所以数组的规模取为5，矩阵的规模取为5*5. （二）、算法设计及效率分析 在进行LU 分解函数中，在计算L 矩阵和U 矩阵时，因为当k=2,3,⋯,n-1时，计算-1=1k ir rk r l u ∑和-1=1k kr rj r l u ∑的循环条件不允许k=1时进入，而正好k=1时，计算1i l 和1j u 不需要-1=1k irrk r lu ∑和-1=1k kr rj r l u ∑，因而对k=1和k=2,3,⋯,n-1,就可以和在一起计算，这样就减少了程序的长度。

而在分别计算2、3、4阶拟合系数方程组的参数时，没有很好的利用前一阶计算的，而每次都要重新计算；而且矩阵是一个堆成矩阵，没有好好利用对称矩阵的特性，导致了重复计算，增加了计算量，降低了程序的效率。

而造成这一结果的原因是自己为了编程的简单而忽视了计算量，在以后的编程时要注意改变这一习惯。

（三）、程序结构程序主要步骤的流程图如下：以上流程图对应的源程序中的函数分别如下：//计算拟合系数方程组中的参数svoid computers(double s[p][p],double x1[],int m) //计算拟合系数方程组中的参数tvoid computert(double t[],double x1[],double y1[],int m) //对拟合系数方程组中的参数s 组成的矩阵进行LU 分解 void LV(double L[p][p],double V[p][p],double s[p][p],int m) //计算拟合多项式的系数void computera(double L[p][p],double V[p][p],double t[],double x[],int m) //由得到的的拟合多项式计算待求点处的函数值 void computerty2(double a[],double x2[],double y2[],int m) //保存得到的拟合多项式和计算处的参数void save(double a[],int m,double y2[])四、编码实现#include <iostream> #include <fstream> #include <string> #include <math.h>#define p 5 //拟合方程的阶次+1#define q 5 //已知点的数目，也是带计算点的数目 using namespace std;ofstream outDatay("G:\\连续系统仿真\\拟合实验\\outy.txt"); //用于保存计算结果 int main() {void computers(double s[p][p],double x1[],int m); //计算拟合方程组的系数svoid computert(double t[],double x1[],double y1[],int m); //计算拟计算出拟合系数方程组中的参数S 计算出拟合系数方程组中的参数t 对列出的拟合方程组阵进行LU 分解根据分解的结果计算拟合多项式的系数，即拟合方程组的解根据得到的拟合多项式计算待求点出函数值保存得到的拟合多项式系数和计算出的函数值合方程组的系数tvoid LV(double L[p][p],double V[p][p],double s[p][p],int m); //LU分解void computera(double L[p][p],double V[p][p],double t[],double x[],int m); //计算拟合多项式的系数void computerty2(double a[],double x2[],double y2[],int m);void save(double a[],int m,double y2[]);ifstream inDatax,outDatax,inDatay;int i;double x1[q],x2[q],y1[q],y2[q],s[p][p]={0},t[p]={0},a[p],L[p][p],V[p][p]; inDatax.open("G:\\连续系统仿真\\拟合实验\\inx.txt"); //已知点的自变量x值i=0;while(!inDatax.eof())inDatax>>x1[i++];inDatax.close();i=0;outDatax.open("G:\\连续系统仿真\\拟合实验\\outx.txt"); //要求的插值点的x值while(!outDatax.eof())outDatax>>x2[i++];outDatax.close();i=0;inDatay.open("G:\\连续系统仿真\\拟合实验\\iny.txt"); //已知点的因变量y值while(!inDatay.eof())inDatay>>y1[i++];inDatay.close();computers(s,x1,2);computert(t,x1,y1,2);LV(L,V,s,2);computera(L,V,t,a,2);computerty2(a,x2,y2,2);save(a,2,y2);computers(s,x1,3);computert(t,x1,y1,3);LV(L,V,s,3);computera(L,V,t,a,3);computerty2(a,x2,y2,3);save(a,3,y2);computers(s,x1,4);computert(t,x1,y1,4);LV(L,V,s,4);computera(L,V,t,a,4);computerty2(a,x2,y2,4);save(a,4,y2);return 0;}void computers(double s[p][p],double x1[],int m){int i,j;for (i=0;i<=m;i++){ for(j=0;j<=m;j++)s[i][j]=0;}for ( i=0;i<=m;i++)for (j=0;j<=m;j++){for(int k=0;k<q;k++){s[i][j]+=pow(x1[k],i+j);}}}void computert(double t[],double x1[],double y1[],int m) {int i;for (i=0;i<=m;i++){t[i]=0;}for (i=0;i<=m;i++)for (int j=0;j<q;j++){t[i]+=y1[j]*pow(x1[j],i);}}void LV(double L[p][p],double V[p][p],double s[p][p],int m) {double sum;for (int i=0;i<=m;i++)for(int j=0;j<=m;j++ ){L[i][j]=0;V[i][j]=0;}for(i=0;i<=m;i++){L[i][i]=1;}for (i=0;i<=m;i++){//计算U矩阵的第i行for (int j=0;j<=m;j++){ sum=0;for (int k=0;k<=i-1;k++){sum+=L[i][k]*V[k][j];}V[i][j]=s[i][j]-sum;}if(i<m){//计算L矩阵的第i列for (j=1;j<=m;j++){sum=0;for (int k=0;k<=i-1;k++){sum+=L[j][k]*V[k][i];}L[j][i]=(s[j][i]-sum)/V[i][i];}}}}void computera(double L[p][p],double V[p][p],double t[],double x[],int m) //计算拟合多项式的系数{double y[p],sum;for (int i=0;i<=m;i++){sum=0;for (int j=0;j<i;j++){sum+=L[i][j]*y[j];}y[i]=t[i]-sum;}for (i=m;i>=0;i--){sum=0;for (int j=m;j>i;j--){sum+=V[i][j]*x[j];}x[i]=(y[i]-sum)/V[i][i];}}void computerty2(double a[],double x2[],double y2[],int m)//由得到的的拟合多项式计算待求点处的函数值{double sum;for (int i=0;i<q;i++){sum=0;for (int j=0;j<=m;j++){sum=sum+a[j]*pow(x2[i],j);}y2[i]=sum;}}void save(double a[],int m,double y2[]) //保存得到的拟合多项式和计算处的参数{outDatay<<m<<"阶拟合模型的系数为："<<endl;for (int i=0;i<=m;i++){outDatay<<a[i]<<" ";}outDatay<<endl;outDatay<<"在待求点处的计算结果为："<<endl;for (i=0;i<q;i++){outDatay<<y2[i]<<" ";}outDatay<<endl;outDatay<<"---------------------------------------"<<endl;}四、实验及分析（一）实验情况2阶拟合模型的系数为：-0.0196404 1.14058 -0.278175在待求点处的计算结果为：0.091636 0.297498 0.481106 0.642461 0.781561---------------------------------------3阶拟合模型的系数为：-0.0013368 1.01202 -0.0330375 -0.136187在待求点处的计算结果为：0.0993986 0.295619 0.47939 0.644177 0.78344---------------------------------------4阶拟合模型的系数为：0.000276 0.99682 0.0132292 -0.192187 0.0233333在待求点处的计算结果为：0.0999004 0.295513 0.479428 0.644214 0.783334---------------------------------------而在上一个插值实验中计算出待求点的值分别为：0.0999004 0.295513 0.479428 0.644214 0.783334从以上计算结果中，我们可以看出4阶拟合模型的计算结果最接近原函数的真实值。