1、3.5 曲线拟合的最小二乘法曲线拟合的最小二乘法 3.5.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 xxxs
2、pann)(*xS使误差平方和miiimiiyxS02*0222)(,)(min02)(miiixSyxS(5.1)这里)()()()(1100 xaxaxaxSnn).(mn(5.2)3 这个问题称为最小二乘逼近,几何上称为曲线拟合的最小二乘法.用最小二乘求拟合曲线时,首先要确定 的形式.)(xS 确定 的形式问题不仅是数学问题,还与问题的实际背景有关.)(xS 通常要用问题的运动规律及给定的数据进行数据描图,确定 的形式,然后通过实际计算选出较好的结果.)(xS4 为了使问题的提法更有一般性,通常在最小二乘法中考虑加权平方和)()()()(1100 xaxaxaxSnn)(mn(5.2),
3、)()(0222miiiiyxSx(5.3)这里 是 上的权函数,它表示不同点 处的数据比重不同.0)(x,ba)(,(iixfx就是 次多项式.)(xSn 若 是 次多项式,)(xkk 的一般表达式为(5.2)表示的线性形式.)(xS5 这样,最小二乘问题就转化为求多元函数),(10naaaIminjiijjixfxax002)()()((5.4)的极小点 问题.),(*1*0naaa 用最小二乘法求拟合曲线的问题,就是在形如(5.2)的 中求一函数 ,)(xS)(*xSy 由求多元函数极值的必要条件,有 minjikiijjikxxfxaxaI000)()()()(2).,1,0(nk使(
4、5.3)取得最小.)()()()(1100 xaxaxaxSnn)(mn(5.2).)()(0222miiiiyxSx(5.3)6若记,)()()(),(0miikijikjxxx(5.5)kmiikiikdxxfxf0)()()(),().,1,0(nk上式可改写为 knjjjkda0),().,1,0(nk(5.6)这个方程称为法方程法方程,可写成矩阵形式7其中,),(,),(T10T10nndddaaada.),(),(),(),(),(),(),(),(),(101110101000nnnnnnG(5.7),dGa 要使法方程(5.6)有惟一解,就要求矩阵 非奇异,G而 在 上线性无关
5、不能推出)(,),(),(10 xxxn,ba矩阵 非奇异,必须加上另外的条件.Gknjjjkda0),().,1,0(nk(5.6)8 显然 在任意 个点上满足哈尔条件.nxx,1)(nmm哈尔条件,则法方程(5.6)的系数矩阵(5.7)非奇异,如果 在 上满足,)(,),(),(10baxxxnmix0函数 的最小二乘解为)(xf 定义定义1010设 的任意线,)(,),(),(10baxxxn性组合在点集 上至多只有 个)(,1,0,nmmixin不同的零点,则称 在点集 )(,),(),(10 xxxn,1,0,mixi上满足哈尔哈尔(Haar)条件条件.,1,0,*nkaakk方程(
6、5.6)存在惟一的解从而得到于是knjjjkda0),().,1,0(nk(5.6)9,)()()()()()(0202*miiiimiiiixfxSxxfxSx这样得到的 ,)(*xS对任何形如(5.2)的 ,)(xS).()()()(*1*10*0*xaxaxaxSnn都有故 确是所求最小二乘解.)(*xS)()()()(1100 xaxaxaxSnn)(mn(5.2)10一般可取 ,但这样做当 时,,1nxxspan3n通常对 的简单情形都可通过求法方程(5.6)得到 1n).(*xS 给定 的离散数据 ,,1,0),(miyxii)(xf 例如,bxaxSe)(,ln)(lnbxaxS
7、求解法方程(5.6)将出现系数矩阵 为病态的问题,G 有时根据给定数据图形,其拟合函数 表面上)(xfy 不是(5.2)的形式,但通过变换仍可化为线性模型.若两边取对数得knjjjkda0),().,1,0(nk(5.6))()()()(1100 xaxaxaxSnn)(mn(5.2)11 例例7 7113125.8865.4454321iiifx这样就变成了形如(5.2)的线性模型.此时,若令 则,ln),(ln)(bBaAxSxS,)(BxAxS已知一组实验数据如下,求它的拟合曲线.12 解解 从图中看到各点在一条直线附近,故可选择线性函数作拟合曲线,将所给数据在坐标纸上标出,见图3-4.
8、图3-413 令,)(101xaaxS,8),(4000ii,22),(),(400110iiix,74),(40211iiix,47),(400iiiff.5.145),(401iiiifxf,1)(,1,40 xnm这里故,)(1xx14.5.1457422,472281010aaaa解得.13.1,77.210aa.13.177.2)(*1xxS由(5.6)得方程组 于是所求拟合曲线为knjjjkda0),().,1,0(nk(5.6)15 关于多项式拟合,Matlab中有现成的程序 ),(polyfitmyxa 其中输入参数 为要拟合的数据,为拟合多项式的次数,yx,m输出参数 为拟合
9、多项式的系数.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 例例8 8设数据 由表3-1给出,)4,3,2,1,0)(,(iyxii,ebxay 用最小二乘法确定 及 .ab46.845.753.679.510.500.275.150.125.100.143210iiyxi1表3 解解,lniiyy表中第4行为通过
10、描点可以看出数学模型为它不是线性形式.,ebxay 用给定数据描图可确定拟合曲线方程为两边取对数得19 若令,ln,lnaAyy先将 转化为),(iiyx),(iiyx为确定 ,bA,根据最小二乘法,取,1)(,)(,1)(10 xxxx.lnlnbxay.,1,xbxAy则得数据表见表3-1.得,5),(00,5.7),(4010iix,875.11),(40211iix135.2008.2876.1756.1629.146.845.753.679.510.500.275.150.125.100.143210iiiyyxi1表3 20.422.14),(401iiiyxy,404.9),(4
11、00iiyy故有法方程.422.14875.1150.7,404.950.75bAbA解得.071.3e,505.0,122.1AabA.e071.3505.0 xy 于是得最小二乘拟合曲线为 21 利用下面的程序,可在Matlab中完成曲线拟合.x=1.00 1.25 1.50 1.75 2.00;y=5.10 5.79 6.53 7.45 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);22结果如
12、下:23 3.5.2 用正交多项式做最小二乘拟合用正交多项式做最小二乘拟合 如果 是关于点集)(,),(),(10 xxxn,0,0)()()(),(0kmiikijikjAxxx,kj,kj(5.8)用最小二乘法得到的法方程组(5.6),其系数矩阵 是病态的.G带权 正交的),1,0(mixi),1,0()(mixi函数族,即knjjjkda0),().,1,0(nk(5.6)24miikimiikiikkkkxxxxfxfa020*)()()()()(),(),().,1,0(nk(5.9)则方程(5.6)的解为 且平方误差为.)(02*2222nkkkaAfknjjjkda0),().,
13、1,0(nk(5.6)25 接下来根据给定节点 及权函数 mxxx,10,0)(x构造带权 正交的多项式 .)(x)(xPn 注意 ,用递推公式表示 ,即mn)(xPk)()()()(),()()(,1)(1110110 xPxPxxPxPxxPxPkkkkk).1,2,1(nk(5.10)这里 是首项系数为1的 次多项式,)(xPkk根据 的)(xPk正交性,得26),(),(11kkkkPPPP(5.11)).1,2,1(nkmiikimiikikxPxxPx02102)()()()(miikimiikiikxPxxPxx02021)()()()(下面用归纳法证明这样给出的 是正交的.)(
14、xPk)(),()(),(xPxPxPxxPkkkk),(),(kkkkPPPxP27),(),(),(0010010PPxPPPP 假定 对 及)(0),(slPPsl1,1,0ls,1,0kl要证 对 均成立.0),(1skPPks,1,0由(5.10)有),(),)(),(111skkskkskPPPPxPP 由(5.10)第二式及(5.11)中 的表达式,有 1),(),(),(),(00000000PPPPPxPxPP.0nk 均成立,(5.12)).,(),(),(11skkskkskPPPPPxP)()()()(),()()(,1)(1110110 xPxPxxPxPxxPxPk
15、kkkk).1,2,1(nk(5.10))()()()(),()()(,1)(1110110 xPxPxxPxPxxPxPkkkkk).1,2,1(nk(5.10)28 而 ,11ks,0),(),(skskxPPPxP于是由(5.12),当 时,2 ks.0),(1skPP 另外,是首项系数为1的 次多项式,它可由)(xxPs1s由归纳法假定,当 时20ks,0),(slPP.0),(1skPP110,sPPP的线性组合表示.由归纳法假定又有),(),)(),(111skkskkskPPPPxPP(5.12)).,(),(),(11skkskkskPPPPPxP29由假定有),(),(11k
16、kkkxPPPxP 再考虑(5.13)),(),(),(),(1111111kkkkkkkkkkPPPPPxPPP),(10kjjjkkPcPP).,(kkPP利用(5.11)中 表达式及以上结果,得 k),(),(),(11111kkkkkkkPPPxPPP.0),(),(kkkkPPPP),(),(11kkkkkPPPP30),(),(),(),(111kkkkkkkkkkPPPPPxPPP至此已证明了由(5.10)及(5.11)确定的多项式 )(xPk组成一个关于点集 的正交系.ix 用正交多项式 的线性组合作最小二乘曲线拟合,)(xPk只要根据公式(5.10)及(5.11)逐步求 的同
17、时,)(xPk相应计算出系数*ka.0),(),(),(),(kkkkkkkkPPPPPxPPxP最后,由 和 的表达式(5.11)有 kk),(),(11kkkkkPPPP),(),(1kkkkkPPPxP31),(),(*kkkkPPPfa),1,0(nk并逐步把 累加到 中去,最后就可得到所求的)(*xPakk)(xS).()()()(*1*10*0 xPaxPaxPaxSynn 用这种方法编程序不用解方程组,只用递推公式,并且当逼近次数增加一次时,只要把程序中循环数加1,其余不用改变.这里 可事先给定或在计算过程中根据误差确定.n)()()()()(200ikmiimiikiixPxx
18、Pxfx拟合曲线323.6 最佳平方三角逼近与快速傅里叶变换最佳平方三角逼近与快速傅里叶变换 当 是周期函数时,显然用三角多项式逼近 比用代数多项式更合适,本节主要讨论用三角多项式做最小平方逼近及快速傅里叶变换,快速傅里叶变换,简称FFT算法.)(xf)(xf33 3.6.1 最佳平方三角逼近与三角插值最佳平方三角逼近与三角插值 设 是以 为周期的平方可积函数,用三角多项式)(xf2nxbnxaxbxaaxSnnnsincossincos21)(110(6.1)作为最佳平方逼近函数.由于三角函数族 kx,kx,x,x,sincossincos1在 上是正交函数族,于是 在 上的最小平方三角逼近
19、多项式 的系数是 2,02,0)(xf)(xSn34 称为傅里叶系数.kkba,函数 按傅里叶系数展开得到的级数)(xf10)sincos(21kkkkxbkxaa(6.3)就称为傅里叶级数.20dcos)(1xkxxfak),1,0(nk(6.2)),1,0(nk20dsin)(1xkxxfbk35 只要 在 上分段连续,则级数(6.3)一致收敛到 .)(xf 2,0)(xf 对于最佳平方逼近多项式(6.1)有.)()()()(222222xSxfxSxfnn由此可以得到相应于(4.11)的贝塞尔不等式.d)(1)(2120212220 xxfbaankkk因为右边不依赖于 ,左边单调有界,
20、所以级数 n10)sincos(21kkkkxbkxaa(6.3)nxbnxaxbxaaxSnnnsincossincos21)(110(6.1).)()(22122*xfxankkk(4.11)36 当 只在给定的离散点集 )(xf1,1,0,2NjjNxj上已知时,则可类似得到离散点集正交性与相应的离散傅里叶系数.下面只给出奇数个点的情形.12220)(21kkkbaa收敛,并有.0limlimkkkkba37122mjxj),2,1,0(mj可以证明对任何 成立 mlk,0令.,0,0sincos;0,12,0212,0coscos;0,212,0,0sinsin202020mjkkxl
21、xklmklmklkxlxklmklklkxlxmjjjmjjjmjjj38这表明函数族 在点集sincossincos1mxmx,x,x,122mjxj上正交.若令),2,1,0()(mjxffjj,),sincos(21)(10mnkxbkxaaxSnkkkn其中 则 的最小二)(xf乘三角逼近为),1,0(122cos12220mkmjkfmamjjk39当 时 mn,)2,1,0()(mjfxSjjm于是(6.4)).,1(122sin12220nkmjkfmbmjjkmkkkmkxbkxaaxS10)sincos(21)(就是三角插值多项式,系数仍由(6.4)表示.40由于),1i,
22、1,1,0()sin(i)cos(eiNjjxjxjx所以函数族 在区间 上是正交的.e,e,1)1(iixNx2,0 一般情形,假定 是以 为周期的复函数,给定)(xf2在 个等分点 上的值)1,1,0(2NjjNxj,2jNffjN.)e,e,1(T)1(2i2iNNjNjj函数 在等距点集 上的值jxie)1,1,0(2NkkNxkkjxie组成的向量记作41102i2ilee),(NkkNskNls当 时,个复向量 具有如下正交性:1,1,0Nj110,NN102)i(eNkkNsl(6.5).,;,0slNsl42事实上,令,e2)(iNslr,0)1(,10sNNl于是,1)1(N
23、slN即;1111NNNslNN若,0 sl.1e2)(i slNr,1,0Nsl若则有,1r则从而43于是 10l),(NkksrrrN11.0若,sl 10),(Nkkssr这就证明了(6.5)成立.即 是正交的.110,N,1r则于是.N 因此,在 个点 上的最小二乘傅里叶逼近为)(xf)1,1,0(2NjjNxjN),(ls(6.5).,;,0slNsl44,e)(10iNncxSnkkxk(6.6)其中.1,1,0,e1102i-nkfNcNjNkjjk(6.7)在(6.6)中,若 ,Nn 则 为 在点)(xS)(xf)1,1,0(Njxj上的插值函数,于是由(6.6)得),()(j
24、jxfxS即.1,1,0,e102iNjcfNkjNkkj(6.8)45 而(6.8)是由 求 的过程,称为反变换反变换.kcjf(6.7)是由 求 的过程,jfkc称为 的离散离散)(xf傅里叶变换傅里叶变换.简称DFT,.1,1,0,e1102i-nkfNcNjNkjjk(6.7).1,1,0,e102iNjcfNkjNkkj(6.8)46 3.6.2 快速傅氏变换(快速傅氏变换(FFT)不论是按(6.7)式由 求 ,jfkc由 求 ,kcjf,1,1,0,10NjwxcNkkjkj(6.9)其中 (正变换)/2iexp(Nw或 (反变换),)/2iexp(Nw,kkba还是由(6.4)计
25、算傅里叶逼近系数都可归结为计算)1,1,0(Nkxk是已知复数序列.或是按(6.8)mjjkmjkfma20122cos122(6.4)mjjkmjkfmb20122sin12247 当 较大且处理数据很多时,就是用高速的电子计算机,很多实际问题仍然无法计算,N 如直接用(6.9)计算 ,需要 次复数乘法和 次jcNN复数加法,称为 个操作,计算全部 共要 个操作.Njc2N 直到20世纪60年代中期产生了FFT算法,大大提高了运算速度,从而使傅氏变换得以广泛应用.FFT算法的基本思想就是尽量减少乘法次数.,1,1,0,10NjwxcNkkjkj(6.9)48用(6.9)计算全部 ,jc表面看
26、要做 个乘法,2N实际上所有 中,1,1,0,),/2iexp(NkjNkj只有 个不N,110Nwww同的值特别当 时,只有 个不同的值.pN22/N 因此可把同一个 对应的 相加后再乘 ,这就能大量减少乘法次数.rwkxrw,1,1,0,10NjwxcNkkjkj(6.9)49 设正整数 除以 后得商 及余数 ,Nmqr则 ,rqNm称为 的 同余数,以 表示.rmNrmN 由于,1),/2iexp(2iewNwN 因此计算 时可用 的 同余数 代替 ,从而推出FFT算法.mwwNrm 以 为例.说明FFT的计算方法.32N 由于 则(6.9)的和是,121,03Njk.7,1,0,70j
27、wxckkjkj(6.10).)(rrqNmwwww故有50将 用二进制表示为 jk,),(222012001122kkkkkkk其中 只能取0或1.)2,1,0(,rjkrr 例如).110(20226022根据 表示法,有jk,).(),(012012kkkxxjjjcckj公式(6.10)可表示为);(222012001122jjjjjjj,70kkjkjwxc(6.10)51 101010)222)(012012012001122012)()(kkkjjjkkkwkkkxjjjc(6.11).)()00(10)0(1010)(012020011120120kjkkkjkkkkkjwww
28、kkkx若引入记号),()(0120120kkkxkkkA,)()(10)00(01020123002kkjwjjkAjjjA(6.12),)()(10)(0120001120120kkkkjwkkkAjkkA,)()(10)0(001101021011kkkjwjkkAjjkA52则(6.11)变成).()(0123012jjjAjjjc 它说明利用 同余数可把计算 分为 步,用公式(6.12)计算.Njcp 每计算一个 只用2次复数乘法,计算一个 用 qAjcp2次复数乘法,计算全部 共用 次复数乘法.jcpN2 若注意 公式(6.12)还可进一步简化为,)1(2/2010jNjjwwp5
29、310)(0120001120120)()(kkkkjwkkkAjkkA),1()0()0(010010011kkAkkAkkA)0(2010)0(01001020010)1()0(kkjjkkjwwkkAwkkA,)1()1()0()0(0100100100kkjjwkkAkkA.)1()0()1()0(01001001101kkwkkAkkAkkA将这表达式中二进制表示还原为十进制表示:,22)0(001101kkkkk54),2()()2(2001kAkAkA).3,2,1,0()2()()12(2001kwkAkAkAk(6.13)同样(6.12)中的 也可简化为 2A,)1()1()
30、0()()00(0010010102011kjjwjkAjkAjjkA即),1()0()0(001001002jkAjkAjkA即 得,3,2,1,0k55.)1()0()1()00(0000010020kwjkAjkAjkA把二进制表示还原为十进制表示,得),22()2()2(21122jkAjkAjkA).1,0;1,0()22()2()22(221122jkwjkAjkAjkAk(6.14)同理(6.12)中 可简化为 3A),1()1()0()(01201201232jjAjjAjjjAj即 56),1()0()0(012012013jjAjjAjjA表示为十进制,有).1()0()1
31、(012012013jjAjjAjjA),2()()(2223jAjAjA(6.15)).3,2,1,0()2()()2(22223jjAjAjA57 根据公式(6.13),(6.14),(6.15),由7,1,0k,)()(0kxkxkA),7,1,0()(3jcjAj逐次计算到见表3-2(略).上面推导的 的计算公式可类似地推广到 的情形.32NpN2 根据公式(6.13),(6.14),(6.15),一般情况的FFT计算公式如下:58(6.16),)22()2(1211111qkpqqqqwjkAjkA),22()2()2(11111pqqqqqqjkAjkAjkA其中.12,1,0;1
32、2,1,0;,11qqpjkpq从 出发,由 到 算到)1,1,0)(0NmmAq1p 一组 占用 个复数单元,计算时需给出两组单元,qAN)22(1qqqjkAqA括号内的数代表它的位置,在计算机中代表存放数的地址.),1,1,0()(NjcjAjp即为所求.59 这个计算公式除了具有不倒地址的优点外,计算只有两重循环,计算过程中只要按地址号存放 则最后得到的,qA)(jAp就是所求离散频谱的次序.外循环 由 计算到 ,内循环 由 计算到q1pk0,12 qp由 计算到j0,121q更重要的是整个计算过程省计算量.由公式看到算一个 共做 次复数乘法,qA2/221/Nqqp而最后一步计算 时
33、,由于pAkNkwwp)(2/21k)1(1)1(060 当 时比值是 它比一般FFT的计算量(次乘法)也快一倍.102N,1:2305.4:1024pN(注意 时 故 ),因此,总共要算pq,012 qp0k次复数乘法,它比直接用(6.9)需 次乘法2/)1(Np 2N.2/)1(:pN快得多,计算量比值是 我们称(6.16)的计算公式为改进的改进的FFT算法算法.613.7 有有 理理 逼逼 近近 3.7.1 有理逼近与连分式有理逼近与连分式 有理函数逼近是指用形如)()()(xQxPxRmnnm的函数逼近 ).(xf 与前面讨论一样,如果 最小就可得到最佳有理一致逼近最佳有理一致逼近.)
34、()(xRxfnm(7.1)mkkknkkkxbxa0062 如果 最小则可得到最佳有理平方逼近最佳有理平方逼近函数函数.2)()(xRxfnm 本节主要讨论利用函数的泰勒展开获得有理逼近函数的方法.对函数 用泰勒展开得)1ln(x.1,1,)1()1ln(11xkxxkkk(7.2)取部分和 nkkknkxxS11)1()().1ln(x63 另一方面,若对(7.2)式用辗转相除可得到 的)1ln(x一种连分式展开 524231211)1ln(22xxxxxx(7.3).52423121122xxxxx.1,1,)1()1ln(11xkxxkkk(7.2)64.612054084042025
35、260630420)(43243244xxxxxxxxxR(7.4)(7.3)右端为 的无穷连分式的前5项,最后式子)1ln(x 若取(7.3)的前2,4,6,8项,则可分别得到 的以下有理逼近:)1ln(x是它的紧凑形式.,22)(11xxxR,6636)(2222xxxxxR,3369060116060)(323233xxxxxxxR65 若用同样多项的泰勒展开部分和 逼近 )(2xSn)1ln(x并计算 处的值 及 ,计算结果见表3-2.1x)1(2 nS)1(nnR00000076.069314642.0058.0634.04000025.0693122.0076.0617.03000
36、84.069231.011.058.02026.0667.019.050.01)1()2ln()1()1()2ln()1(222nnRnnnsnRRSSn表32ln,69314718.0的准确值为从表3-1可以看出,,69314642.0)1(44R,634.0)1(8S66 但它们的计算量是相当的,这说明用有理逼近比多项式逼近好得多.由此看出 的精度比 高出近10万倍,)1(44R)1(8S 例例9 9,4091572115111353381452)(2323443xxxxxxxxR用辗转相除法将它化为连分式并写成紧凑形式.解解给出有理函数用辗转相除可逐步得到6740915721284644
37、32)(23243xxxxxxxR 本例中用连分式计算 的值只需3次除法,1次乘法和7次加法.)(43xR7116)9(654322xxxxx98765432xxxx.98765432xxxx68 若直接用多项式计算的秦九韶算法则需6次乘法和1次除法及7次加法.可见将 化成连分式可节省计算乘除法次数.)(xRnm 对一般的有理函数,(7.1)可转化为一个连分式.)()(121llnmdxcdxcxPxR它的乘除法运算只需 次.),max(nm而直接用有理函数(7.1)计算乘除法次数为 次.mn69 3.7.2 帕德逼近帕德逼近 利用函数 的泰勒展开可以得到它的有理逼近.)(xf 设 在 的泰勒
38、展开为)(xf0 x.)!1()()0(!1)(1)1(0)(NNNkkkxNfxfkxf(7.5)它的部分和记作 NkkkxfkxP0)()0(!1)((7.6).0Nkkkxc70 定义定义1111设,),()(1mnNaaCxfNmmnnnmxbxbxaxaaxR1101)(其中 无公因式,且满足条件)(),(xQxPmn),1,0()0()0()()(NkfRkknm(7.8)则称 为函数 在 处的 阶帕德逼近帕德逼近,)(xRnm)(xf0 x),(mn记作 ,简称 的帕德逼近.),(mnR),(mnR如果有理函数(7.7),)()(xQxPmn71 根据定义,若令),()()()(
39、xPxQxPxhnm则满足条件(7.8)等价于.,1,0,0)0()(Nkhk即.,1,0,0)()()()0(0)()(NkxPxQxPhxknmk由于 应用莱布尼茨求导公式得,!)0()(kknakPkkjjkjxknmakbckxPxQxP!)()()(00)(,1,0Nk,0),1,0()0()0()()(NkfRkknm(7.8)72这里 是由(7.6)得到的,)0(!1)(jjfjc上式两端除 ,!k并由 可得),(0,10时当mjbbjnkcbcakkjjkjk,1,0,10(7.9)及 mnnkcbckkjjkj,1,10(7.10)注意当 时mj,0jb故(7.10)可写成N
40、kkkxfkxP0)()0(!1)((7.6).0Nkkkxc73,11222112211211mnmnmnmnnnnmmnnnnmmncbcbcbccbcbcbccbcbcbc(7.11)其中 时 ,0j,0jc若记,121211mnmnnnnmnnnmncccccccccH(7.12),),(T11bbbmmb.),(T21mnnncccc74则方程组(7.11)的矩阵形式为.cbH 定理定理1010(7.7)的有理函数 是 的 阶帕德逼近的)(xRnm)(xf),(mn充分必要条件是多项式 的系数 )(),(xQxPmnnaaa,10及 满足方程组(7.9)及(7.11).mbbb,10
41、,),()(1mnNaaCxfN设则形如75)4,4()3,4()2,4()1,4()0,4(4)4,3()3,3()2,3()1,3()0,3(3)4,2()3,2()2,2()1,2()0,2(2)4,1()3,1()2,1()1,1()0,1(1)4,0()3,0()2,0()1,0()0,0(0432103-3nm表 根据定理10,求 的帕德逼近时,)(xf首先要由(7.11)解出 的系数 ,)(xQmmbbb,10再由(7.9)直接算出 的系数 .)(xPnnaaa,10 的各阶帕德逼近可列成一张表,称为帕德表(见表3-3).)(xf,11222112211211mnmnmnmnnn
42、nmmnnnnmmncbcbcbccbcbcbccbcbcbc(7.11)),1,0(10nkcbcakkjjkjk(7.9)76 例例1010求 的帕德逼近 及 .)1ln()(xxf)2,2(R)3,3(R 解解由 的泰勒展开)1ln(x432413121)1ln(xxxxx得.,41,31,21,1,043210ccccc 当 时,由(7.11)得 2 mn.413121,31211212bbbb求得,61,121bb再由(7.9)得77,21,1,0210aaa于是得 222261121)(xxxxxR当 时,由(7.11)得 3 mn,61514131,51413121,413121
43、123123123bbbbbbbbb.663622xxxx78代入(7.9)得.6011,1,1,03210aaaa解得.201,53,23321bbb于是得 323233201532316011)(xxxxxxxR.33690601160603232xxxxxx79 为了求帕德逼近 的误差估计,由(7.9)及(7.11)求得的 系数 及 ,直接代入则得)(xRnm)(),(xQxPmnnaaa,10mbbb,10,)()()()(0011llmkklmnkmnnmxcbxxPxQxf 将 除上式两端,即得)(xQm可以看到这里得到的 及 与 的前面)(22xR)(33xR)1ln(x连分式展开得到的有理逼近(7.4)结果一样.80,)()()(01xQxrxxRxfmlllmnnm(7.13)其中.01mkklmnklcbr 当 时可得误差近似表达式 1x.,)()(01010mkkmnkmnnmcbrxrxRxf81