实验题目用正交多项式做小二乘曲线拟合

实验题目：用正交多项式做小二乘曲线拟合实验题目： 用正交多项式做最小二乘的曲线拟合 学生组号： 6 完成日期： 2011/11/27 1 实验目的针对给定数据的煤自燃监测数据中煤温与NO 22,之间的非线性关系，用正交多项式做最小二乘曲线拟合。

2 实验步骤2.1 算法原理设给定n+1个数据点：（yx kk,），k=0,1,···,n ，则根据这些节点作一个m 次的最小二乘拟合多项式pm（x ）=a+x a x a a mm x +++ (2)21=x a jmj j ∑=0①其中，m ≤n,一般远小于n.。

若要构造一组次数不超过ｍ的在给定点上正交的多项式函数系｛)(x Qj（ｊ＝０，１，．．．，ｍ）｝，则可以首先利用｛)(x Qj（ｊ＝０，１，．．．，ｍ）｝作为基函数作最小二乘曲线的拟合，即pm（x ）＝)(...)()(11x x x Qq Q q Q q mm+++ ②根据②式，其中的系数qj（ｊ＝０，１，．．．，ｍ）为∑∑===nk kjnk kjkjx Q x Q y q2)()(，ｊ＝０，１，．．．，ｍ ③将④代入③后展开就成一般的多项式。

构造给定点上的正交多项式)(x Qj（ｊ＝０，１，．．．，ｍ）的递推公式如下：⎪⎪⎩⎪⎪⎨⎧-=--=-==-+1,...,2,1),()()()()()(1)(11010m j x x x x x x x QQ Q Q Q j jj j j βαα ④其中αj＝dx x jk j=0，ｊ＝０，１，．．．，ｍ－１ ⑤βj＝dd j j1-，ｊ＝１，２，．．．，ｍ－１ ⑥∑==nk k jjx Q d2)(，ｊ＝０，１，．．．，ｍ－１ ⑦则实际计算过程中，根据⑤式逐步求出个正交多项式)(x Qj，并用公式④计算出q j，并将每次计算展开后累加到拟合多项式①中。

2.2 算法步骤用三个向量Ｂ,Ｔ,Ｓ,存放多项式)(1x Qj -，)(x Q j，)(1x Qj +的系数。

（１） 构造)(0x Q ，设)(0x Q ＝b，根据④式，得b＝１。

再根据⑦③⑤式，计算：d＝ｎ＋１dy q nk k000∑==dx nk k00∑==α最后将)(0x Q q 项展开后累加到拟合多项式中，则q 0b 0＝a（２）构造)(1x Q ，设)(1x Q ＝t＋t1x ，根据递推公式④，则可知，t 0＝－α0，t 1＝１。

由公式⑦③⑤⑥求得，∑==nk k x Q d 0211)(dx Q y q nk kk111)(∑==dx x k kk111∑==α dd11=β最后将)(11x Q q 展开累加到拟合多项式①中，有at q a t q a 111010⇒⇒+(3)对于j=2,3,……,m,逐步递推Qj（x ）根据递推公式④有Qj（x ）=（x-α1-j ）Q j 1-（x ）-β1-j Qj 2-（x ）=（x-α1-j ）(tj 1-xj 1-+….+t t x 01+)-β1-j (b j 2-xj 2-+….+b b x 01+)假设Qj（x ）=x s jj+xs j j 11--+…..+x s 1+s则可以得到计算sk（k=0,1,…,j ）的公式：ts j j1-=t t sj j j j 2111----+-=α,111b t t s k j k k j kβα----+-= k=j-2,….,2,1,011010b t t sj k j βα----+-=然后分别根据 ⑦式③式⑤式与⑥式计算下列量：)(02x Q dk nk jj∑==)(0x Q y q kJnk kj∑==/dj)()(2x Q x x k jk nk k j∑==α/d jdd j jj1-=β再将)(x Qq Jj项展开后累加到拟合多项式①中，即有,a sq a k kjk=>+ k=j-1,….,1,0a sq j jj=>最后，，为了便于循环使用向量B,T,与S,应将向量T 传递给B ，向量S 传递送给T,即 b tk k=>，k=j-1,…,1,0t sk k=>，k=j,…,1,0在实际计算中，为了防止运算溢出，将xk用xk—x -来代替，其中x -=∑=+n k kx n 011在这种情况下，拟合多项式的形式为 )()(10x x a a px m--+=+)(22x x a --+…..+)(x x a mm--2.3 程序流程图3 实验结果分析温度数据：51.1 52.8 54.8 57.2 58.3 62.7 65.2 67.7 70.6 73.5 75.7 78.680.8 84.8 87.8 89.5 92.1 96.4 103.1 112.5 120.8 134.7 152.4 159.1氧气数据：20.13 19.9 19.61 18.98 19.61 19.32 19.73 19.59 19.38 20.18 19.98 19.58 19.74 19.26 19.59 18.77 18.66 17.47 17.01 16.33 15.95 13.76 5.91 5.43 氮气数据：79.5 79.67 79.59 80.24 80 80.31 79.35 80.1 80.37 79.7 79.78 80.1580 80.46 80.09 80.9 80.99 82.06 82.33 83.02 83.31 84.84 90.85 93.96 实验结果：注：dt(0)为误差的平方和，dt(1)为误差的绝对值之和，dt(2)误差绝对值最大值。

氧气与温度的拟合曲线：氮气与温度的拟合曲线：4 实验结论这次实验我们的拟合曲线是选择三次拟合曲线，虽说更高次的拟合曲线可以达到更好的效果，但是由于在计算机计算的过程中舍入误差的存在，使得次数高的项系数为零。

由于数据是观测得来，而我们的误差最大不超过1.4，所以误差在允许范围内，故从误差以及各方面的考虑，我们最终选择三次来拟合曲线。

为了能够达到视觉上的效果，我们在实验结果处附上了用matlab 所拟合得到的曲线，从图中可以看出，随着温度的不断增加氧气的含量在降低，且降低率随温度的增加而增加，但是对于氮气却是随着温度的增加而增加，且增加率随温度的增加而增加。

由于温度的不断增加，经过一定的化学变换，放出氮气，同时消耗氧气，而且在温度相对高时（在一定的温度范围内），其化学反应的速率快，这就是对上述结果的一个解释。

5 问题归纳与总结在本次试验中，组员在分析问题的过程中主要遇到了一下几方面的问题，下面一一表述并给出解决办法。

问题一：αj，βj分别Qj（x ）之间的关系解决办法：观察公式，找出关系，将公式分解，分步求出分子分母，将分母用)(02x Q dk nk jj∑==表示。

问题2：Qj（x ）的迭代问题解决办法：分别用三个向量B,T,S 分别存在)(1x Qj -，)(x Q j，)(1x Qj +的系数，用b tk k=>， t sk k=>，依次迭代问题3：系数拟合问题解决办法：将上一次正交多项式次数相等的对应系数加到下一次的系数， 即: ,a sq a k kjk=>+ k=j-1,….,1,0a sq j jj=>问题4：Qj（xk）的算法解决方法：Qj（xk）是一个j 次多项式，从j 次到1次，递减乘x 来增加次数，再加上s[k].Qj（x i）= (s s sj j jx x 21)(--++)Qj（xi）= (s s sx s j j j j x x 3212)(---+++)……..Qj （x）=xs j j+xs jj11--+…..+ xs1+s0问题5：溢出问题解决方法：由于给定温度数据过大，次数大时会出现溢出，所以采用平移思想用x--x来表示x，图像不变。

参考文献李庆扬，王能超，易大义.数值分析.北京：清华大学出版社，2008附录（源代码）#include "stdio.h"#include "math.h"double spir(double x[],double y[],int n,double a[],int m,double dt[]) {int i,j,k;double z,p,c,g,q,d1,d2,s[20],t[20],b[20];for(i=0;i<=m-1;i++)//系数的初始化a[i]=0.0;if(m>n)m=n;if(m>20)m=20;z=0.0;for(i=0;i<=n-1;i++)z=z+x[i]/(1.0*n);//温度平均值b[0]=1.0;//Q0(x)多项式系数d1=1.0*n;p=0.0;c=0.0;for(i=0;i<=n-1;i++){p=p+(x[i]-z);c=c+y[i];}c=c/d1;p=p/d1;a[0]=c*b[0];if(m>1) //构造Q1(x){t[1]=1.0;t[0]=-p;d2=0.0;c=0.0;g=0.0;for(i=0;i<=n-1;i++){q=x[i]-z-p;d2=d2+q*q;c=c+y[i]*q;g=g+(x[i]-z)*q*q;}c=c/d2;p=g/d2;q=d2/d1;d1=d2;a[1]=c*t[1];a[0]=c*t[0]+a[0];}for(j=2;j<=m-1;j++) //构造Qj(x) {s[j]=t[j-1];s[j-1]=-p*t[j-1]+t[j-2];if(j>=3)for(k=j-2;k>=1;k--)s[k]=-p*t[k]+t[k-1]-q*b[k];s[0]=-p*t[0]-q*b[0];d2=0.0;c=0.0;g=0.0;for(i=0;i<=n-1;i++){q=s[j];for(k=j-1;k>=0;k--)q=q*(x[i]-z)+s[k];d2=d2+q*q;c=c+y[i]*q;g=g+(x[i]-z)*q*q;}c=c/d2;p=g/d2;q=d2/d1;d1=d2;a[j]=c*s[j];t[j]=s[j];for(k=j-1;k>=0;k--){a[k]=c*s[k]+a[k];b[k]=t[k];t[k]=s[k];}}//计算误差的平方和、误差的绝对值之和与误差绝对值的最大值dt[0]=0.0;dt[1]=0.0;dt[2]=0.0;for(i=0;i<=n-1;i++){q=a[m-1];for(k=m-2;k>=0;k--)q=a[k]+q*(x[i]-z);p=q-y[i];if(fabs(p)>dt[2])dt[2]=fabs(p);dt[0]=dt[0]+p*p;dt[1]=dt[1]+fabs(p);}return z;}void main(){int i;double a[8],dt[3],z;double x[24]={51.1,52.8,54.8,57.2,58.3,62.7,65.2,67.7,70.6,73.5,75.7,78.6,80.8,//温度数据84.8,87.8,89.5,92.1,96.4,103.1,112.5,120.8,134.7,152.4,159.1};double y1[24]={20.13,19.9,19.61,18.98,19.61,19.32,19.73,19.59,19.38,20.18,19.98,//氧气数据19.58,19.74,19.26,19.59,18.77,18.66,17.47,17.01,16.33,15.95,13.76,5.91,5.43};double y2[24]={79.5,79.67,79.59,80.24,80,80.31,79.35,80.13,80.37,79.7,79.78,80.15,80,80.46,80.09,80.9,80.99,82.06,82.33,83.02,83.31,84.84,90.85,93.96};//氮气数据printf("----------------正交多项式拟合-------------\n");printf("输出温度数据:\n");for(i=0;i<=23;i++)printf("%lf\t",x[i]);printf("\n输出氧气数据:\n");for(i=0;i<=23;i++)printf("%lf\t",y1[i]);z=spir(x,y1,24,a,4,dt);//函数的调用来求拟合函数系数、误差printf("\n");for(i=0;i<=3;i++)printf("a(%2d)=%lf\n",i,a[i]);printf("\n");for(i=0;i<=2;i++)printf("dt(%2d)=%lf ",i,dt[i]);//误差的输出printf("\n\n");printf("输出氧气与温度拟合的正交多项式:\n\n");printf("p(x)=%lf",a[0]);for(i=1;i<=3;i++)printf("+(%lf)*(x-%lf)^%d",a[i],z,i);printf("\n\n");printf("输出温度数据:\n");for(i=0;i<=23;i++)printf("%lf\t",x[i]);printf("\n输出氮气数据:\n");for(i=0;i<=23;i++)printf("%lf\t",y2[i]);z=spir(x,y2,24,a,4,dt);//函数的调用来求拟合函数系数、误差printf("\n");for(i=0;i<=3;i++)printf("a(%2d)=%lf\n",i,a[i]);//函数系数的输出printf("\n");for(i=0;i<=2;i++)printf("dt(%2d)=%lf ",i,dt[i]);//误差的输出printf("\n\n");printf("输出氮气与温度拟合的正交多项式:\n\n");printf("y(x)=%lf",a[0]);for(i=1;i<=3;i++)printf("+(%lf)*(x-%lf)^%d",a[i],z,i);printf("\n");}。