1、1 3.5.2 3.5.2 用正交多项式做最小二乘拟合用正交多项式做最小二乘拟合 如果 是关于点集)(,),(),(10 xxxn带权 正交的), 1 ,0(mixi), 1 ,0()(mixi,0,0)()()(),(0kmiikijikjAxxx,kj ,kj （5.8） 用最小二乘法得到的法方程组（5.6），其系数矩阵 是病态的. G函数族，即2miikimiikiikkkkxxxxfxfa020*)()()()()(),(),()., 1 ,0(nk（5.9）则方程（5.6）的解为 且平方误差为 .)(02*2222nkkkaAf3 接下来根据给定节点 及权函数 mxxx,10, 0)
2、(x构造带权 正交的多项式 .)( x)(xPn 注意 ，用递推公式表示 ，即mn )(xPk)()()()(),()()(, 1)(1110110 xPxPxxPxPxxPxPkkkkk).1,2, 1(nk（5.10）根据 的)(xPk这里 是首项系数为1的 次多项式，)(xPkk正交性，得4miikimiikiikxPxxPxx02021)()()()(),(),(11kkkkPPPP（5.11）).1,2, 1(nk 下面用归纳法证明这样给出的 是正交的. )(xPk)(),()(),(xPxPxPxxPkkkk),(),(kkkkPPPxPmiikimiikikxPxxPx02102
3、)()()()(5),(),(),(0010010PPxPPPP 假定 对 及)(0),(slPPsl1, 1 ,0ls;,1 ,0kl要证 对 均成立. 0),(1skPPks,1 ,0由（5.10）有 ),(),)(),(111skkskkskPPPPxPP 由（5.10）第二式及（5.11）中 的表达式，有 1),(),(),(),(00000000PPPPPxPxPP.0nk 均成立，（5.12）).,(),(),(11skkskkskPPPPPxP6 而 ，11ks,0),(),(skskxPPPxP于是由（5.12），当 时， 2 ks.0),(1skPP 另外， 是首项系数为1的
4、 次多项式，它可由)(xxPs1s由归纳法假定，当 时20ks,0),(slPP.0),(1skPP110,sPPP的线性组合表示.由归纳法假定又有7由假定有 ),(),(11kkkkxPPPxP 再考虑 （5.13）),(),(),(),(1111111kkkkkkkkkkPPPPPxPPP),(10kjjjkkPcPP).,(kkPP利用（5.11）中 表达式及以上结果，得 k),(),(),(11111kkkkkkkPPPxPPP.0),(),(kkkkPPPP8),(),(),(),(111kkkkkkkkkkPPPPPxPPP至此已证明了由（5.10）及（5.11）确定的多项式 )(
5、xPk组成一个关于点集 的正交系.ix 用正交多项式 的线性组合作最小二乘曲线拟合，)(xPk只要根据公式（5.10）及（5.11）逐步求 的同时，)( xPk相应计算出系数*ka.0),(),(),(),(kkkkkkkkPPPPPxPPxP最后，由 和 的表达式（5.11）有 kk9),(),(*kkkkPPPfa),1 ,0(nk并逐步把 累加到 中去，最后就可得到所求的 )(*xPakk)( xS).()()()(*1*10*0 xPaxPaxPaxSynn 用这种方法编程序不用解方程组，只用递推公式，并且当逼近次数增加一次时，只要把程序中循环数加1，其余不用改变. 这里 可事先给定或
6、在计算过程中根据误差确定. n)()()()()(200ikmiimiikiixPxxPxfx拟合曲线103.6 3.6 最佳平方三角逼近与快速傅里叶变换最佳平方三角逼近与快速傅里叶变换 当 是周期函数时，显然用三角多项式逼近 比用代数多项式更合适，本节主要讨论用三角多项式做最小平方逼近及快速傅里叶变换，快速傅里叶变换，简称FFT算法. )( xf)( xf11 3.6.1 3.6.1 最佳平方三角逼近与三角插值最佳平方三角逼近与三角插值 设 是以 为周期的平方可积函数，用三角多项式 )( xf2nxbnxaxbxaaxSnnnsincossincos21)(110（6.1）做最佳平方逼近函数
7、. 由于三角函数族 kx,kx,x,x,sincossincos1在 上是正交函数族，于是 在 上的最小平方三角逼近多项式 的系数是 2,02,0)( xf)(xSn12 称为傅里叶系数. kkba, 函数 按傅里叶系数展开得到的级数 )( xf10)sincos(21kkkkxbkxaa（6.3）就称为傅里叶级数.20cos)(1kxdxxfak),1 ,0(nk（6.2）),1 ,0(nk20sin)(1kxdxxfbk13 只要 在 上分段连续，则级数（6.3）一致收敛到 . )(xf 2,0)(xf 对于最佳平方逼近多项式（6.1）有 .)()()()(222222xSxfxSxfnn
8、由此可以得到相应于（4.11）的贝塞尔不等式 .)(1)(2120212220dxxfbaankkk因为右边不依赖于 ，左边单调有界，所以级数 n14 当 只在给定的离散点集 )( xf1, 1 ,0,2NjjNxj上已知时，则可类似得到离散点集正交性与相应的离散傅里叶系数. 下面只给出奇数个点的情形. 12220)(21kkkbaa收敛，并有 .0limlimkkkkba15122mjxj),2, 1 ,0(mj可以证明对任何 成立 mlk,0令.,0,0sincos;0, 12;0212;,0coscos;0,212;0,0sinsin202020mjkkxlxklmklmklkxlxkl
9、mklklkxlxmjjjmjjjmjjj16这表明函数族 在点集sincossincos1mxmx,x,x,122mjxj上正交. 若令),2, 1 ,0()(mjxffjj则 的最小二)( xf,),sincos(21)(10mnkxbkxaaxSnkkkn其中 乘三角逼近为), 1 ,0(122cos12220mkmjkfmamjjk17当 时 mn ,)2,1 ,0()(mjfxSjjm于是 （6.4）).,1(122sin12220nkmjkfmbmjjkmkkkmkxbkxaaxS10)sincos(21)(就是三角插值多项式，系数仍由（6.4）表示. 18由于 ),1i, 1,
10、1 , 0()sin(i)cos(eiNjjxjxjx所以函数族 在区间 上是正交的. e,e,1)1(iixNx2,0 一般情形，假定 是以 为周期的复函数，给定 )( xf2在 个等分点 上的值)1, 1 ,0(2NjjNxj,2jNffjN函数 在等距点集 上的值jxie)1,1 ,0(2NkkNxk.)e,e,1()1(2i2ijTNNjNjkjxie组成的向量记作19102i2ilee),(NkkNskNls当 时， 个复向量 具有如下正交性： 1,1 ,0Nj110,NN102)i(eNkkNsl（6.5）.,;,0slNsl20事实上，令,e2)(iNslr,0)1(,10sNN
11、l于是 ,1)1(NslN即 ;1111NNNslNN若,0 sl;1e2)(i slNr,1,0Nsl若则有,1r则从而21于是 10l),(NkksrrrN11.0若,sl 10s),(Nkksr这就证明了（6.5）成立. 即 是正交的. 110,N,1r则于是.N 因此， 在 个点 上的最小二乘傅里叶逼近为 )( xf)1, 1 ,0(2NjjNxjN22,e)(10iNncxSnkkxk（6.6）其中 ).1,1 ,0(,e1102i -nkfNcNjNkjjk（6.7）在（6.6）中，若 ，Nn 则 为 在点)( xS)( xf)1,1 ,0(Njxj上的插值函数，于是由（6.6）得
12、 ),()(jjxfxS即).1,1 ,0(,e102iNjcfNkjNkkj（6.8）23（6.7）是由 求 的过程，jfkc称为 的离散离散)( xf 而（6.8）是由 求 的过程，称为反变换反变换.kcjf傅里叶变换傅里叶变换. 简称DFT，24 3.6.2 3.6.2 快速傅氏变换（快速傅氏变换（FFT） 不论是按（6.7）式由 求 ，jfkc由 求 ，kcjf),1,1 ,0(,10NjwxcNkkjkj（6.9）其中 (正变换) )/2iexp(Nw或 (反变换)，)/2iexp(Nw ,kkba还是由（6.4）计算傅里叶逼近系数都可归结为计算)1,1 ,0(Nkxk是已知复数序列
13、.或是按（6.8）25 当 较大且处理数据很多时，就是用高速的电子计算机，很多实际问题仍然无法计算， N 如直接用（6.9）计算 ，需要 次复数乘法和 次jcNN复数加法，称为 个操作，计算全部 共要 个操作. Njc2N 直到20世纪60年代中期产生了FFT算法，大大提高了运算速度，从而使傅氏变换得以广泛应用. FFT算法的基本思想就是尽量减少乘法次数.26用（6.9）计算全部 ，jc表面看要做 个乘法，2N实际上所有 中，1, 1 , 0,),/2iexp(NkjNkj只有 个不N,110Nwww同的值特别当 时，只有 个不同的值.pN22/N 因此可把同一个 对应的 相加后再乘 ，这就能
14、大量减少乘法次数.rwkxrw27 设正整数 除以 后得商 及余数 ，Nmqr则 ，rqNm称为 的 同余数，以 表示. rmNrmN 由于, 1),/2iexp(2iewNwN 因此计算 时可用 的 同余数 代替 ，从而推出FFT算法. mwwNrm 以 为例. 说明FFT的计算方法. 32N 由于 则（6.9）的和是 ,121,03Njk).7,1 ,0(,70jwxckkjkj（6.10）.)(rrqNmwwww故有28将 用二进制表示为 jk ,),(222012001122kkkkkkk其中 只能取0或1， )2,1 ,0(,rjkrr 例如).110(20226022根据 表示法，
15、有jk ,).(),(012012kkkxxjjjcckj公式（6.10）可表示为 );(222012001122jjjjjjj29 101010)222)(012012012001122012)()(kkkjjjkkkwkkkxjjjc（6.11）)00(10)0(1010)(012020011120120)(kjkkkjkkkkkjwwwkkkx若引入记号 ),()(0120120kkkxkkkA.)()(10)00(01020123002kkjwjjkAjjjA（6.12）,)()(10)(0120001120120kkkkjwkkkAjkkA,)()(10)0(001101021011
16、kkkjwjkkAjjkA30则（6.11）变成 ).()(0123012jjjAjjjc 它说明利用 同余数可把计算 分为 步，用公式（6.12）计算. Njcp 每计算一个 只用2次复数乘法，计算一个 用 qAjcp2次复数乘法，计算全部 共用 次复数乘法.jcpN2 若注意 公式（6.12）还可进一步简化为 ,) 1(2/2010jNjjwwp3110)(0120001120120)()(kkkkjwkkkAjkkA),1()0()0(010010011kkAkkAkkA)0(2010)0(01001020010)1()0(kkjjkkjwwkkAwkkA,)1()1()0()0(010
17、0100100kkjjwkkAkkA.)1()0()1()0(01001001101kkwkkAkkAkkA将这表达式中二进制表示还原为十进制表示：,22)0(001101kkkkk32),2()()2(2001kAkAkA).3,2,1 ,0()2()()12(2001kwkAkAkAk（6.13）同样（6.12）中的 也可简化为 2A,)1()1()0()()00(0010010102011kjjwjkAjkAjjkA即 ),1()0()0(001001002jkAjkAjkA即 得,3,2, 1 ,0k33.)1()0()1()00(0000010020kwjkAjkAjkA把二进制表示
18、还原为十进制表示，得 ),22()2()2(21122jkAjkAjkA).1 ,0;1 ,0()22()2()22(221122jkwjkAjkAjkAk（6.14）同理（6.12）中 可简化为 3A),1()1()0()(01201201232jjAjjAjjjAj即 34),1()0()0(012012013jjAjjAjjA表示为十进制，有 ).1()0()1(012012013jjAjjAjjA),2()()(2223jAjAjA（6.15）).3,2,1 ,0()2()()2(22223jjAjAjA35 根据公式(6.13)，(6.14)，(6.15)，由)7,1 ,0(kkxk
19、xkA)()(0),7,1 ,0()(3jcjAj逐次计算到见表3-2(略）. 上面推导的 的计算公式可类似地推广到 的情形. 32NpN2 根据公式(6.13)，(6.14)，(6.15)，一般情况的FFT计算公式如下： 36（6.16）,)22()2(1211111qkpqqqqwjkAjkA),22()2()2(11111pqqqqqqjkAjkAjkA其中 .12,1 ,0;12,1 ,0;,11qqpjkpq从 出发， 由 到 算到 )1,1 ,0)(0NmmAq1p 一组 占用 个复数单元，计算时需给出两组单元，qAN)22(1qqqjkAqA括号内的数代表它的位置，在计算机中代表
20、存放数的地址.),1,1 ,0()(NjcjAjp即为所求.37 这个计算公式除了具有不倒地址的优点外，计算只有两重循环， 计算过程中只要按地址号存放 则最后得到的 ,qA)( jAp就是所求离散频谱的次序.外循环 由 计算到 ，内循环 由 计算到q1pk0,12 qp由 计算到j0,121q更重要的是整个计算过程省计算量. 由公式看到算一个 共做 次复数乘法，qA2/221/Nqqp而最后一步计算 时，由于pAkNkwwp)(2/21k)1(1)1(038（注意 时 故 ），因此，总共要算pq ,012 qp0k次复数乘法，它比直接用（6.9）需 次乘法2/)1(Np 2N 当 时比值是 它
21、比一般FFT的计算量（ 次乘法）也快一倍. 102N,1:2305.4:1024pN.2/)1(:pN快得多，计算量比值是 我们称（6.16）的计算公式为改进的改进的FFT算法算法 .393.7 3.7 有有 理理 逼逼 近近 3.7.1 3.7.1 有理逼近与连分式有理逼近与连分式 有理函数逼近是指用形如 )()()(xQxPxRmnnm的函数逼近 ).( xf 与前面讨论一样，如果 最小就可得到最佳有理一致逼近最佳有理一致逼近. . )()(xRxfnm（7.1）mkkknkkkxbxa0040 如果 最小则可得到最佳有理平方逼近最佳有理平方逼近函数函数. 2)()(xRxfnm 本节主要
22、讨论利用函数的泰勒展开获得有理逼近函数的方法. 对函数 用泰勒展开得 )1ln(x.1 , 1)1()1ln(11xkxxkkk（7.2）取部分和 nkkknkxxS11)1()().1ln(x41 另一方面若对（7.2）式用辗转相除可得到 的)1ln(x一种连分式展开 524231211)1ln(22xxxxxx（7.3）.52423121122xxxxx42.612054084042025260630420)(43243244xxxxxxxxxR（7.4）（7.3）右端为 的无穷连分式的前5项，最后式子)1ln(x 若取（7.3）的前2，4，6，8项，则可分别得到 的以下有理逼近 )1ln
23、(x是它的紧凑形式. ,22)(11xxxR,6636)(2222xxxxxR,3369060116060)(323233xxxxxxxR43 若用同样多项的泰勒展开部分和 逼近 )(2xSn)1ln(x并计算 处的值 及 ，计算结果见表3-3.1x)1(2 nS)1(nnR00000076.069314642.0058.0634.04000025.0693122.0076.0617.0300084.069231.011.058.02026.0667.019.050.01)1()2ln()1()1()2ln()1(22nnRnnnsnRRSSn3表32ln,69314718.0的准确值为从表3
24、-1可以看出，,69314642.0)1 (44R,634.0)1 (8S44 但它们的计算量是相当的，这说明用有理逼近比多项式逼近好得多. 由此看出 的精度比 高出近10万倍，)1(44R)1(8S 例例9 94091572115111353381452)(2323443xxxxxxxxR用辗转相除法将它化为连分式并写成紧凑形式. 解解给出有理函数用辗转相除可逐步得到454091572128464432)(23243xxxxxxxR 本例中用连分式计算 的值只需3次除法，1次乘法和7次加法. )(43xR7116)9(654322xxxxx98765432xxxx.98765432xxxx4
25、6 若直接用多项式计算的秦九韶算法则需6次乘法和1次除法及7次加法. 可见将 化成连分式可节省计算乘除法次数. )( xRnm 对一般的有理函数（7.1）可转化为一个连分式.)()(121llnmdxcdxcxPxR它的乘除法运算只需 次.),max(nm而直接用有理函数（7.1）计算乘除法次数为 次. mn47 3.7.2 3.7.2 帕德逼近帕德逼近 利用函数 的泰勒展开可以得到它的有理逼近. )( xf 设 在 的泰勒展开为 )( xf0 x.)!1()()0(!1)(1)1(0)(NNNkkkxNfxfkxf（7.5）它的部分和记作 NkkkxfkxP0)()0(!1)(（7.6）.0
26、Nkkkxc48 定义定义1111设,),()(1mnNaaCxfNmmnnnmxbxbxaxaaxR1101)(其中 无公因式，且满足条件)(),(xQxPmn),1 ,0()0()0()()(NkfRkknm（7.8）则称 为函数 在 处的 阶帕德逼近帕德逼近，)( xRnm)( xf0 x),(mn记作 ,简称 的帕德逼近.),(mnR),(mnR如果有理函数（7.7）,)()(xQxPmn49 根据定义，若令 ),()()()(xPxQxPxhnm则满足条件（7.8）等价于 .,1 ,00)0()(Nkhk即 ., 1 ,0,0)()()()0(0)()(NkxPxQxPhxknmk由
27、于 应用莱布尼兹求导公式得 ,!)0()(kknakPkkjjkjxknmakbckxPxQxP!)()()(00)(,1 ,0Nk050这里 是由（7.6）得到的，)0(!1)( jjfjc上式两端除 ，!k并由 可得 ),(0,10时当mjbbj),1 ,0(10nkcbcakkjjkjk（7.9）及 ),1(10mnnkcbckkjjkj（7.10）注意当 时mj ,0jb故（7.10）可写成51,11222112211211mnmnmnmnnnnmmnnnnmmncbcbcbccbcbcbccbcbcbc（7.11）其中 时 ，0j,0jc若记,121211mnmnnnnmnnnmnc
28、ccccccccH（7.12）,),(11Tmmbbbb.),(21Tmnnncccc52则方程组（7.11）的矩阵形式为 .cbH 定理定理1010（7.7）的有理函数 是 的 阶帕德逼近的)( xRnm)( xf),(mn充分必要条件是多项式 的系数 )(),(xQxPmnnaaa,10及 满足方程组（7.9）及（7.11）. mbbb,10,),()(1mnNaaCxfN设则形如53 根据定理10, 求 的帕德逼近时，首先要由（7.11）)( xf解出 的系数 ，)( xQmmbbb,10的系数 .)( xPnnaaa,10)( xf的各阶帕德逼近可列成再由（7.9）直接算出一张表，称为
29、帕德表（见表3-4）. )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(0432104-3nm表54 例例1010求 的帕德逼近 及 . )1ln()(xxf)2,2(R)3,3(R 解解由 的泰勒展开)1ln(x432413121)1ln(xxxxx得 .,41,31,21,1,043210ccccc 当 时，由（7.11）得 2 mn.413121,31211212bbbb
30、求得,61,121bb再由（7.9）得55,21,1,0210aaa于是得 222261121)(xxxxxR当 时，由（7.11）得 3 mn,61514131,51413121,413121123123123bbbbbbbbb.663622xxxx56代入（7.9）得 .6011, 1, 1,03210aaaa解得 .201,53,23321bbb于是得 323233201532316011)(xxxxxxxR.33690601160603232xxxxxx57可以看到这里得到的 及 与 的前面 )(22xR)(33xR)1ln(x 为了求帕德逼近 的误差估计，由(7.9)及(7.11)求得的 系数 及 ，直接代入则得 )( xRnm)(),(xQxPmnnaaa,10mbbb,10,)()()()(0011llmkklmnkmnnmxcbxxPxQxf 将 除上式两端，即得 )( xQm连分式展开得到的有理逼近（7.4）结果一样.58,)()()(01xQxrxxRxfmlllmnnm（7.13）其中 .01mkklmnklcbr 当 时可得误差近似表达式 1x.,)()(01010mkkmnkmnnmcbrxrxRxf
