第3章 函数逼近与快速傅里叶变换 3.1 函数逼近简介 多项式插值要求插值函数与被插值函数在节点处相等,在非插值节点处误差有可能比较大.函数逼近与曲线拟合则是从整体上考虑,强调按照某种范数度量一个函数逼近另一个函数,它们都是近似计算的基础,在实际问题的计算中亦有重要应用.在函数逼近理论中,最简单的是最佳平方逼近,包括连续最佳平方逼近和离散最佳平方逼近,后者亦称为曲线拟合的最小二乘法. 设f(x),g(x)∈C\,ρ(x)为\上的权函数,定义f(x),g(x)的内积为(f,g)=∫baρ(x)f(x)g(x)dx,内积范数为‖f‖2=(f,f),‖f‖2也称为f(x)的平方范数.假设φ0(x),φ1(x),…,φn(x)∈C\且在\上线性无关,令φ=span{φ0(x),φ1(x),…,φn(x)},则函数的最佳平方逼近问题定义为求S(x)∈φ,使得‖f-S‖2=minS∈φ‖f-S‖2.设最佳平方逼近问题的解S(x)=∑nj=0ajφj(x),则aj(j=0,1,…,n)满足法方程组∑nj=0(φk,φj)aj=(f,φk), k=0,1,…,n.(3.1)直接计算可得最佳平方逼近误差‖f-S‖22=‖f‖22-∑nj=0aj(φj,f).(3.2) 如果{φ0(x),φ1(x),…,φn(x)}为\上带权ρ(x)的正交函数序列,则法方程组的解aj=(f,φj)/‖φj‖22,j=0,1,…,n.常用的正交多项式有勒让德(Legendre)正交多项式(区间\,权函数ρ(x)=1)和切比雪夫(Chebyshev)正交多项式区间\,权函数ρ(x)=11-x2,分别定义为P0(x)=1,Pn(x)=12nn!dndxn(x2-1)n,n≥1和Tn(x)=cosnarccosx. 若f(x)∈C\,使用勒让德正交多项式作基底的n次最佳平方逼近多项式Sn(x)=∑nk=0akPk(x),其中ak=2k+12∫1-1f(x)Pk(x)dx,平方逼近误差为‖f(x)-Sn(x)‖22=∫1-1f(x)2dx-∑nk=022k+1a2k. (3.3)另外,使用切比雪夫正交多项式作基底的n次最佳平方逼近多项式为S~n(x)=c02+∑nk=1ckTk(x),其中ck=2π∫π0f(cosθ)coskθdθ. 若仅知道f(x)在m+1个节点a=x0,x1,…,xm=b上的函数值,则这时的最佳平方逼近称为离散最佳平方逼近,只需将上面的内积定义和范数定义中的求积分改为求和,如(f,g)=∑mi=0ρ(xi)f(xi)g(xi),则只要拟合基函数满足哈尔(Harr)条件\,上面有关连续最佳平方逼近的结论仍成立,求解法方程组得到的S(x)称为给定离散数据的最小二乘拟合曲线.曲线拟合的最小二乘法是一类重要的数学方法,在生产实践中有重要应用,该方法在统计中称为回归分析,常用来对某些现象进行预测和分析,在本章及后续章节中,我们多次用到了这一方法. 函数的多项式最佳平方逼近计算简单,其收敛性有理论保证,在实际应用中对于大部分函数可以得到良好的逼近效果,但对于一些具有特殊性质的函数,使用其他形式的逼近效果更好.实验3.6使用有理函数逼近平方根函数得到了非常好的计算效果.如果函数在给定区间上具有周期性,则使用三角函数逼近效果较好.在三角函数逼近中,经常需要计算给定数据的离散傅里叶(Fourier)变换,通常采用快速傅里叶变换(FFT)方法计算.实验3.7给出了周期函数的三角插值公式和离散傅里叶变换公式,实验3.8直观演示了离散傅里叶变换的一个良好性质,即数据列经离散傅里叶变换后非常有规律,使得FFT在数字信号处理中起着举足轻重的作用. Mathematica提供了一些命令用以实现函数逼近,主要有: LegendreP\,n次勒让德正交多项式,x表示函数自变量. ChebyshevT\,n次切比雪夫正交多项式,x表示函数自变量. FromCharacterCode\,将整数n转换为对应的字符. Fit\,用数据data,以vars为自变量,按基函数funs的形式构造拟合函数,它是一个线性最小二乘拟合函数. FindFit\,非线性拟合函数,data为数据,vars为自变量,expr为拟合函数表达式,其中的参数在pars中体现.若expr关于参数为线性形式,则该函数返回全局最优值,否则,结果为局部最优值. FindFit\,同上,非线性拟合函数,cons对expr中的参数进行限制. Fourier\,求由复数构成的表list的离散傅里叶变换. InverseFourier\,求表list的离散傅里叶逆变换. 3.2 例题选解 例3.1 设f(x)=1+x2,求f(x)在\ 上的一次和二次最佳平方逼近多项式. 解 求最佳平方逼近多项式的关键在于构造法方程组(3.1),使用单项式基底的一次最佳平方逼近多项式程序为 u={1, x}; (u为基向量) m=Outer\ (向量的外积,即列向量乘以行向量) {{1, x}, {x, x^2}} A=Integrate\ (计算法方程组系数矩阵A) {{1, 1/2},{1/2, 1/3}} f\:=Sqrt\x^2\]; (定义函数) d=Integrate\, {x, 0, 1}\];(计算右端向量d) d=N\ (d的小数形式,结果{1.14779,0.609476}) c=LinearSolve\(解向量为c={0.93432, 0.426947}) S\=c.u(求最佳平方逼近,S\=0.93432+0.426947x ) delt2=Integrate\^2, {x, 0, 1}\]-c.d; (公式(3.2)) delt=Sqrt\(计算最佳平方逼近平方范数,为0.0267007) deltinf=NMaximize\-S\\], 0<= x<= 1}, x\] {0.0300434, {x -> 0.472142}} (计算最佳平方逼近最大模范数,结果为区间内的局部极值,真正的最大误差在左端点达到,为0.06568,参见图3.1.) Plot\,S\},{x,0,1},PlotStyle->{{Red,Thick},{Blue,Thick,Dashed}}, AxesLabel->{x,y},AxesStyle->Thick,BaseStyle->{FontSize->16}\] (图3.1,虚线为最佳平方逼近多项式)图3.1 函数(实线)及其一次最佳平方 逼近多项式(虚线)图形 由图3.1知,函数的一次最佳平方逼近多项式精度不高.至于二次最佳平方逼近多项式,只需将上面程序中的基向量u改为u={1,x,x^2}即可,最后得二次最佳平方逼近多项式 S(x)=0.993768+0.0702575x+0.35669x2, 该最佳平方逼近的平方范数误差为0.00247163,逼近效果见图3.2. 进一步绘出f(x)-S(x)的误差图形,如图3.3所示.图3.2 函数(实线)及其二次最佳平 方逼近多项式(虚线)图形 图3.3 f(x)-S(x)的误差图形 f(x)的二次最佳平方多项式逼近效果已相当好,在\ 区间上的误差分布比较均匀. 例3.2 求f(x)=ex在\上的三次最佳平方逼近多项式. 解 此题有两种方法,一是使用例3.1中的方法,需求解法方程组;二是使用勒让德正交多项式逼近,无需形成法方程组.我们采用第2种方法,Mathematica程序为 f\:=Exp\ (定义函数) c=Table\1)/2Integrate\f\,{x,-1,1}\],{k,0,3}\]//N {1.1752,1.10364,0.357814,0.0704556}(3次最佳平方逼近多项式系数) S\=c.Table\,{n,0,3}\]//Expand 0.996294+0.997955x+0.536722x2+0.176139x3(正交逼近多项式使用计算机计算非常简单) delt2=Integrate\^2,{x,-1,1}\]-(c^2).Table\1),{n,0,3}\];(公式(3.3)) delt=Sqrt\ (平方误差为0.00472111) deltinf=NMaximize\-S\\],-1<=x<=1},x\] (求最大绝对误差,结果不准确,应为0.0111723,在右边界点达到)逼近函数和逼近误差图形分别见图3.4和图3.5,其中图3.5还画出了f(x)在x=0的三次泰勒展开式的误差图形(虚线). 图3.4 函数(实线)及其三次最佳平方 逼近多项式(虚线)图形 图3.5 误差图形,实线为f(x)-S(x), 虚线为泰勒展开式的误差图形 总体而言,f(x)的三次最佳平方逼近多项式已有比较高的整体精度表现在图3.4中,函数(实线)与其最佳平方逼近多项式(虚线)的图形完全重合,而由图3.5知,函数的泰勒展开式具有局部性质,仅在展开点附近精度较高. 例3.3 已知一组实验数据如表3.1所示,求它的最小二乘拟合曲线.表3.1 实验数据(表中wi为权数)xi12345yi44.5688.5wi21311 解 最小二乘拟合一般分两部分,一是画出离散点图形,二是根据这些点的走向选择合适的拟合函数类.图3.6 数据点及其线性拟合图形 xp=Range\; fp={4,4.5,6,8,8.5}; wp={2,1,3,1,1} (准备数据) pts=Transpose\(绘制离散点图形要求的数据格式) {{1,4},{2,4.5},{3,6},{4,8},{5,8.5}} p1=ListPlot\> {Red,PointSize\},AxesStyle->Thick, AxesLabel->{x,y},BaseStyle-> {Bold,FontSize->15}\] (绘数据点散点图,见图3.6中的离散点分布)由图3.6知,这些点差不多在一条直线上,所以,选择直线函数类作拟合基底.由于本题中节点权数不为1,不能直接使用函数拟合命令Fit,我们从法方程组的形成求拟合函数. n=Length\; (下面命令形成法方程组的系数矩阵) A=Table\\]xp\\]^ixp\\]^j,{k,1,n}\], {i,0,1},{j,0,1}\] (查看法方程组的具体形式\) {{8,22},{22,74}} b=Table\\]xp\\]^ifp\\],{k,1,n}\], {i,0,1}\] (右端项计算,结果{47.,145.5}) sol=LinearSolve\ (结果为{2.56481,1.2037}) y\=sol.{1,x} (定义线性拟合函数) 2.56481+1.2037x p2=Plot\,{x,xp\\],xp\\]},PlotStyle->{Blue, Thick},AxesStyle->Thick, BaseStyle->Medium\]; Show\(数据点及其线性拟合图形,见图3.6)例3.4 已知一组实验数据如表3.2,求它的拟合曲线.表3.2 实验数据xi11.251.51.752yi5.15.796.537.458.46 解 首先利用给定数据,画出离散点图形,见图3.7,由图形知,这些点差不多在一条指数曲线上,所以,选择y=aexp(bx)作拟合函数.图3.7 数据点及其指数函数拟合图形 对于指数函数,先做对数变换,将其变换为lny=lna+bx,然后再进行线性最小二乘拟合. xp=Range\; yp={5.1,5.79,6.53,7.45,8.46}; pts=Transpose\; (用于绘制散点图) ybar=Log\(对数变换) {1.62924,1.75613,1.87641,2.00821,2.13535} n=Length\; A=Table\\]^ixp\\]^j,{k,1,n}\], {i,0,1},{j,0,1}\](计算法方程组系数矩阵) {{5,7.5},{7.5,11.875}} b=Table\\]^iybar\\],{k,1,n}\],{i,0,1}\] {9.40534,14.4241}(法方程组右端项结果) sol=LinearSolve\ (结果为{1.12249,0.50572}) f\=Exp\\]\]Exp\\]x\] (定义指数拟合函数) 3.07249e0.50572x (指数型拟合公式,见图3.7) err=yp-f\; emax=Sqrt\\] (计算得平方误差0.034727)注 此题可以使用Mathematica提供的FindFit函数进行拟合. sol=FindFit\,{a,b},x\] {a->3.06658,b->0.506955} (输出结果为参数替换形式) f=aExp\/.sol (变量替换) 3.06658e0.506955x (算法不同,结果与上面的拟合函数略有不同)该函数的平方误差0.0341225,略小于上面得到的结果. 3.3 函数的最佳逼近3.3 函数的最佳逼近 实验3.1 函数的伯恩斯坦(Bernstein)多项式逼近. 设f(x)∈C\,构造Bn(x)=∑nk=0fknCknxk(1-x)n-k,其中Ckn为二项式展开系数,则有limn→∞Bn(x)=f(x)在\上一致成立.这里Bn(x)称为f(x)的伯恩斯坦逼近多项式. 实验内容 设f(x)=ex,x∈\,构造其1到20次的伯恩斯坦逼近多项式,计算它们与f(x)的绝对误差,画出误差散点图,并拟合该图形. 函数的伯恩斯坦逼近多项式是一个重要的数学分析工具,在理论上非常完美,该多项式用来计算效果又将如何呢?本实验将回答这一个问题. f\:=Exp\ (定义函数) bf=Table\\]Binomial\x^k(1-x)^(n-k),{k,0,n}\], {n,1,20}\] //Expand//Chop上面语句的作用是计算函数的前20个伯恩斯坦多项式,Binomial\计算二项式系数,Chop的作用是舍去绝对值非常小的项.前三个多项式以表的形式表示如下: {1.+1.71828x,1.+1.29744x+0.420839x2,1.+1.18684x+0.469528x2+0.061917x3,…} err=Table\-bf\\]\],0<=x<=1},x\],{k,Length\}\] (计算前20个多项式与f(x)在\上的最大绝对误差,经绘图知, 最大绝对误差在区间内达到,计算正确) {{0.211867,{x->0.541325}},{0.108049,{x->0.579553}},{0.0724064,{x->0.592396}},…} (误差表,每个子表的第一项为绝对误差最大值,接下来是最大值点) t=Table\\]},{i,Length\}\] (将各次多项式的最大值构成一个用于绘图的表) {{1,0.211867},{2,0.108049},{3,0.0724064},…} p1=ListPlot\>{Red,PointSize\},Frame->True,AxesStyle->Thick, BaseStyle->{FontSize->16}\] (绘制误差散点图) sub=FindFit\(对数据点进行拟合) g=a/n^b/.sub (变量替换,得拟合公式g=0.212244/n0.984094)将数据点和拟合曲线在同一坐标系中绘出,图3.8 函数的1到20次伯恩斯坦逼近 多项式误差拟合图形见图3.8. 由误差图形知,伯恩斯坦多项式逼近收敛趋势明显,但速度比较慢,仅是线性收敛,20次多项式的最大误差为0.01.注意到例3.2中函数的三次最佳平方逼近最大绝对误差为0.0111723,说明伯恩斯坦多项式不适合用来计算.这是一个理论上非常完美但实际计算效果不佳的典型例子. 实验3.2 在同一坐标系中绘出若干个勒让德多项式Pn(x)的图形. 此实验主要熟悉Mathematica的绘图功能,并增强对勒让德正交多项式Pn(x)的认识.首先以表格形式给出前5个勒让德正交多项式. TableForm\,{n,5}\]\] (结果略)为了在图形上显示Pn(x),编写一个函数ntos,将任一个整数转换为字符型整数,程序为 ntos\:=Module\IntegerDigits\;FromCharacterCode\48\]\]绘图程序 ns=2;ne=6;step=2; (初值,终值和步长,可以任意设定) st=Table\,Thickness\},{n,ns,ne,step}\]; (设置曲线颜色和线宽) txt=Table\,Text\>ntos\<>"(x)", {(n1.4)/(ne)-0.8,0.9}\]},{n,ns,ne,step}\];(文字) Plot\,{n,ns,ne,step}\]\], {x,-1,1}, PlotStyle->st,Frame->True,FrameLabel->{x,y},Epilog->txt\] (图3.9)图3.9 勒让德正交多项式P2(x),P4(x),P6(x)的图形(曲线由细到粗) 程序解释 st是一个表,用来限定曲线的颜色和粗细,它们随勒让德多项式的次数而变化,次数升高,曲线变粗.txt用于在图形上写文字,符号<>是字符连接运算符,功能是将两个字符串连接起来,所以,要用ntos函数先将整数n转换为字符n.在txt中,"P"<>ntos\<>"(x)"可用Subscript\\代替.Evaluate的作用是强制求值,请比较不用该函数的效果.另外,Epilog->txt的作用是补充绘图,即在绘完图形后写文字. 该程序只需给定勒让德正交多项式Pn(x)中n的起始值、终止值和步长,则可以绘出给定范围内的Pn(x),且曲线颜色和文字颜色一致,便于识别不同的曲线. 课堂练习 调整ns,ne和step的值,绘制奇数阶勒让德正交多项式. 实验3.3 用切比雪夫多项式对y(x)=ex,x∈\进行四次逼近. (1) 对y(x)=ex进行幂级数展开,展开到x5,记为P(x),并记P(x)舍去最高次项的多项式为Pc(x). (2) 将P(x)中x的幂次用切比雪夫多项式代替,舍去T5(x)项,并换回变量x,记所得的多项式为Tc(x),求出Tc(x). (3) 在\区间上绘出y(x)-Pc(x)和y(x)-Tc(x)的误差曲线,观察哪种逼近更精确? (4) 求y(x)的四次切比雪夫展开式和勒让德展开式,即y(x)的以切比雪夫多项式和勒让德多项式为基底的最佳平方逼近多项式S(x)和Sp(x),比较S(x),Sp(x)和Tc(x)对y(x)的逼近哪个更精确? (5) 重复以上过程,用切比雪夫多项式对y(x)=ex进行更高次(如10次)逼近. 切比雪夫多项式的一个重要性质是可以使用切比雪夫多项式降低幂级数的次数而不严重损失幂级数精度.本实验首先对切比雪夫多项式这一良好性质给予充分展示,另外,用这种方法得到的多项式与同次数的最佳逼近多项式的关系也将在实验中给予探讨. 本实验要用到n次单项式与切比雪夫多项式之间的转换,它们的关系式是: 令m=n2〗不超过n2的最大整数,则有xn=21-n∑mk=0CknTn-2k(x), n为奇数, xn=21-n∑m-1k=0CknTn-2k(x)+12CmnT0(x)〗, n为偶数. 对于以上两个函数可以如下定义,其中T\表示n次切比雪夫多项式. f\OddQ\]:=2^(1-n)Sum\T\2k,x\],{k,0,Floor\}\] f\EvenQ\]:=2^(1-n)(Sum\T\2k,x\],{k,0,Floor\-1}\]+ Binomial\\]/2 T\)以上我们使用了带条件的函数定义,OddQ为奇数问询函数,表示若传递的值为奇数,则执行该函数; EvenQ为偶数问询函数,表示若传递的值为偶数,则执行该函数.对于xn,需要采用变量替换方法将其化为切比雪夫多项式的形式,例如 {x^2/.x^2->f\,x^3/.x^3->f\} {1/2(T\+T\),1/4(3T\+T\)}这种替换方法显然缺乏通用性,Mathematica的模式替换有更强的功能,如,x^n_表示x的任意幂次,幂次为n.但这种方法不适用于x本身,若使这种方法适用于x本身,则应写为x^n_.,例如 g=4+x+x^4+(x+1)^2//Expand (结果 g=5+3x+x2+x4 ) g/.{x^n_.->f\} (x的任意n次幂替换为f\ ) 5+3T\+1/2(T\+T\)+1/8(3T\+4T\+T\) (注意输出结果中常数项并没有被替换, 若要使常数项也被替换, 则应使用下面命令) g/.{Coefficient\->Coefficient\f\, x^n_.->f\} (Coefficient\表示多项式g的常数项) 5T\+3T\+1/2(T\+T\)+1/8(3T\+4T\+T\) (替换结果)(1) 对y(x)=ex进行幂级数展开,展开到x5. p=Series\,{x,0,5}\]//Normal 1+x+x2/2+x3/6+x4/24+x5/120(2) 将P(x)中x的幂次用切比雪夫多项式代替,舍去T5(x)项,并换回变量x,所得的多项式记为Tc(x),求出Tc(x). tchy=p/.{Coefficient\->Coefficient\f\, x^n_.->f\}; (变量替换) tchy=Expand\ (表达式展开) 81/64T\+217/192T\+13/48T\+17/384T\+1/192T\+T\/1920注意tchy的最后一项为 T\/1920,由于 |T\|≤1,所以,这一项与p的最后一项x5/120相比,减小很多,舍去tchy的最后一项. tc=Drop\1\] (Drop选项中的-1表示从后向前舍去一项) 81/64T\+217/192T\+13/48T\+17/384T\+1/192T\将以上多项式变回关于x的多项式.由于T\仅是一个符号函数,需将其变换为切比雪夫多项式. tc=tc/.T->ChebyshevT;(ChebyshevT:切比雪夫多项式名字) tc=tc//Expand 1+(383x)/384+x2/2+(17x3)/96+x4/24 (注意tc与p的区别)图3.10 y(x)-Pc(x)(实线)和y(x)- Tc(x)(虚线)的误差曲线(3) 在\区间上绘制y(x)-Pc(x)和y(x)-Tc(x)的误差曲线. pc=Drop\1\] (删除p的最后一项,即去掉p中x5/120) Plot\-pc,Exp\-tc},{x,-1,1}, AxesLabel->{x,error},AxesStyle->Thick, BaseStyle->{FontSize->16},PlotRange-> All,PlotStyle->{{Red,Thick}, {Blue,Thick,Dashed}}\]结果见图3.10,其中实线为exp(x)的四次泰勒展开式的误差曲线,虚线为切比雪夫多项式降阶后得到的Tc(x)的误差曲线,可见后者的误差在区间上分布均匀,最大误差小于同阶的泰勒展开式的误差.