2d Chebyshev polynomial interpolation

Hello Community,

may i ask you for your help.

I have a full matrix M with the N+1 number of elments in each direction.

The x and y coordinates are the roots of the chebyshev polynomial (of the first kind).

My goal is to find an exact polynomail in x and y

I have achieved this in 1D with the following line:

[P, S, mu]=polyfit(x,y,N);

Here is an example of the data:

N = 8;
% define x and y as N+1 chebychev 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);

How could I interpolate a 2D polynom in x and y?

If N = 1 the size of M is [2,2]. And I could youe a bilinear interpolation polynom:

a11 = f(1,1);
a12 = f(1,2)-f(1,1);
a21 = f(2,1)-f(1,1);
a22 = f(2,2)-(f(2,1)+f(1,2));
fBilinear = a11 +a12*x +a21*y +a22*x*y; 

How could i extend this formula to N = 48 ?

Any help would be greatly appreciated.

Answer 1

[Kautuk Raj]

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.