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 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