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

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

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

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

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

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

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

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

由公式⑦③⑤⑥求得，
∑==n
k k x Q d 0
2
1
1)(
d
x Q y q n
k k
k
1
1
1
)
(∑==
d
x x k k
k
1
1
1
∑==
α d
d
11
=β
最后将
)(1
1
x Q q 展开累加到拟合多项式①中，有
a
t q a t q a 1
1
10
10
⇒⇒+
(3)对于j=2,3,……,m,逐步递推Q
j
（x ）
根据递推公式④有
Q
j
（x ）=（x-
α
1
-j ）
Q j 1
-（x ）-
β1
-j Q
j 2
-（x ）
=（x-α
1
-j ）(t
j 1
-x
j 1
-+….+
t t x 01+)-
β
1
-j (
b j 2
-x
j 2
-+….+
b b x 0
1
+)
假设
Q
j
（x ）=
x s j
j
+
x
s j j 1
1
--+…..+
x s 1
+s
则可以得到计算
s
k
（k=0,1,…,j ）的公式：
t
s j j
1
-=
t t s
j j j j 2111
----+-=α
,1
11b t t s k j k k j k
βα----+-= k=j-2,….,2,1
,0
1
1010
b t t s
j k j β
α----+-=
然后分别根据 ⑦式③式⑤式与⑥式计算下列量：
)(02
x Q d
k n
k j
j
∑==
)(0
x Q y q k
J
n
k k
j
∑==/d
j
)()(2
x Q x x k j
k n
k k j
∑==α
/d j
d
d j j
j
1
-=
β
再将
)(x Q
q J
j
项展开后累加到拟合多项式①中，即有
,a s
q a k k
j
k
=>+ k=j-1,….,1,0
a s
q j j
j
=>
最后，，为了便于循环使用向量B,T,与S,应将向量T 传递给B ，向量S 传递送给T,即 b t
k k
=>，k=j-1,…,1,0
t s
k k
=>，k=j,…,1,0
在实际计算中，为了防止运算溢出，将
x
k
用
x
k
—
x -
来代替，其中
x -
=∑=+n k k
x n 011
在这种情况下，拟合多项式的形式为 )()(10x x a a p
x m
-
-+=+)
(
2
2x x a -
-+…..+
)
(x x 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
x Q d
k n
k j
j
∑==
表示。

问题2：
Q
j
（x ）的迭代问题
解决办法：分别用三个向量B,T,S 分别存在
)(1
x Q
j -，)(x Q j
，)(1
x Q
j +的系数，用
b t
k k
=>， t s
k k
=>，依次迭代
问题3：系数拟合问题
解决办法：将上一次正交多项式次数相等的对应系数加到下一次的系数， 即: ,a s
q a k k
j
k
=>+ k=j-1,….,1,0
a s
q j j
j
=>
问题4：
Q
j
（
x
k
）的算法
解决方法：
Q
j
（
x
k
）是一个j 次多项式，从j 次到1次，递减乘x 来增加次数，再加上
s[k].
Q
j
（x i
）= (s s s
j j j
x x 21)(--++)
Q
j
（
x
i
）= (s s s
x s j j j j x x 321
2
)(---+++)
……..
Q
j （x）=x
s j j+x
s j
j
1
1
-
-
+…..+ x
s1+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");
}。