1、 4 曲线拟合的最小二乘法曲线拟合的最小二乘法 1 最小二乘法及其计算最小二乘法及其计算 在函数的最佳平方逼近中 如果 只在一组离散点集 上给定，这就是科学实验中经常见到的实验数据 的曲线拟合.,)(baCxf)(xf,1,0,mixi,1,0),(miyxii1 记误差,1,0,)(*miyxSiii则 的各分量分别为 个数据点上的误差.T10),(mm 问题为利用 求出一个函数,1,0),(mixfyii)(*xSy 与所给数据 拟合.,1,0),(miyxii2 设 是 上线性无关函数族，)(,),(),(10 xxxn,baC在 中找一函数 ，)(,),(),(10 xxxspann)
2、(*xS使误差平方和miiimiiyxS02*02)(,)(min02)(miiixSyxS这里)()()()(1100 xaxaxaxSnn).(mn 3 这个问题称为最小二乘逼近，几何上称为曲线拟合的最小二乘法.用最小二乘求拟合曲线时,首先要确定 的形式.)(xS 确定 的形式问题不仅是数学问题,还与问题的实际背景有关.)(xS 通常要用问题的运动规律及给定的数据进行数据描图,确定 的形式,然后通过实际计算选出较好的结果.)(xS4 为了使问题的提法更有一般性，通常在最小二乘法中考虑加权平方和)()()()(1100 xaxaxaxSnn)(mn,)()(),(02miiiiyxSx这里
3、是 上的权函数，它表示不同点 处的数据比重不同.0)(x,ba)(,(iixfx就是 次多项式.)(xSn 若 是 次多项式，)(xkk 的一般表达式为线性形式.)(xS5 这样，最小二乘问题就转化为求多元函数),(10naaaIminjiijjixfxax002)()()(的极小点 问题.),(*1*0naaa 用最小二乘法求拟合曲线的问题，就是在 中求一函数 ，)(xS)(*xSy 由求多元函数极值的必要条件，有 minjikiijjikxxfxaxaI000)()()()(2).,1,0(nk使误差取得最小.)()()()(1100 xaxaxaxSnn)(mn 6若记,)()()(),
4、(0miikijikjxxxkmiikiikdxxfxf0)()()(),().,1,0(nk上式可改写为 knjjjkda0),().,1,0(nk这个方程称为法方程法方程，可写成矩阵形式7其中,),(,),(T10T10nndddaaada.),(),(),(),(),(),(),(),(),(101110101000nnnnnnG,dGa 要使法方程有唯一解，就要求矩阵 非奇异，G而 在 上线性无关不能推出)(,),(),(10 xxxn,ba矩阵 非奇异，必须加上另外的条件.Gknjjjkda0),().,1,0(nk8 显然 在任意 个点上满足哈尔条件.nxx,1)(nmm哈尔条件，
5、则法方程 的系数矩阵 非奇异，如果 在 上满足,)(,),(),(10baxxxnmix0函数 的最小二乘解为)(xf定义定义1010设 的任意线,)(,),(),(10baxxxn性组合在点集 上至多只有 个)(,1,0,nmmixin不同的零点，则称 在点集 )(,),(),(10 xxxn,1,0,mixi上满足哈尔哈尔(Haar)条件条件.,1,0,*nkaakk方程存在唯一的解从而得到于是knjjjkda0),().,1,0(nk9,)()()()()()(0202*miiiimiiiixfxSxxfxSx这样得到的 ，)(*xS对任何的 )(xS).()()()(*1*10*0*x
6、axaxaxSnn都有故 确是所求最小二乘解.)(*xS)()()()(1100 xaxaxaxSnn)(mn 10一般可取 ，但这样做当 时，,1nxxspan3n通常对 的简单情形都可通过求法方程得到 1n).(*xS 给定 的离散数据 ，,1,0),(miyxii)(xf求解法方程时将出现系数矩阵 为病态的问题，我们在下面考虑用正交多项式的方法解决。用正交多项式的方法解决。Gknjjjkda0),().,1,0(nk11 例例1 1113125.8865.4454321iiifx已知一组实验数据如下，求它的拟合曲线.12 解解 从图中看到各点在一条直线附近，故可选择线性函数作拟合曲线，将
7、所给数据在坐标纸上标出，见图3-4.图3-413 令,)(101xaaxS,8),(4000ii,22),(),(400110iiix,74),(40211iiix,47),(400iiiff.5.145),(401iiiifxf,1)(,1,40 xnm这里故,)(1xx14.5.1457422,472281010aaaa解得.20.1,58.210aa.20.158.2)(*1xxS可得方程组 于是所求拟合曲线为knjjjkda0),().,1,0(nk15 关于多项式拟合，Matlab中有现成的程序 ),(polyfitmyxa 其中输入参数 为要拟合的数据，为拟合多项式的次数，yx,m
8、输出参数 为拟合多项式的系数.a 利用下面的程序，可在Matlab中完成上例的多项式拟合.16x=1 1 2 3 3 3 4 5;f=4 4 4.5 6 6 6 8 8.5;aa=poly(x,f,1);y=polyval(aa,x);plot(x,f,r+,x,y,k)xlabel(x);ylabel(y);gtext(y=s1(x)）17结果如下：18有时根据给定数据图形，其拟合函数 表面上不是线性模型的形式，但通过变换仍可化为线性模型.,ln)(lnbxaxS 例如，bxaxSe)(若两边取对数得此时，若令,ln),(ln)(bBaAxSxS,)(BxAxS这样就变成了线性模型.19 例
9、例2 2设数据 由表3-1给出，)4,3,2,1,0)(,(iyxii,ebxay 用最小二乘法确定 及 .ab46.845.753.679.510.500.275.150.125.100.143210iiyxi1表3 解解,lniiyy表中第4行为通过描点可以看出数学模型为它不是线性形式.,ebxay 用给定数据描图可确定拟合曲线方程为两边取对数得20 若令,ln,lnaAyy先将 转化为),(iiyx),(iiyx为确定 ，bA,根据最小二乘法，取,1)(,)(,1)(10 xxxx.lnlnbxay.,1,xbxAy则得数据表见表3-1.得,5),(00,5.7),(4010iix,87
10、5.11),(40211iix135.2008.2876.1756.1629.146.845.753.679.510.500.275.150.125.100.143210iiiyyxi1表3 21.422.14),(401iiiyxy,404.9),(400iiyy故有法方程.422.14875.1150.7,404.950.75bAbA解得.071.3e,505.0,122.1AabA.e071.3505.0 xy 于是得最小二乘拟合曲线为 22 利用下面的程序，可在Matlab中完成曲线拟合.x=1.00 1.25 1.50 1.75 2.00;y=5.10 5.79 6.53 7.45
11、8.46;y1=log(y);aa=poly(x,y1,1);a=aa(1);b=exp(aa(2);y2=b*exp(a*x);plot(x,y,r+,x,y2,k)xlabel(x);ylabel(y);gtext(y=a*exp(bx);23结果如下：24 2 用正交多项式做最小二乘拟合用正交多项式做最小二乘拟合 如果 是关于点集)(,),(),(10 xxxn,0,0)()()(),(0kmiikijikjAxxx,kj,kj 用最小二乘法得到的法方程组，其系数矩阵 是病态的.G带权 正交的),1,0(mixi),1,0()(mixi函数族，即knjjjkda0),().,1,0(nk
12、（5.6）25miikimiikiikkkkxxxxfxfa020*)()()()()(),(),().,1,0(nk则方程的解为 且平方误差为.)(),(),(02*nkkkaAffknjjjkda0),().,1,0(nk26 接下来根据给定节点 及权函数 mxxx,10,0)(x构造带权 正交的多项式 .)(x)(xPn 注意 ，用递推公式表示 ，即mn)(xPk)()()()(),()()(,1)(1110110 xPxPxxPxPxxPxPkkkkk).1,2,1(nk这里 是首项系数为1的 次多项式，)(xPkk根据 的)(xPk正交性，得27),(),(11kkkkPPPP).1
13、,2,1(nkmiikimiikikxPxxPx02102)()()()(miikimiikiikxPxxPxx02021)()()()(下面用归纳法证明这样给出的 是正交的.)(xPk)(),()(),(xPxPxPxxPkkkk),(),(kkkkPPPxP28),(),(),(0010010PPxPPPP 假定 对 及)(0),(slPPsl1,1,0ls,1,0kl要证 对 均成立.0),(1skPPks,1,0有),(),)(),(111skkskkskPPPPxPP 由 的表达式，有 1),(),(),(),(00000000PPPPPxPxPP.0nk 均成立，).,(),(),
14、(11skkskkskPPPPPxP)()()()(),()()(,1)(1110110 xPxPxxPxPxxPxPkkkkk).1,2,1(nk29 而 ，11ks,0),(),(skskxPPPxP 当 时，2 ks.0),(1skPP 另外，是首项系数为1的 次多项式，它可由)(xxPs1s由归纳法假定，当 时20ks,0),(slPP.0),(1skPP110,sPPP的线性组合表示.由归纳法假定又有),(),)(),(111skkskkskPPPPxPP).,(),(),(11skkskkskPPPPPxP30由假定有),(),(11kkkkxPPPxP 再考虑),(),(),()
15、,(1111111kkkkkkkkkkPPPPPxPPP),(10kjjjkkPcPP).,(kkPP利用 表达式及以上结果，得 k),(),(),(11111kkkkkkkPPPxPPP.0),(),(kkkkPPPP),(),(11kkkkkPPPP31),(),(),(),(111kkkkkkkkkkPPPPPxPPP至此已证明了此多项式 )(xPk组成一个关于点集 的正交系.ix 用正交多项式 的线性组合作最小二乘曲线拟合，)(xPk只要根据公逐步求 的同时，)(xPk相应计算出系数*ka.0),(),(),(),(kkkkkkkkPPPPPxPPxP最后，由 和 的表达式（5.11）
16、有 kk),(),(11kkkkkPPPP),(),(1kkkkkPPPxP32),(),(*kkkkPPPfa),1,0(nk并逐步把 累加到 中去，最后就可得到所求的)(*xPakk)(xS).()()()(*1*10*0 xPaxPaxPaxSynn 用这种方法编程序不用解方程组，只用递推公式，并且当逼近次数增加一次时，只要把程序中循环数加1，其余不用改变.这里 可事先给定或在计算过程中根据误差确定.n)()()()()(200ikmiimiikiixPxxPxfx拟合曲线33以上述例1为例 先求正交系 20.15.1325.16),(),(88.5847),(),(75.2;75.28
17、22),(),(;1111*1000*0100010PPPfaPPPfaxPPPPxxPP令.20.158.2)2.75-20.188.5)(*1xxxS（34例3设X=1.00,1.25,1.50,1.75,2.00,在X上定义内积 5 (f,g)=xi f(xi)g(xi)i=1 1)在函数系 1,x2中求一个X上的正交函数系.2)用最小二乘法求一个形如y=a+bx2的经验公式，使它与下列数据拟合.xi 1.00 1.25 1.50 1.75 2.00yi 3.00 4.50 5.50 7.00 9.0035625.21X625.2;625.20)(0),(;112215151325110
18、210 xxPxxxxPPxPPiiiiiii，上正交系为）令36.93.121.1)2.625-93.128.6)(93.1),(),(28.6),(),()625.2()(S)222*1111*1000*02*1*0*xxxSPPPfaPPPfaxaax（设375 最佳一致逼近多项式最佳一致逼近多项式 1.基本概念及其理论基本概念及其理论 设,baCf 在 中求多项式,1nnxxspanH),(*xPn)()(max*xPxfnbxa这就是最佳一致逼近或切比雪夫逼近问题.使其误差.|)()(|maxmin,xpxfnbaxHPnn38使误差 对连续函数 ，它不能用有限个线性无关的,)(ba
19、Cxf函数表示，故 是无限维的，但它的任一元素 ,baC)(xf均可用有限维的 逼近，nHxp)()()(maxxpxfbxa(为任给的小正数)，这就是著名的魏尔斯特拉斯定理.39使 定理定理1 1总存在一设 ，,)(baCxf则对任何 ，0个代数多项式 ，)(xp|)()(|max,xpxfbax在 上一致成立.,ba 伯恩斯坦1912年给出的证明是一种构造性证明.他根据函数整体逼近的特性构造出伯恩斯坦多项式,)()(),(0nkknxPnkfxfB40 为二项式展开系数，并证明了!)1()1(kknnnkn)(),(limxfxfBnn在 上一致成立；1,0 若 在 上 阶导数连续，则)(
20、xf1,0m其中,)1()(knkkxxknxP).(),(lim)()(xfxfBmmnn 这个结果不但证明了定理1，而且给出了 的一个逼近多项式.)(xfnkknxPnkfxfB0)()(),(41定理定理2(最佳一致逼近的存在性）最佳一致逼近的存在性）设设f(x)在在a,b上连续，则存在上连续，则存在pn*(x)Hn使使.|)()(|maxmin,xpxfnbaxHPnn)()(max*xPxfnbxa下面研究求下面研究求pn*(x)的方法定义：设定义：设f(x)Ca,b,p(x)Hn,若若x=x0时时则称则称x0为为p(x)的偏差点的偏差点.|)()(|max|)()(|,00 xfx
21、pxfxpbax.,)()(,)()(000000为负偏差点若为正偏差点，若xxfxpxxfxp42要证明的是其中,)1()()(kkkxfxP这样的点组称为切比雪夫交错点组切比雪夫交错点组.证明证明 假定在 上有 个点使上式成立，,ba2n 定理定理3 3即有 个点 ，bxxxan 2212n“负”的偏差点，在 上至少有 个轮流为“正”、)(xP,ba2n是 的最佳逼近多项式nHxP)(,baCf 的充分必要条件是使是 在 上的最佳逼近多项式.)(xf,ba)(xP只证充分性.|)()(|max,xfxpbax43 用反证法，若存在 ，)()(,)(xPxQHxQn|,)()(|max|)(
22、)(|max,xPxfxQxfbaxbax由于)()()()()()(xfxQxfxPxQxP故 也在 个点上轮流取“+”、“-”号.)()(xQxP2n使 由连续函数性质，它在 内有 个零点，但因,ba1n0)()(xQxP是不超过 次的多项式，n不能超过 .n所以它的零点个数在点 上的符号与221,nxxx)2,1)()(nkxfxPkk一致，44 这说明假设不对，故 就是所求最佳逼近多项式.)(xP 必要性证明略.推论推论1 1 若 ，,baCf 充分性得证.则在 中存在唯一的最佳逼近nH多项式.45 零偏差最小问题.1|0|maxmin|0|maxmin|maxmin|max|)(|m
23、axmin|)(|max)(,)()(1,1,1,1,)(1-n011p*111-n0111-n0*11111*1111-n011-n0*1*111*1*1111次多项式的为首项系数为其中既设的最佳一致逼近。为使求设nxaxpppxaxxaxxpxxpxxaxpxaxpxxfpHHpxxxfjjjnnnxHnxHpjjjnxHpjjjnxnnxHpnnxjjjnjjjnnnnnnnnnnnnnn46 证明证明)(21)(1xTxnnn)(max21)(max11111xTxnxnnx且点 是 的切比雪夫交错点组，),1,0(cosnknkxk)(Txn 定理定理4 4在区间 上所有最高次项系数
24、为1的 次多1,1n项式中，)(T21)(1xxnnn与零的偏差最小，.211n其偏差为由于),(*1xPxnn,211n47由定理3可知，即 是与零的偏差最小的多项式.)(xn区间 上 在 中最佳逼近多项式1,1nx1nH),(*1xPn为定理得证.48.min)()(max*211xPxfx由定理6可知，)(21)()(),(212)()(3*232*2xTxPxfxTxPxf既多项式 与零偏差最小，)()(*2xPxf 解解由题意，所求最佳逼近多项式 应满足)(*2xP当时，故例例3 3求 在 上的最佳2次逼122)(23xxxxf1,1近多项式.49)(21)()(3*2xTxfxP就
25、是 在 上的最佳2次逼近多项式.)(xf1,11272xx50 2 最佳一次逼近多项式最佳一次逼近多项式 定理3给出了 的特性，这里讨论具体求法.)(xP 先讨论 的情形.1n 假定,2baCf 且 在 内不变号，)(xf ),(ba 根据定理3可知,至少有3个点,321bxxxa.)()()()(max)1()()(111的最大（小）值点）也就是求xfxPxQxfxPxfxPbxakkk).3,2,1,1(k求最佳一次逼近多项式 .xaaxP101)(我们要51,0)()()(2122xfaxfxP即 .12)(axf 由于 在 上不变号，)(xf ,ba故 单调，)(xf)(1xfaQ在
26、内只有一个零点，记为 ，),(ba2x 另外两个偏差点必是区间端点，即 且,21bxax)()()()(11bfbPafaP由此得到 于是满足).()(221xfxP52解出 代入得).()()();()(2102101010 xaaxfafaaabfbaaafaaa),()()(21xfabafbfa这就得到最佳一次逼近多项式 ，其几何意义如图3-3.)(1xP.2)()(2)()(220 xaabafbfxfafa53 直线 与弦MN平行，且通过MQ的中点D，)(1xPy).2()()(21212xaxaxfafy图3-3其方程为54由 可算出 例例1 1求 在 上的最佳一次逼近多项式.2
27、1)(xxf1,0 解解,414.0121a又 ,1)(2xxxf,4551.02122x由 得,121222 xx故解得.0986.11)(222xxf),()()(21xfabafbfa.2)()(2)()(220 xaabafbfxfafa55,414.0955.0)(1xxP即;10,414.0955.012xxx误差限为.045.0)(1max1210 xPxx于是得 的最佳一次逼近多项式为 21x,955.0221121220 xaxa56在上式中若令,1abx.414.0955.022baba则.414.0955.012abab从而可得一个求根式的公式;10,414.0955.012xxx57
