It is my understanding that you are seeking a method to interpolate a 2D polynomial over a grid of Chebyshev points, extending the concept of 1D polynomial fitting in MATLAB to two dimensions for a matrix M, and require guidance on how to do this for a high order N, such as 48. I assume that you are using MATLAB R2023b.
To interpolate a 2D polynomial given a matrix of values M at Chebyshev points in both the x and y directions, you can use a tensor product of Chebyshev polynomials. In 1D, you used polyfit to find the coefficients of the polynomial that fits your data. In 2D, the process is more involved because you need to fit a polynomial in two variables.
Here is a MATLAB function that uses the Chebyshev polynomial basis to find the coefficients for a 2D polynomial interpolation:
function P2D = chebyshevInterpolation2D(M, N)
% M is the matrix of function values at Chebyshev points
% N is the order of the Chebyshev polynomial
% Get the Chebyshev points
x = cos((0:N) * pi / N);
y = x;
% Create Chebyshev polynomials at these points
T = zeros(N+1, N+1);
for k = 0:N
T(:, k+1) = cos(k * acos(x));
End
% Calculate the coefficients for the 2D polynomial
% We use the orthogonality of Chebyshev polynomials
% and compute the coefficients using a discrete cosine transform (DCT)
a = dct2(M) / (N+1);
a(1,:) = a(1,:) / 2;
a(:,1) = a(:,1) / 2;
% Construct the 2D polynomial using the tensor product
P2D = zeros(size(M));
for i = 0:N
for j = 0:N
P2D = P2D + a(i+1, j+1) * T(:,i+1) * T(:,j+1)';
end
end
end
To use this function with your example data, we would call it like this:
N = 8;
% define x and y as N+1 Chebyshev points
x = cos((0:N) * pi / N);
y = x;
[X, Y] = meshgrid(x, y);
% define matrix M as dummy
M = sin(X * pi) + cos(pi * Y.^2);
% Interpolate a 2D polynomial
P2D = chebyshevInterpolation2D(M, N);
% Now P2D contains the interpolated 2D polynomial values
Note that as N increases, the complexity of the polynomial increases, and the calculation becomes more computationally intensive.