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
____________________________________________________________________