Chebyshev 展开

2023-12-28 23:53| 来源: 网络整理| 查看: 265

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

【本文地址】

公司简介

联系我们