Chebyshev 展开

Chebyshev展开是将有限区间上的光滑函数以Chebyshev多项式为基做展开。和Taylor展开不同的是，它对展开函数的光滑性要求较低，只需连续即可。著名的Chebfun系统基础之一就是Chebyshev展开。下面是Mathematica上的一个简单的Chebyshev展开，展开系数使用Gauss-Chebyshev积分计算，积分的代数精度是2*M+1，这里M是展开的阶数。

(************************************************************************)(*                         Chebyshev approximation                      *)(************************************************************************)ChebyshevApproximation[f_, x_, a_, b_, M_] := Module[{tmp, i, j, z, F, A},(*M:{1,x,x^2,...,x^M}*)(* Comment: A[j] is accurate for polynomials of degree less than 2M+1 *)If[FreeQ[f,x],Return[f]];tmp[0] = f /. x -> a + (b - a)*(z + 1)/2;If[Length[nodes] != M + 1, (* global nodes and T to avoid reduplicate computation*)nodes = Table[Cos[(2*i + 1)*Pi/(2*(M + 1))] // GetDigit, {i, 0, M}];T = Table[ChebyshevT[i, nodes[[j]]], {i, 0, M}, {j, 1, M + 1}]];F = tmp[0] /. z -> nodes // Expand;For[j = 0, j <= M, j++,A[j] = 2/(M + 1)*F.T[[j + 1]] // GetDigit;];tmp[1] = A[0]/2 + Sum[A[j]*ChebyshevT[j, z], {j, 1, M}];tmp[2] = tmp[1] /. z -> -1 + 2*(x - a)/(b - a) (*// Expand//GetDigit*)];

这里GetDigit是自定义的一个函数，其作用是保持数字的有效位，默认是100位。

(************************************************************************)(*                             Define  GetDigit[]                       *)(************************************************************************)(* 2012.01.05, 12:21, happy!! *)Naccu=100;GetDigit[p_Plus] := Map[GetDigit, p];GetDigit[p_List] := Map[GetDigit, p];GetDigit[p_Complex]:=GetDigit[Re[p]]+I*GetDigit[Im[p]];GetDigit[c_Real]    := 0 /;Abs[c] < 10^(-Naccu+1);GetDigit[c_*f_] := GetDigit[c]*f /; NumericQ[c];GetDigit[f_] := f /; !NumericQ[f];GetDigit[c_] := SetPrecision[Apply[Rationalize[#1]*10^#2&,MantissaExponent[c]],Naccu]/;NumericQ[c]

Maple中的Chebyshev 级数展开，展开长度M由eps自适应决定，帮助文件说“The series computed is the ‘infinite’ Chebyshev series, truncated by dropping all terms with coefficients smaller than eps muliplied the largest coefficient”，函数chebyshev的源代码如下：

> with(numapprox);

>interface(verboseproc = 2);

>print(chebyshev);

proc(f::{algebraic,procedure},eqn::{name,name=algebraic..algebraic,algebraic..algebraic},eps::numeric,size::integer)option `Copyright (c) 1992 by the University of Waterloo. All rights reserved.`; local f_is_operator,x,r,epsilon,oldDigits,a,b,evalhfOK,fproc,err,nofun,c,intf,tol,maxcoef,k,K,s,y,flags,dim,g; if nargs<2 or 4<nargs then error "wrong number (or type) of arguments" elif type(f,'series') then error "expecting an algebraic expression not of type series" end if; if type(eqn,`=`) then f_is_operator:=false; x,r:=op(eqn) elif type(eqn,'name') then f_is_operator:=false; x:=eqn; r:=-1..1 else f_is_operator:=true; r:=eqn end if; if type(f,'procedure') and type(eqn,{`=`,'name'}) then f_is_operator:=true end if; if nargs<3 then epsilon:=Float(1,-Digits) elif Float(1,-Digits)<=eps then epsilon:=evalf(eps) else error "eps is less than 10^(-Digits)" end if; oldDigits:=Digits; Digits:=Digits+length(Digits+9)+1; a:=evalf(op(1,r)); b:=evalf(op(2,r)); if not type([a,b],'[numeric,numeric]') then error "non-numeric end-points of interval" elif b<a then a,b:=b,a end if; try if f_is_operator then fproc,evalhfOK:=`evalf/int/CreateProc`(f,_Operator,a,b,epsilon) else fproc,evalhfOK:=`evalf/int/CreateProc`(f,x,a,b,epsilon) end if catch "function has a pole in the interval": error "singularity in or near interval" catch: error  end try; userinfo(1,'`numapprox:-chebyshev`',printf("procedure for evaluation is:\\n%a\\n",eval(fproc))); if nargs=4 then dim:=size else dim:=487 end if; flags:=array(1..3); Assert(not assigned(intf)); if evalhfOK and Digits<=evalhf(Digits) then try c:=hfarray(1..dim); intf:=evalhf(`evalf/int/ccquad`(fproc,a,b,epsilon,dim,var(c),var(flags),true)) catch: userinfo(1,'`numapprox:-chebyshev`',nprintf("evalhf mode unsuccessful - try using software floats")) end tryend if; if not assigned(intf) then try c:=array(1..dim); intf:=`evalf/int/ccquad`(fproc,a,b,epsilon,dim,c,flags,false) catch: error  end try end if; if type(intf,{'infinity','undefined'}) then error "singularity in or near interval" end if; Digits:=oldDigits; err:=flags[1]; nofun:=round(flags[2]); if flags[3]<>0 or max(epsilon*abs(intf),0.001*epsilon)+Float(1,-Digits)<err then if dim<=487 then userinfo(2,'`numapprox:-chebyshev`',nprintf("try again using more function evaluations")); s:=procname(fproc,convert(a,'rational')..convert(b,'rational'),epsilon,4375); return if(type(eqn,'range'),eval(s),s(x))else error "singularity in or near interval" end if end if; c:=map((v,d)->evalf(v/d),[seq(c[k],k=1..nofun)],nofun - 1); maxcoef:=max(1,op(c)); tol:=1/5*epsilon*maxcoef; for K from nofun by -1 to 1 while abs(c[K])<tol do  end do; for k by 2 to K while abs(c[k])<tol do  end do; if K<k and 0<K then if f_is_operator then error "even coefficients are zero - try f(x)/x" elif a+1.0<>0. or b - 1.0<>0. then userinfo(2,'`numapprox:-chebyshev`',nprintf("even coefficents are zero -- transform to interval -1..1")); a:=op(1,r); b:=op(2,r); g:=eval(f,x=1/2*(b - a)*y+1/2*a+1/2*b); s:=procname(g,y=-1..1,epsilon,dim); s:=subs(y=2*x/(b - a) - (a+b)/(b - a),s) else userinfo(2,'`numapprox:-chebyshev`',nprintf("function is odd -- try f(x)/x")); Digits:=Digits+2; s:=procname(f/x,x=r,1/100*epsilon); s:=numapprox:-chebmult('T'(1,x),s); Digits:=Digits - 2; if type(s,`+`) then s:=map(proc(t,tol) if abs(lcoeff(t))<tol then 0 else t end if end proc,s,tol) elif abs(lcoeff(s))<tol then s:=0 end if end if; s:=evalf(s); s:=numapprox:-chebsort(s); return send if; c:=[seq(if(abs(c[k])<tol,0.,c[k]),k=1..max(K,1))]; a:=op(1,r); b:=op(2,r); y:=2*x/(b - a) - (a+b)/(b - a); s:=1/2*c[1]*'T'(0,y)+add(c[k]*'T'(k - 1,y),k=2..K); s:=numapprox:-chebsort(s); if type(eqn,'range') then unapply(s,x) else s end if end proc