实验题目用正交多项式做小二乘曲线拟合.docx
《实验题目用正交多项式做小二乘曲线拟合.docx》由会员分享，可在线阅读，更多相关《实验题目用正交多项式做小二乘曲线拟合.docx（17页珍藏版）》请在冰豆网上搜索。
实验题目用正交多项式做小二乘曲线拟合
实验题目：
用正交多项式做小二乘曲线拟合
学生组号：
_6_完成日期：
2011/11/27
1实验目的
针对给定数据的煤自燃监测数据中煤温与O2,之间的非线性关系，用正交多项式做最小二乘曲线拟合。
2实验步骤
2.1算法原理
设给定n+l个数据点：
（血，）［），k=0,l,・・・,ii,则根据这些节点作一个m次的最
小二乘拟合多项式化（x）=偽+卧+0才+・・・+6?
』”=工°/‘①
7=0
其中，—般远小于n.。
若要构造一组次数不超过m的在给定点上正交的多项式函数系｛。
‘切（j=0,
1,...,m）｝,则可以首先利用｛g（a）（j=0,1m）｝作为基函数作最小二乘
曲线的拟合，即p（X）=Qo（x）+q^Q（x）+...+qQ（x）②根据②式，其中的系数q.（j=o,1,...,m）为
f儿0（儿）
j=0,1八・・,m
k=OJ
将④代入③后展开就成一般的多项式。
构造给定点上的正交多项式2（X）（j=0,
1m）的递推公式如下:
0。
（归
迄⑴弋-久）④
0+1（x）=（x~a）Ql（x）-p.（x），j=1,2,•••，用-1
其中
刃2
—k=0J
d,
m—1
ar
d广丈2（Q'j=0‘I，…’m—l
k=0J
则实际计算过程中，根据⑤式逐步求出个正交多项式Q\x）,并用公式④计算出q,并将每次计算展开后累加到拟合多项式①中。
2.2算法步骤
用三个向量B,T,s,存放多项式0，T（x），0,3，2,+1W的系数。
（1）构造Qq（x），设°o（x）=bo,根据④式，得bo=1。
再根据⑦®⑤式，计算:
do=n+1
n
Sy,
~d7
k=O
a°~~d7
最后将^020W项展开后累加到拟合多项式中，则q°bo=a。
（2）构造Q（x）,设Q（Q=/o+/\x,根据递推公式④，则可知，仏=一Q。
，t=
io由公式⑦③⑤⑥求得，
d产丈03）
20
tykQSxk）
k=0
最后将展开累加到拟合多项式①中，有
Qo+qi。
*a。
⑶对于j=2,3,逐步递推Q（x）
根据递推公式④有
0（X）=（X-%T）0鬥（X）-0鬥Q,-2（X）
=（X■図J（匚tx~'+…•+Ax+?
o>0j_AbH£2+•+bm+bo）
Qj（x）=SjX+Sj-lXJ+……+51X+So
则可以得到计算5；（k=O,l,…J）的公式：
S产-OCj-W卩厂\bik=j-2,….,2,1
然后分别根据⑦式③式⑤式与⑥式计算卜•列屋:
d广丈05）
q厂土几Q'xjdj
“=工尢3）0（竝）a
再将q,QSx）项展开后累加到拟合多项式①中，即有
ak+QjSk=>QjJ,…•丄°
QjSj^aj
最后，，为了便于循坏使用向量BT,与S,应将向量T传递给E,向量S传递送给T,即tk=>bk»k=j-l,…,1,0
mk=j,…,1,0
在实际计算中，为了防止运算溢出，将兀，用Xk~X来代替，其中
Jt=o
在这种情况下，拟合多项式的形式为
_2/w
pmM=a0+al（x-x）+a2（x—兀）+…+am（x—兀）
2.3程序流程图
灼卜二＞
c=t/d2；p=g/d2；q=d2/dlL；
dl=d2；3[j]=c-s[j];1[J]=4jl；
d2+=q*qx-»-=v[i]*q;
2[k]4=c"（町；t>[k]=t[k];t[k]=s[k];
K=j-l;k>=O;k-
dt（O）=O.O;dt[l]=a.O;dt[2]=0.0
q=q*（x[i]-zl*s[k]
i=0;iq=a[ml]
K=m-2;k>=0；k-
q=ak]+q*（xil-z）;
p=q-vi'i;
3实验结果分析
温度数据:
51.1
52.8
54.8
57.258.3
62.7
65.2
67.7
70.673.575.77&6
80.8
84.8
87.8
89.592.1
96.4
103.1
112.5
120.8134.7
152.4159.1
氧气数据：
20.13
19.9
19.61
18.9819.61
19.32
19.73
19.5919.38
20.1819.9819.58
19.74
19.26
；19.5918.77
18.66
17.47
17.01
16.3315.95
13.765.915.43
氮气数据:
79.579.6779.5980.248080.3179.3580.180.3779.779.7880.15
8080.4680.0980.980.9982.0682.3383.0283.3184.8490.8593.96
实验结果:
1-
■正交多项式拟合
输出温度数据：
51.100000
52・800000
54.800000
5?
・200000
58・300000
62.700000
65.200000
67.700000
70.600000
73.500000
75.700000
78.600000
80.800000
84.800000
87.800000
89.500000
92.100000
96.400000
103.100000
112.500000
120.800000
输出氧气数据：
134.700000
152.400000
159.100000
20.130000
仅900000
19.610000
18・980000
19.610000
19.320000
19.730000
19.590000
19.380000
20.180000
19.980000
19.580000
19.740000
19.260000
19.590000
18.770000
18.660000
17.470000
17.010000
16.330000
15.950060
a<0>=19.060178a=-0.050815a<2>=-0.001271
13・760000
5.910000
5.430000
3>=-0.000010
dt=8.06955?
dt<1>=10.859434dt<2>=1.385219
输出氧气与温度拟合的正交多项式:
输岀温度数据：
51.100000
52.800000
54.800000
57.200000
58.300000
62.700000
65.200000
6?
・700000
70.600000
73.500000
75.700000
78.600000
80.800000
84.S00000
87.S00000
89.500000
92.100000
96・400000
103.100000
112.500000
120.800000
134・700000
152.400000
159.100000
输出氮气数据，
79.500000
79・670000
79.590000
80.240000
80・000000
80.310000
79.350000
80.130000
80.370000
79.700000
79.780000
80.150000
80・000000
80.460000
80・090000
80.900000
80.990000
82.060000
82.330000
83.020000
83.310008
a<0>=80.628611a<1）=0.038346
a<2>=0・000776
a<3>=0・0000丄6
84.840008
90.850008
93.960008
dtC0>=5.624939
dt<1>=9.117542dt<2>=1
.167073
输出氮气与温度拟合的正交多项式:
y=80.628611+C0.038346>*^l+<0.000?
76>*^2+*^3
Pressanykeytocontinue
注：
dt（0）为误差的平方和，dt⑴为误差的绝对值之和，dt
（2）误差绝对值最人值。
氧气与温度的拟合曲线：
20
*拟合点
——三次拟合-
18
16
14
12
10
8
6
440
60
80
氮气与温度的拟合曲线:
1
100
120
140
160
实验结论
这次实验我们的拟合曲线是选择三次拟合曲线，虽说更高次的拟合曲线可以达到更好的效果，但是由于在计算机计算的过程中舍入误差的存在，使得次数高的项系数为零。
由于数据是观测得来，而我们的误差最人不超过1.4,所以误差在允许范閑内，故从误差以及各方面的考虑，我们最终选择三次来拟合曲线。
为了能够达到视觉上的效果，我们在实验结果处附上了用matlab所拟合得到的曲线，从图中可以看出，随着温度的不断增加氧气的含屋在降低，且降低率随温度的增加而增加,但是对于氮气却是随着温度的增加而增加，且增加率随温度的增加而增加。
由于温度的不断增加，经过一定的化学变换，放出氮气，同时消耗氧气，而且在温度相对高时（在一定的温度范围内），其化学反应的速率快，这就是对上述结果的一个解释。
5问题归纳与总结
在本次试验中，组员在分析问题的过程中主要遇到了一下几方面的问题，下面一一表述并给出解决办法。
问题一：
0,分别Q,（x）之间的关系
解决办法:
观察公式，找出关系，将公式分解，分步求出分子分母，将分母用d广丈0（竝）
表示。
问题2：
Q（X）的迭代问题
解决办法：
分别用三个向量E,T,S分别存在0/t⑴，QSX）,0鬥⑴的系数，用
tk=>brsk=>tr依次迭代
问题3：
系数拟合问题
解决办法：
将上一次正交多项式次数相等的对应系数加到卞一次的系数，
即：
久+乞$=>久，k=j-l,….,1,0
问题4：
Q（x）的算法
解决方法：
Q（Xk）是一个j次多项式，从j次到1次，递减乘x来增加次数，再加上
s[k].
Qj（兀）=（（$/+$」+$』
Qj（兀）=（（工X+5；-1X+Sj-2）x+57-3）
Qj（x）=SjX^Sj-,X'+……+S1X+So
问题5：
溢出问题
解决方法：
由于给定温度数据过犬，次数大时会出现溢出，所以采用平移思想用X-X
来表示X,图像不变。
参考文献
李庆扬，王能超，易大义.数值分析.北京：
清华大学岀版社.200S
附录（源代码）
#include"stdio.h"
#include”math.h"
doublespu（doublex[].doublev[].mtn.doublea[],iiitm.doubledt[]）
intij,k;
doublez,p,c,g,q,d1,d2,s[20],t[20],b[20];for（i=0;i<=m-l;i++）//系数的初始化a[i]=0.0;
m=n;
if20）
m=20;
z=0.0;
for（i=0;i<=n-l;i++）
z=z+x[i]/（1.0*n）;//温度平均值
b[0]=1.0y/Q0（x）多项式系数dl=1.0*n;
P=0.0;
c=0.0;
for（i=0;i<=n-l;i++）
p=p+（x[i]-z）;c=c+v[i];
}
c=c/dl;
p=p/dl;
a[0]=c*b[0];
〃构造Ql（x）
t[l]=l・O；
t[O]=-p；
d2=0.0;
c=0.0;
g=0.0;
for（i=0;i<=n-l;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/dl;
dl=d2;
a[l]=c*t[l];a[O]=c*t[O]+a[O];
}for（j=2j<=m-l;j++）〃构造Qj（x）
{
s[j]=t|j-l]；s[j-l]=-p*t[j-l]4-t[j-2]；珂>=3）
fbr（k=j-2;k>=l;k-）
s[k]=-p*t[k]+t[k-l]-q*b[k];
s[O]=-p*t[O]-q*b[O];d2=0.0;c=0.0;g=0.0;fbr（i=O;i<=n-l;i++）
{
q=s[j];
fbr（k=j-l;k>=O;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;
dl=d2;
a[j]=c*s[j];
fbr（k=j-l;k>=O;k-）
{
a[k]=c*s[k]+a[k];
b[k]=t[k];
t[k]=s[k];
}
〃计算误差的平方和、误差的绝对值之和与误差绝对值的最人值
dt[O]=O.O;
dt[l]=0.0;dt[2]=0.0;
fbr（i=O;i<=n-l;i++）
{
q=a[m-l];
for（k=m-2;k>=0;k-）
q=a[k]+q*（x[i]-z）;
p=q-yW;
if（fabs（p）>dt[2]）
dt[2]=fabs（p）;
dt[O]=dt[O]+p*p;
dt[l]=dt[l]+fabs（p）;
}
returnz;
}
voidmain（）
{
inti;
doublea[8],dt[3],z;
doublex[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・&87・8,89・乂92丄96・4」03丄112・5」20・&134・7,152・4,159・1};
doubleyl[24]={20.13J9.9J9.61J8.98,19.614932,19.73,19.59,1938,20.18,19.98//氧气数据
I9・58,19・74,19・26」9・59,18・77,18・66」7・47,17・0I」6・33」5・95,13・76,5・91,5・43};
doubley2[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};//氮/[数据
pnntf（"正交多项式拟合\n”）；
pmitf（”输出温度数据:
\n”）；
fbr（i=0;i<=23;i++）
pimtf（"\n输出氧气数据:
\n”）；
fbr（i=0;i<=23;i++）
z=spir（x,vl,24,a,4,dt）;〃函数的调用来求拟合函数系数、误差
pnntf（H\nM）;
for（i=0;i<=3;i-H-）
pnntf（H\nM）;
for（i=0;i<=2;i-H-）
prmtf（ndt（%2d）=%lff\i,dt[i]）;//误差的输出
pmitf（”输出氧气与温度拟合的正交多项式:
\n\n”）；pnntf（Mp（x）=%lf\a[O]）;
for（i=l;i<=3;i-H-）
pnntf（”+（%lf）*（x・％lf）A%d”Q[i],z，i）;
pnntf（H\n\QM）;
printf（”输出温度数据:
\n“）;
for（i=0;i<=23;i++）
pimtf（"\n输出氮气数据:
\n”）；
for（i=0;i<=23;i++）
pnntff%lf\L，y2[i]）;
z=spir（x,v2,24,a,4,dt）;//函数的调用来求拟合函数系数、误差pnntf（H\nM）;
for（i=0;i<=3;i-H-）pimtf（"a（%2d）=%lfn",i,a[i]）y/函数系数的输出pnntf（H\nM）;
for（i=0;i<=2;i-H-）
prmtf（"dt（%2d）=%lf”,i,dt[i]）;〃误差的输出
pnntf（”输出氮气与温度拟合的正交多项式:
\n\n"）;piintf（Mv（x）=%lf\a[O]）;
for（i=l;i<=3;i-H-）
pnntf（”+（%lf）*（x・％lf）A%d”Q[i],z，i）;
pnntf（H\nM）;
升级会员 