Re: 2-D Chebyshev Polynomial Regression

*To*: mathgroup at smc.vnet.net*Subject*: [mg13321] Re: 2-D Chebyshev Polynomial Regression*From*: Paul Abbott <paul at physics.uwa.edu.au>*Date*: Mon, 20 Jul 1998 02:49:42 -0400*Organization*: University of Western Australia*References*: <6nske5$fc8@smc.vnet.net>*Sender*: owner-wri-mathgroup at wolfram.com

Chris Farr wrote: > Has anyone created a function in Mathematica to approximate a function > of two variables using 2-D Chebyshev Polynomial Regression? > > That is, has someone created a Mathematica algorithm which takes as its > input a real valued function f(x,y) defined on [a,b] X [c,d] and > returns a Chebyshev polynomial approximation p(x,y)? The following code may do what you want. It takes a real valued function f[x,y] defined on [-1,1] X [-1,1] and returns a Chebyshev polynomial approximation p[f][x,y]. Extending the one-dimensional Chebyshev formula to two dimensions in an obvious fashion, In[1]:= c[i_, j_, m_, n_][f_] := c[i, j, m, n][f] = Chop[N[4/(n m) Sum[f[Cos[Pi (k-1/2)/m], Cos[Pi (l-1/2)/n]]* Cos[Pi i (k-1/2)/m] Cos[Pi j (l-1/2)/n],{k,1,m},{l,1,n}]]] and defining a short-hand notation for the Chebyshev T polynomials (with the correction for the zeroth coefficient incorporated into the polynomial itself to simplify the notation), In[2]:= t[n_][x_]=ChebyshevT[n,x]; In[3]:= t[0][x_] = 1/2; then the (m,n)-th order polynomial approximation for f reads In[4]:= p[m_,n_][f_]:= p[m,n][f] = Function[{x,y},Sum[c[k,l,m,n][f] t[k][x] t[l][y],{k,0,m-1},{l,0,n-1}]] For example, with the (pure) function In[5]:= f = Function[{x,y},(x y^2)/(x + y + 1)]; the (3,3)-approximation reads In[6]:= p[3,3][f][x,y] Out[6]= 2 2 0.375 y x + 0.125 (2 y - 1) x + 0.125 x + 0.125 (2 x - 1) - 2 2 2 0.375 (2 x - 1) y - 0.375 y + 0.125 (2 x - 1) (2 y - 1) + 2 0.125 (2 y - 1) + 0.125 This code uses dynamic programming (i.e., g[x_]:= g[x] = ...) to store intermediate computations and pure functions (so that computation for new functions f is immediate). The transformation from [a,b] X [c,d] to [-1,1] X [-1,1] is x -> (2x-a-b)/(b-a) y -> (2y-c-d)/(d-c) with inverse x -> ((b-a)x+a+b)/2 y -> ((d-c)y+c+d)/2 so only a simple modification of the above code is required to generalize it to [a,b] X [c,d]. Note also that Mathematica includes a number of packages for related computations including NumericalMath`Approximations` and Calculus`Pade`. Hope that this helps. Cheers, Paul ____________________________________________________________________ Paul Abbott Phone: +61-8-9380-2734 Department of Physics Fax: +61-8-9380-1014 The University of Western Australia Nedlands WA 6907 mailto:paul at physics.uwa.edu.au AUSTRALIA http://www.pd.uwa.edu.au/~paul God IS a weakly left-handed dice player ____________________________________________________________________