scatteredInterpolant: what is linear interpolation in 2d?

%coords
supportPts = [3 3; 3.3 3; 3 3.25; 3.3 3.25; 3.6 3; 3.6 3.25; 3 3.5; 3.3 3.5; 3.6 3.5];
%values
Fval = [0 0.1121 0.064604 0.184942 0.233029 0.352622 0.121444 0.206496 0.375685]';
%interpolation object
Interp = scatteredInterpolant(supportPts(:,1),supportPts(:,2),Fval);

%evaluate at center of bottom left element
Interp(3.15, 3.125)

ans = 0.0884

See the above example with nine points that represent four axis-parrallel elements.

I would have expected that the value of the interpoland at the center of the bottom left element is the mean value of the coresponding four corner values, i.e.,

%expected value at (3.15,3.125)
0.25*(0+0.064604+0.184942+0.1121)

ans = 0.0904

but this is clearly not so.

So what interpolation scheme is scatteredInterpoland using (linear) to calculate values in the interior of the elements? At the element boundaries, I double-checked that there is just 1d interpolation involved.

Answer 1

[Matt J]

No, the support points are being divided into triangular elements, not rectangles. It is puzzling that you would be using scatteredInterpolant instead of griddedInterpolant for data like this.

[SA-W]

Thank you! :)

Well, I would like to use griddedInterpolant, but it throws an error

Interpolation requires at least two sample points for each grid dimension.

How can I transform my 7x,y,z} pairs from my question into the structure for griddedInterpolant?

[Matt J]

supportPts = [3 3; 3.3 3; 3 3.25; 3.3 3.25; 3.6 3; 3.6 3.25; 3 3.5; 3.3 3.5; 3.6 3.5];
%values
Fval = [0 0.1121 0.064604 0.184942 0.233029 0.352622 0.121444 0.206496 0.375685]';

[~,is]=sortrows(supportPts); Fval=reshape(Fval(is),3,[]);

Interp=griddedInterpolant({[3,3.3,3.6], [3,3.25,3.5]}, Fval);

Interp(3.15, 3.125)

ans = 0.0904

[SA-W]

Is there a solution without hard-coding the line

Interp=griddedInterpolant({[3,3.3,3.6], [3,3.25,3.5]}, Fval);

?

[Matt J]

You can put whatever you want there that generates the x and y grid coordinates.

[SA-W]

Say, supportPts is my input (nx2) matrix which contains n pairs of x and y coordinates in arbitrary order that define a rectangular grid. The goal should be to derive a griddedInterpolant without hard-coding it for a specific n.

Would

Interp=griddedInterpolant(unique(sort(supportPts(:,1)),  unique(sort(supportPts(:,2))}, Fval);

work always?

[Bruno Luong]

"Say, supportPts is my input (nx2) matrix which contains n pairs of x and y coordinates in arbitrary order that define a rectangular grid. The goal should be to derive a griddedInterpolant without hard-coding it for a specific n."

Use sortrows command to order you data in grid, then call griddedInterpolant .

[Matt J]

Say, supportPts is my input (nx2) matrix which contains n pairs of x and y coordinates in arbitrary order that define a rectangular grid.

As Bruno says, you would sort them using sortrows() and then unique().

[P,is]=sortrows(supportPts);

x=unique(P(:,1));

y=unique(P(:,2));

Fval=reshape(Fval(is),numel(x),numel(y));

Interp=griddedInterpolant({x,y},Fval)

But why would you even be starting with data in such a form? If it's a rectangular grid, you would usually be given the sampling of x and y in separate vectors from the beginning, which is the input format that meshgrid, ndgrid, and various other grid operation commands require.

[X,Y]=ndgrid(3:0.3:3.6  , 3:0.25:3.5), Z=X.^2+Y.^2

X = 3×3

3.0000 3.0000 3.0000 3.3000 3.3000 3.3000 3.6000 3.6000 3.6000

Y = 3×3

3.0000 3.2500 3.5000 3.0000 3.2500 3.5000 3.0000 3.2500 3.5000

Z = 3×3

18.0000 19.5625 21.2500 19.8900 21.4525 23.1400 21.9600 23.5225 25.2100

[SA-W]

But why would you even be starting with data in such a form? If it's a rectangular grid, you would usually be given the sampling of x and y in separate vectors from the beginning, which is the input format that meshgrid, ndgrid, and various other grid operation commands require.

Well, I come from a finite-element-like context. There, the degrees of freedom (the Fval's) and the corresponding points where they are defined (the supportPts) are sorted to reduce the bandwidth of a linear system. That's why I do not start with vectors x and y as you suggested but have pairs of coordinates in a random order.

GriddedInterpolant does not divide the elements into triangles but rather performs a bilinear interpolation?

[Matt J]

Yes, all griddedInterpolant's interpolaiton methods operate on rectangular data neighborhoods and the default is bilienar interpolation.

[SA-W]

and the default is bilienar interpolation.

I expected a different output from griddedInterpolant when comparing the results to those of my visualization tool (paraview).

In my finite-element framework, I have bilinear interpolation functions

N_0 = 0.25*(1-u_1)*(1-u_2)

N_1 = 0.25*(1+u_1)*(1-u_2)

N_2 = 0.25*(1+u_1)*(1+u_2)

N_3 = 0.25*(1-u_i)*(1+u_2)

where u_1,u_2 \in [-1,1] are reference coordinates. If you think of a rectangle with corners [-1,-1], [1,-1], [1,1], [-1,1] , then the N_i are interpolation functions. For instance, N_0 is one at [-1,-1] and zero at the remaining three vertices. So if we have a point [x,y] of a rectangular element, this point is mapped into the reference coordinates [u_1,u_2] and the interpolated value is obtained as

val = N_i * E_i (i=1,...,4),

where the E_i are the local vertex values of the element where the point [x,y] lives.

This is clearly a bilinear interpolation and I expected this to be the same what griddedInterpolant does with linear method. However, if I evaluate the interpoland from my question at the point [3.2,3.2], I obtain 0.13714 in matlab and 0.13081 in paraview. At the element centers, the output is the same, also along the element "boundaries", only in the interior of the elements there is a mismatch.

Given that, is griddedInterpolant doing something different than what I described?

[Bruno Luong]

"Given that, is griddedInterpolant doing something different than what I described?"

No. Your test is wrong.

This is the shape of interpolation fuction with linear method

x=[0 1];
y=[0 1];
xq=linspace(0,1,41);
yq=linspace(0,1,41);
[X,Y]=ndgrid(x,y);
[XQ,YQ]=ndgrid(xq,yq);
Z=rand(2)

Z = 2×2

0.6745 0.8370 0.1259 0.8627

f=griddedInterpolant(X,Y,Z,'linear');
% any black line (// to x/y axis) is a staight line
% The 4 corner values are (2x2) Z avlues
surf(xq,yq,f(XQ,YQ))

format long

This is the value of mlidle point

f(0.5,0.5)

ans = 
   0.625023461369597

and the mean of Z

mean(Z,"all")

ans = 
   0.625023461369597

It does match.

[Matt J]

If you think of a rectangle with corners [-1,-1], [1,-1], [1,1], [-1,1] , then the N_i are interpolation functions. For instance, N_0 is one at [-1,-1] and zero at the remaining three vertices.

If so, then your triangle basis functions should fall off with slope 1/2, whereas in your formulae they appear to be falling off with slope 1.

[Bruno Luong]

Here is the 4 bilinear functions that has 0s at three corners and 1 in one of them

x=[-1 1];
y=[-1 1];
xq=linspace(-1,1,41);
yq=linspace(-1,1,41);
[X,Y]=ndgrid(x,y);
[XQ,YQ]=ndgrid(xq,yq);
figure
for k=1:4
    Z = zeros(2,2);
    Z(k) = 1;
    fk=griddedInterpolant(X,Y,Z);
    subplot(2,2,k)
    surf(xq,yq,fk(XQ,YQ))
end 

[SA-W]

@Bruno Luong

I think my basis functions

N_0 = 0.25*(1-u_1)*(1-u_2)

N_1 = 0.25*(1+u_1)*(1-u_2)

N_2 = 0.25*(1+u_1)*(1+u_2)

N_3 = 0.25*(1-u_i)*(1+u_2)

are correct. If I plot them individually over

x=[-1 1];
y=[-1 1];
[X,Y]=ndgrid(x,y);

I obtain the four plots you showed me.

It does match.

Yes. As I said, the values match at the center of the element and on its boundary. But they do not match at other points than the center in the interior of the element.

[SA-W]

@Matt J

If so, then your triangle basis functions should fall off with slope 1/2, whereas in your formulae they appear to be falling off with slope 1.

The basis functions are rectangular. So I think the factor 0.25 is correct.

[Bruno Luong]

"I obtain the four plots you showed me. "

It does match.

"Yes. As I said, the values match at the center of the element and on its boundary. But they do not match at other points than the center in the interior of the element. "

If it matches the functions (basis), it must match everywhere inside the rectangle because (MATLAB) bilinear function just the sum :

f(x,y) = sum Z(x_i,y_j) * B_i,j(x,y), i=1,2, j=1,2

where B_i,j are the 4 basis functions (that match, and you denote by N_0, ... N_3) and Z(x_i,y_j) are value at the 4 corners (given).

You make a wrong assessment somewhere.

[Bruno Luong]

"The basis functions are rectangular. So I think the factor 0.25 is correct. "

This I agree since the size of your rectangle is 2x2, 0.25 is just 1/(2*2), and your formula is

N_0 = w1*w2

N_1 = (w1)*(1-w2)

N_2 = (1-w1)*(1-w2)

N_3 = w1*(1-w2)

with w1 := (1-u1)/2, w2 := (1-u2)/2. Both w1 and w2 are in (0,1) inside the rectangle.

This is bilinear formula.

[SA-W]

@Bruno So are the N_i that I provided wrong? They sum up to one everywhere on the reference rectangle, which is why I think they are correct. *f(x,y) = sum Z(x_i,y_j) * B_i,j(x,y), i=1,2, j=1,2 where B_i,j are the 4 basis functions (that match) and Z(x_i,y_j) are value at the 4 corners (given).* What are the analytical formulae for the B_i,j here?

[Matt J]

The basis functions are rectangular.

Bilinear interpolation basis functions have a rectangular footprint, but triangular profiles. I think this is what you want:

f=@(x) (1-abs(x)) .*(abs(x)<=1); %half-width=1
f=@(x) f(x/2));  %half-width=2

N_0 = f(u_1-1)*f(u_2-1)
N_1 = f(u_1-1)*f(u_2+1)
N_2 = f(u_1+1)*f(u_2+1)
N_3 = f(u_1+1)*f(u_2-1)

[Bruno Luong]

What are the analytical formulae for the B_i,j here?

You already give the formula inside the rectangle (-1,1) x (-1,1)

N_0 = 0.25*(1-u_1)*(1-u_2)
N_1 = 0.25*(1+u_1)*(1-u_2)
N_2 = 0.25*(1+u_1)*(1+u_2)
N_3 = 0.25*(1-u_1)*(1+u_2)

N_0 correspond to Z(-1,-1)

N_1 corresponds to Z(+1,-1)

N_2 corresponds to Z(+1,+1)

N_3 corresponds to Z(-1,+1)

Put together the formula for (x,y) in (-1,1)x(-1,1) is

f(x,y)  = 0.25*(Z(-1,-1)*(1-x)*(1-y) + ...
                Z(+1,-1)*(1+x)*(1-y) + ...
                Z(+1,+1)*(1+x)*(1+y) + ...
                Z(-1,+1)*(1-x)*(1+y));

If you find MATLAB does not return this let us know.

[Bruno Luong]

You can read this thread where I post formula (using matrix form) of bilinear interpolation

[Bruno Luong]

Code to illustrate the bilinear formula in (-1,1) x (-1,1)

x=[-1 1];
y=[-1 1];

xq=linspace(-1,1,41);
yq=linspace(-1,1,31);

[X,Y]=ndgrid(x,y);
[XQ,YQ]=ndgrid(xq,yq);

Z=rand(2);
f=griddedInterpolant(X,Y,Z,'linear');

FXY  = 0.25*(Z(1,1)*(1-XQ).*(1-YQ) + ...
             Z(2,1)*(1+XQ).*(1-YQ) + ...
             Z(2,2)*(1+XQ).*(1+YQ) + ...
             Z(1,2)*(1-XQ).*(1+YQ));


figure
subplot(2,1,1);
surf(yq,xq,f(XQ,YQ));
title('griddedInterpolant')
subplot(2,1,2);
surf(yq,xq,FXY);
title('formula')

[SA-W]

@Bruno Luong @Matt J

If you find MATLAB does not return this let us know.

I managed to create an example that illustrates the differences I see between griddedInterpolant and the formulae (sum of basis functions times nodal values):

%coordinates
supportPts = [3 3; 3.3 3; 3 3.25; 3.3 3.25; 3.6 3; 3.6 3.25; 3 3.5; 3.3 3.5; 3.6 3.5];
%nodal values
Fval = [0 0.1121 0.064604 0.184942 0.233029 0.352622 0.121444 0.206496 0.375685]';
%gridded interpolation object
[P,is] = sortrows(supportPts); 
Inv1 = unique(P(:,1));
Inv2 = unique(P(:,2));
Fval = reshape(Fval(is),numel(Inv1),numel(Inv2));
Interp = griddedInterpolant({Inv1,Inv2},Fval);
%vertices and corresponding values of bottom left element
vertices = [3 3; 3.3 3; 3.3 3.25; 3 3.25];
vertexValues = [0 0.1121 0.184942 0.064604];
%point where interpoland should be evaluated (currently center)
ptReal = [3.15 3.125];
%determine the coordinates on the reference element (-1,1)x(-1,1)
fun = @(x)myFun(x, ptReal, vertices);
ptRef = fsolve(fun, [1 1]);

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

disp('value formulae with basis functions: ');

value formulae with basis functions: 

fEval(ptRef, vertexValues)

ans = 0.0904

disp('value griddedInterpolant: ');

value griddedInterpolant: 

Interp(ptReal(1), ptReal(2))

ans = 0.0904

%function to determine reference coordinates of ptEval
function F = myFun(x, ptEval, vertexValues)
N_0 = 0.25*(1-x(1))*(1-x(2));
N_1 = 0.25*(1+x(1))*(1-x(2));
N_2 = 0.25*(1+x(1))*(1+x(2));
N_3 = 0.25*(1-x(1))*(1+x(2));
F = N_0.*vertexValues(1,:) + N_1.*vertexValues(2,:) + ...
    N_2.*vertexValues(3,:) + N_3.*vertexValues(4,:) - ptEval;
end
%function to interpolate using local basis functions
function fval = fEval(x, values)
N_0 = 0.25*(1-x(1))*(1-x(2));
N_1 = 0.25*(1+x(1))*(1-x(2));
N_2 = 0.25*(1+x(1))*(1+x(2));
N_3 = 0.25*(1-x(1))*(1+x(2));
N = [N_0 N_1 N_2 N_3];
fval = N*values';
end

As you can see, both methods return the same value at the element center. But if you define

ptReal = [3.2 3.2];

we get

%value formulae with basis functions: 
ans =
0.056050010755214
%value griddedInterpolant: 
ans =
0.032302000000000

0.05605 is also what I see in my visualization tool. Note that et the center (ptReal = [3.15 3.125]), both methods give the same value.

So where is the "error" in my calculation?

[Bruno Luong]

"So where is the "error" in my calculation? "

Your basis functions are wrong.

N_0 = 0.25*(1-x(1))*(1-x(2));
N_1 = 0.25*(1+x(1))*(1-x(2));
N_2 = 0.25*(1+x(1))*(1+x(2));
N_3 = 0.25*(1-x(1))*(1+x(2));

If you know formula for generic grid points (not, (-1, 1)) see the the thread in the link I provide.

[Stephen23]

You have accidentally transposed your data. This is because SORTROWS does not return the data in NDGRID format as GRIDDEDINTERPOLANT requires when using grid vectors as you did:

https://www.mathworks.com/help/matlab/ref/griddedinterpolant.html#bvh2f0p-5

First lets have a look at an NDGRID example, to see what order it has:

[X,Y] = ndgrid([3,3.3,3.6],[3,3.25,3.5,Inf])

X = 3×4

3.0000 3.0000 3.0000 3.0000 3.3000 3.3000 3.3000 3.3000 3.6000 3.6000 3.6000 3.6000

Y = 3×4

3.0000 3.2500 3.5000 Inf 3.0000 3.2500 3.5000 Inf 3.0000 3.2500 3.5000 Inf

Note the linear order of elements in X is: 3, 3.3, 3.6, 3, 3.3... etc. Now lets try SORTROWS:

supportPts = [3,3;3.3,3;3,3.25;3.3,3.25;3.6,3;3.6,3.25;3,3.5;3.3,3.5;3.6,3.5];
Fval = [0;0.1121;0.064604;0.184942;0.233029;0.352622;0.121444;0.206496;0.375685];
[P,is] = sortrows(supportPts) % is P in NDGRID order? (hint: no)

P = 9×2

3.0000 3.0000 3.0000 3.2500 3.0000 3.5000 3.3000 3.0000 3.3000 3.2500 3.3000 3.5000 3.6000 3.0000 3.6000 3.2500 3.6000 3.5000

is = 9×1

1 3 7 2 4 8 5 6 9

From this point on all results are rubbish, because P is clearly not in NDGRID order and therefore neither will FVAL be when you use those indices. You have effectively swapped the dim1 and dim2, transposing your data. This is easy to see when you plot the data (which you should always do!):

scatter3(supportPts(:,1),supportPts(:,2),Fval(:),99,'g','filled')
hold on
Inv1 = unique(P(:,1));
Inv2 = unique(P(:,2));
Fone = reshape(Fval(is),numel(Inv1),numel(Inv2))

Fone = 3×3

0 0.1121 0.2330 0.0646 0.1849 0.3526 0.1214 0.2065 0.3757

[M1,M2] = ndgrid(Inv1,Inv2) % expected order of GRIDDEDINTERPOLANT

M1 = 3×3

3.0000 3.0000 3.0000 3.3000 3.3000 3.3000 3.6000 3.6000 3.6000

M2 = 3×3

3.0000 3.2500 3.5000 3.0000 3.2500 3.5000 3.0000 3.2500 3.5000

scatter3(M1(:),M2(:),Fone(:),42,'r','filled') % Ooops, swapped dimensions!
hold off

You can see that the data are reflected X-Y. A simple solution is to specify the SORTROWS column order:

[P,is] = sortrows(supportPts,[2,1]);
Inv1 = unique(P(:,1));
Inv2 = unique(P(:,2));
Ftwo = reshape(Fval(is),numel(Inv1),numel(Inv2));
Interp = griddedInterpolant({Inv1,Inv2},Ftwo);

Lets plot those data:

scatter3(supportPts(:,1),supportPts(:,2),Fval(:),99,'g','filled')
hold on
scatter3(M1(:),M2(:),Ftwo(:),42,'r','filled') % Aaaah, much better!
hold off

Now try your various interpolation points:

vertices = [3 3; 3.3 3; 3.3 3.25; 3 3.25];
vertexValues = [0 0.1121 0.184942 0.064604];
pt1 = [3.15,3.125];
fun = @(x)myFun(x, pt1, vertices);
ptRef = fsolve(fun, [1 1]);

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

fEval(ptRef, vertexValues)

ans = 0.0904

val1 = Interp(pt1(1), pt1(2))

val1 = 0.0904

pt2 = [3.2,3.2];
fun = @(x)myFun(x, pt2, vertices);
ptRef = fsolve(fun, [1 1]);

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

fEval(ptRef, vertexValues)

ans = 0.1308

val2 = Interp(pt2(1), pt2(2))

val2 = 0.1308

Alternatively you could avoid those grid vectors entirely:

obj2 = griddedInterpolant(...
    reshape(P(:,1),3,3),... or M1 & M2
    reshape(P(:,2),3,3),Ftwo);
obj2(pt1(1),pt1(2))

ans = 0.0904

obj2(pt2(1),pt2(2))

ans = 0.1308

And finally checking against another routine:

interpn(M1,M2,Ftwo,pt1(1),pt1(2))

ans = 0.0904

interpn(M1,M2,Ftwo,pt2(1),pt2(2))

ans = 0.1308

"0.05605 is also what I see in my visualization tool."

Note that the value for [3.2,3.2] cannot be less than that for [3.15,3.125] because you move closer to larger values along both axes, so stating that the expected value is smaller simply makes no sense. Rather than blindly copying results from a tool, you should always perform some basic sanity checks yourself.

Lets plot them. Try rotating and see how they line up (unfortunately this is not (yet) interactive on this forum). I also added a contour to show where all solutions == 0.05605 lie, i.e. nowhere near [3.2,3.2].

mesh(M1,M2,Ftwo,'EdgeColor','k','FaceColor','none')
hold on
contour3(M1,M2,Ftwo,[0.05605,0.05605],'-r') % == 0.05605
scatter3(supportPts(:,1),supportPts(:,2),Fval(:),99,'g','filled')
scatter3(pt1(1),pt1(2),val1,64,'b','filled')
scatter3(pt2(1),pt2(2),val2,64,'k','filled')

So far everything seems to be working as expected. Unmodified support functions:

function F = myFun(x, ptEval, vertexValues)
N_0 = 0.25*(1-x(1))*(1-x(2));
N_1 = 0.25*(1+x(1))*(1-x(2));
N_2 = 0.25*(1+x(1))*(1+x(2));
N_3 = 0.25*(1-x(1))*(1+x(2));
F = N_0.*vertexValues(1,:) + N_1.*vertexValues(2,:) + ...
    N_2.*vertexValues(3,:) + N_3.*vertexValues(4,:) - ptEval;
end
%function to interpolate using local basis functions
function fval = fEval(x, values)
N_0 = 0.25*(1-x(1))*(1-x(2));
N_1 = 0.25*(1+x(1))*(1-x(2));
N_2 = 0.25*(1+x(1))*(1+x(2));
N_3 = 0.25*(1-x(1))*(1+x(2));
N = [N_0 N_1 N_2 N_3];
fval = N*values';
end

Lesson: always always always look at your data.

[Bruno Luong]

Here is the code to check the bilinear formula. No need for fsolve.

vertices = [3 3; 3.3 3; 3.3 3.25; 3 3.25];
vertexValues = [0 0.1121 0.184942 0.064604];
pt1 = [3.15,3.125];
pt2 = [3.2,3.2];
[XY,is] = sortrows(vertices, [2 1]);
x = unique(XY(:,1));
y = unique(XY(:,2));
nx = length(x);
ny = length(y);
sz = [nx ny];
X = reshape(XY(:,1),sz);
Y = reshape(XY(:,2),sz);
Z = reshape(vertexValues(is),sz);
f=griddedInterpolant(X,Y,Z);
format long
f(pt1(1),pt1(2))

ans = 
   0.090411500000000

bilinearformula(x, y, Z, pt1(1), pt1(2))

ans = 
   0.090411500000000

f(pt2(1),pt2(2))

ans = 
   0.130810133333333

bilinearformula(x, y, Z, pt2(1), pt2(2))

ans = 
   0.130810133333333

%% No check for overflowed so the formula is clearer
function zq = bilinearformula(x, y, Z, xq, yq)
i = discretize(xq,x);
j = discretize(yq,y);
dx = x(i+1)-x(i);
dy = y(j+1)-y(j);
wx2 = (xq-x(i))./dx;
wx1 = 1-wx2;
wy2 = (yq-y(j))./dy;
wy1 = 1-wy2;
B = [wx1.*wy1, wx1.*wy2, wx2.*wy2, wx2.*wy1]; % basis functions evaluate at (xq,yq)
Zij = [Z(i,j); Z(i,j+1); Z(i+1,j+1); Z(i+1,j)];
zq = B * Zij;
end

[SA-W]

@Stephen23

Thank you, transposing the data solved the issue. I compared the output between griddedInterpolant and paraview at several points and they are identical as it should be the case.

Note that the value for [3.2,3.2] cannot be less than that for [3.15,3.125] because you move closer to larger values along both axes, so stating that the expected value is smaller simply makes no sense. Rather than blindly copying results from a tool, you should always perform some basic sanity checks yourself.

Yes, I mixed some coordinates. The value 0.05605 I referred to corresponds to [3.15 3], which also stands out in your graphics.

[SA-W]

@Bruno Luong

Here is the code to check the bilinear formula. No need for fsolve.

Yes, but I wanted to reproduce the workflow of my extermal program, which is to only work on the reference elements and map a point real coordinates into reference coordinates.

You wrote,

Your basis functions are wrong.

N_0 = 0.25*(1-x(1))*(1-x(2));
N_1 = 0.25*(1+x(1))*(1-x(2));
N_2 = 0.25*(1+x(1))*(1+x(2));
N_3 = 0.25*(1-x(1))*(1+x(2));

If you know formula for generic grid points (not, (-1, 1)) see the the thread in the link I provide.

I think the basis functions are not wrong, but finding the reference coordinates of a point in real space and working just with the basis functions on (-1,1)x(-1,1) is an alternative to your bilinearformula(x, y, Z, xq, yq), right?

[Torsten]

Why do you ask for better solutions if you are satisfied with your coding ?

The logical way to proceed when doing bilinear interpolation is like in @Bruno Luong 's function "bilinearformula". It chooses the left lower point of the rectangle as reference point for interpolation while your scheme seems to choose the center of the rectangle.

[Bruno Luong]

@SA-W I think the basis functions are not wrong, but finding the reference coordinates of a point in real space and working just with the basis functions on (-1,1)x(-1,1) is an alternative to your bilinearformula(x, y, Z, xq, yq), right?

Right. When I make the comment I overlook the fsolve part to map the rectangle to (-1,1)x(-1x1).

I would never do that. You want to know the formula for bilinear interpolation to understand how it works, and then you use a numerical function that does not provide any formula. Fortunately @Stephen23 give you a big help to figure out the problem of your code.

But more interesting, I just think you approach could be applied where the grid is topological equivalent to generic quadrilateral mesh, such as pixels of some non-linear camera projection (fisheye) and not grided data.

BTW for rectangular gridded the cross-term x*y of the mapping vanish, and essentially the solution of fsolve can be determined separately for x and y, and if you plug the formula together you find my formula, consist, compact and fast to compute.

For quatrilateral I would look for analytical inversion rather than a big fsolve hammer, if it exists.

[Torsten]

syms x_left x_right y_low y_high
syms u_left_low u_right_low u_left_high u_right_high
syms x y
% First way
% Linearly interpolate between (x_left,y_low,u_left_low) and (x_right,y_low,u_right_low);
f_low(x) = u_left_low*(x-x_right)/(x_left-x_right) + u_right_low*(x-x_left)/(x_right-x_left);
% Linearly interpolate between (x_left,y_high,u_left_high) and (x_right,y_high,u_right_high);
f_high(x) = u_left_high*(x-x_right)/(x_left-x_right) + u_right_high*(x-x_left)/(x_right-x_left);
% Linearly interpolate between f_low(x) and f_high(x)
f1_inter(x,y) = f_low(x)*(y-y_high)/(y_low-y_high) + f_high(x)*(y-y_low)/(y_high-y_low)

f1_inter(x, y) = 

% Second way
% Linearly interpolate between (x_left,y_low,u_left_low) and (x_left,y_high,u_left_high);
f_left(y) = u_left_low*(y-y_high)/(y_low-y_high) + u_left_high*(y-y_low)/(y_high-y_low);
% Linearly interpolate between (x_right,y_low,u_right_low) and (x_right,y_high,u_right_high);
f_right(y) = u_right_low*(y-y_high)/(y_low-y_high) + u_right_high*(y-y_low)/(y_high-y_low);
% Linearly interpolate between f_left(y) and f_right(y)
f2_inter(x,y) = f_left(y)*(x-x_right)/(x_left-x_right) + f_right(y)*(x-x_left)/(x_right-x_left)

f2_inter(x, y) = 

x_left = 3;
x_right = 3.3;
y_low = 3;
y_high = 3.25;
u_left_low = 0;
u_right_low = 0.1121;
u_left_high = 0.064604;
u_right_high = 0.184942;
xq = 3.15;
yq = 3.125;
double(subs(f1_inter(xq,yq)))

ans = 0.0904

double(subs(f2_inter(xq,yq)))

ans = 0.0904

xq = 3.2;
yq = 3.2;
double(subs(f1_inter(xq,yq)))

ans = 0.1308

double(subs(f2_inter(xq,yq)))

ans = 0.1308