Fitting Chebyshev to 2D data

Question

I have complex data points on a 2D square (say $[0,10]^2$). I tried the following example function: $xye^{i(x+y)}$.

I split up the Re and Im parts and just fitted both separately with NonlinearModelFit:

m=15;
f[x_, y_] = x*y*Exp[I  (x + y)];
dataR = Flatten[
   Table[{x, y, Re@f[x, y]}, {x, 0, 10, 0.1}, {y, 0, 10, 0.1}], 1];
dataI = Flatten[
   Table[{x, y, Im@f[x, y]}, {x, 0, 10, 0.1}, {y, 0, 10, 0.1}], 1];

model[m_, x_, y_] := 
  Sum[a[n]  ChebyshevT[n, 1/5  x - 1], {n, 0, m}] Sum[
    b[n]  ChebyshevT[n, 1/5  y - 1], {n, 0, m}];

nlmR[m_Integer?Positive, x_, y_] := 
 NonlinearModelFit[dataR, model[m, x, y], 
  Join[Array[a, m + 1, 0], Array[b, m + 1, 0]], {x, y}, 
  Method -> "Gradient"]
Plot3D[Evaluate@{Abs[Re[f[x, y]] - nlmR[m, x, y] // Normal], 
   Abs@Re[f[x, y]]}, {x, 0, 10}, {y, 0, 10}, PlotRange -> Full]

nlmI[m_Integer?Positive, x_, y_] := 
 NonlinearModelFit[dataI, model[m, x, y], 
  Join[Array[a, m + 1, 0], Array[b, m + 1, 0]], {x, y}, 
  Method -> "Gradient"]
Plot3D[Evaluate@{Abs[Im[f[x, y]] - nlmI[m, x, y] // Normal], 
   Abs[Im[f[x, y]]]}, {x, 0, 10}, {y, 0, 10}, PlotRange -> Full]

and I plotted their absolute difference and the original function. The error I make is of order 1. (Increasing m makes it worse)

Is there a way to fit such data with Chebyshevs in 2D? They should be a complete basis set of polynomials.

(P.s. In 1D, this seems to work decently well.)

Answer 1

This is not MMA fault. The problem you have is evident from (for simplicity I only consider the real part):

Plot3D[Re[f[x, y]], {x, 0, 10}, {y, 0, 10}, PlotRange -> Full]

As you see, the function can not be written by a fx[x] fy[y]. That is, it is not separable into a product of 2 functions.

However, the basis you are using is exactly doing this. You must include the cross terms, like the following, where we plot the error and the fit.:

m = 15;
f[x_, y_] = x*y*Exp[I   (x + y)];
dataR = Flatten[
   Table[{x, y, Re@f[x, y]}, {x, 0, 10, 0.1}, {y, 0, 10, 0.1}], 1];

model[x_, y_] = 
  Flatten[Sum[
    a[n1, n2]   ChebyshevT[n1, 1/5   x - 1]  ChebyshevT[n2, 
      1/5   y - 1], {n1, 0, m}, {n2, 0, m}]] ;
param := Flatten[Table[a[n1, n2] , {n1, 0, m}, {n2, 0, m}]] ;

fit[x_, y_] = 
 NonlinearModelFit[dataR, model[x, y], param, {x, y}][x, y]

Plot3D[Evaluate@{Re[f[x, y]] - fit[x, y]}, {x, 0, 10}, {y, 0, 10}, 
 PlotRange -> Full]
Plot3D[Evaluate@{fit[x, y]}, {x, 0, 10}, {y, 0, 10}, PlotRange -> Full]