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 try
end 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 s
end 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
