实验题目用正交多项式做小二乘曲线拟合.pdf
实验题目： 用正交多项式做小二乘曲线拟合 实验题目：用正交多项式做最小二乘的曲线拟合 学生组号： 6 完成日期： 2011/11/27 1 实验目的 针对给定数据的煤自燃监测数据中煤温与 NO 22 ,之间的非线性关系，用正交多项式 做最小二乘曲线拟合。 2 实验步骤 2.1 算法原理 设给定n+1 个数据点：（ y xk k ,） ，k=0,1,· · ·,n，则根据这些节点作一个m 次的最 小二乘拟合多项式 pm（x）=a0+ xaxaa m m x. 2 21 = xa j m j j 0 其中， mn,一般远小于n.。 若要构造一组次数不超过的在给定点上正交的多项式函数系)(x Q j （， ， ，） ，则可以首先利用)(x Q j （， ，） 作为基函数作最小二乘 曲线的拟合，即 pm （x）)(.)()( 1100 xxx QqQqQq mm 根据式，其中的系数 q j （， ，）为 n k k j n k k jk j x Q x Qy q 0 2 0 )( )( ， ， 将代入后展开就成一般的多项式。 构造给定点上的正交多项式)(x Q j （， ，）的递推公式如下： 1,.,2, 1),()()()( )()( 1)( 11 0 1 0 mjxxxx xx x QQQ Q Q jjj j j 其中 j d x Q x j n k k j k 0 2 )( ， ， j d d j j 1 ， ， n k k j j xQd 0 2 )(， ， 则实际计算过程中，根据式逐步求出个正交多项式)(x Qj ，并用公式计算出 qj ，并将 每次计算展开后累加到拟合多项式中。 2.2 算法步骤 用三个向量 , ,存放多项式)( 1 x Q j ，)(x Q j ，)( 1 x Q j 的系数。 （）构造)( 0 x Q ，设)( 0 x Q b0 ，根据式，得 b0 。再根据式，计算： d0 d y q n k k 0 0 0 d x n k k 0 0 0 最后将)( 00 x Qq 项展开后累加到拟合多项式中，则 q0 b0 a0 （）构造)( 1 x Q ，设)( 1 x Q t0 t1 x，根据递推公式，则可知， t0 0 ，t 1 。由公式求得， n k k x Q d 0 2 1 1 )( d x Qy q n k k k 1 0 1 1 )( d x Q x n k kk 1 0 2 1 1 )( d d 0 1 1 最后将)( 11 x Qq 展开累加到拟合多项式中，有 atq atqa 11 1 00 1 0 (3)对于 j=2,3, ,m,逐步递推 Q j （ x） 根据递推公式有 Q j （x）=（x- 1j ）Q j1 （x）- 1j Q j2 （x） =（x- 1j ） (t j 1x j 1 +.+ tt x 01 )- 1j (b j 2x j2 + .+ bb x 01 ) 假设 Qj （x）= xs j j + xs j j 1 1 + + xs1 + s0 则可以得到计算 sk （k=0,1,j）的公式： ts jj1 tts jjjj2111 , 1 11 btts k j kkjk k=j-2, .,2,1 , 0 1 1010 bttsj kj 然后分别根据式式式与式计算下列量： )( 0 2 x Q d k n k j j )( 0 xQyq k J n k kj /d j )()( 2 0 x Q xx k j k n k kj /d j d d j j j 1 再将)(x Qq Jj 项展开后累加到拟合多项式中，即有 , asqakkjk k=j-1, .,1,0 as qjj j 最后，为了便于循环使用向量B,T,与 S,应将向量 T 传递给 B，向量 S传递送给 T,即 bt kk ，k=j-1, ,1,0 ts kk ，k=j, ,1,0 在实际计算中，为了防止运算溢出，将 xk 用 xk x 来代替，其中 x = n k k x n 0 1 1 在这种情况下，拟合多项式的形式为 )()( 10 xxaa p x m + )( 2 2 xx a + )( xx a m m 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.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 氧气数据： 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.15 80 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 分别 Q j （x）之间的关系 解决办法：观察公式，找出关系，将公式分解， 分步求出分子分母， 将分母用)( 0 2 xQd k n k j j 表示。 问题 2： Q j （x）的迭代问题 解决办法：分别用三个向量B,T,S 分别存在)( 1 x Q j ，)(x Qj ，)( 1 x Qj 的系数，用 bt kk ， ts kk ，依次迭代 问题 3：系数拟合问题 解决办法：将上一次正交多项式次数相等的对应系数加到下一次的系数， 即: , as q a kk j k k=j-1, .,1,0 as qjj j 问题 4： Q j （ xk ）的算法 解决方法： Q j （ xk ）是一个j 次多项式，从j 次到 1 次，递减乘x 来增加次数，再加上 sk. Q j （ xi ） = ( sss jjj xx 21) () Q j （ xi ） = ( sssxs jjjj xx 321 2 )() Q j （x） = xs j j + xs j j 1 1 + + 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,s20,t20,b20; for(i=0;in) m=n; if(m20) m=20; z=0.0; for(i=0;i1) /构造 Q1(x) t1=1.0; t0=-p; d2=0.0; c=0.0; g=0.0; for(i=0;i=3) for(k=j-2;k=1;k-) sk=-p*tk+tk-1-q*bk; s0=-p*t0-q*b0; d2=0.0;c=0.0;g=0.0; for(i=0;i=0;k-) q=q*(xi-z)+sk; d2=d2+q*q;c=c+yi*q; g=g+(xi-z)*q*q; c=c/d2;p=g/d2;q=d2/d1; d1=d2; aj=c*sj; tj=sj; for(k=j-1;k=0;k-) ak=c*sk+ak; bk=tk; tk=sk; /计算误差的平方和、误差的绝对值之和与误差绝对值的最大值 dt0=0.0; dt1=0.0;dt2=0.0; for(i=0;i=0;k-) q=ak+q*(xi-z); p=q-yi; if(fabs(p)dt2) dt2=fabs(p); dt0=dt0+p*p; dt1=dt1+fabs(p); return z; void main() int i; double a8,dt3,z; double x24=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 y124=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 y224=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(“%lft“,xi); printf(“n 输出氧气数据 :n“); for(i=0;i=23;i+) printf(“%lft“,y1i); z=spir(x,y1,24,a,4,dt);/ 函数的调用来求拟合函数系数、误差 printf(“n“); for(i=0;i=3;i+) printf(“a(%2d)=%lfn“,i,ai); printf(“n“); for(i=0;i=2;i+) printf(“dt(%2d)=%lf “,i,dti);/ 误差的输出 printf(“nn“); printf(“ 输出氧气与温度拟合的正交多项式:nn“); printf(“p(x)=%lf“,a0); for(i=1;i=3;i+) printf(“+(%lf)*(x-%lf)%d“,ai,z,i); printf(“nn“); printf(“ 输出温度数据:n“); for(i=0;i=23;i+) printf(“%lft“,xi); printf(“n 输出氮气数据 :n“); for(i=0;i=23;i+) printf(“%lft“,y2i); z=spir(x,y2,24,a,4,dt);/ 函数的调用来求拟合函数系数、误差 printf(“n“); for(i=0;i=3;i+) printf(“a(%2d)=%lfn“,i,ai);/函数系数的输出 printf(“n“); for(i=0;i=2;i+) printf(“dt(%2d)=%lf “,i,dti);/ 误差的输出 printf(“nn“); printf(“ 输出氮气与温度拟合的正交多项式:nn“); printf(“y(x)=%lf“,a0); for(i=1;i=3;i+) printf(“+(%lf)*(x-%lf)%d“,ai,z,i); printf(“n“);
