Bilinear and Trilinear Interpolation

I am writing matlab code for Bilinear and Trilinear Interpolation (Numerical Aalysis/Methods). It uses following formula to approximate a function at reqired values of x, y and z

I have a code which calculates all the Lagrange Coefficients but I'm not able to substitute all the values in above formula. Also, all the values of vectors x, y and z are yto be substituted in f(x,y,z) as can be seen in the formula

syms x y Z
x = [0.05 0.06];
y = [1.23 1.41];
z = [0.1 0.2];
w=length(x);
u=length(y);
t=length(z);
l=w-1;
L=zeros(w,w);
for k=1:l+1
   A=1;
   for j=1:l+1
      if k~=j
         A=conv(A,poly(x(j)))/(x(k)-x(j));
      end
   end
   L(k,:)=A;
end
L
M=zeros(u,u);
m=u-1;
for k=1:m+1
   B=1;
   for j=1:m+1
      if k~=j
         B=conv(B,poly(y(j)))/(y(k)-y(j));
      end
   end
   M(k,:)=B;
end
M
n=t-1;
N=zeros(t,t);
for k=1:n+1
   C=1;
   for j=1:n+1
      if k~=j
         C=conv(C,poly(z(j)))/(z(k)-z(j));
      end
   end
   N(k,:)=C;
end
N

Matrices L, M and N represent Lagrange polynomials in each row.

This is what my code looks like.

Can anyone help me writing the code for substituting all the values required in above formula.

[John D'Errico]

Is there a reason why you are not using interp2, interp3 or just interpn? It is rarely a good idea to write inefficient code to do that which is already provided.

[Jalal Hassan]

I did not know about these. I,m all ears now. Please guide me

[Jalal Hassan]

Also, please help me to use the formula which I put in the question to write the code.

Answer 1

[John D'Errico]

Wow. I guess I've been working with those tools for so many years, as well as building variations of them, that I think of them as second nature.

What is often called bilinear interpolation, when done on a lattice is what interp2 or griddedInterpolant calls linear. And, I guess it may not be obvious that what one tool calls linear, is what someone else might call bilinear (for example, photoshop, as I recall has a bilinear interpolation tool for image interpolation) or someone else might call tensor product linear interpolation. But the formula you have written is essentially that tensor product linear interpolant. So in any dimension, the result will be a piecewise linear function.

I'll show how to use griddedInterpolant, which I prefer myself. This allows you to learn to use one tool, in various numbers of dimensions.

In the help for griddedInterpolant, we see:

F = griddedInterpolant(...,METHOD) specifies the interpolation method:
      'linear'   - (default) linear, bilinear, trilinear,... interpolation

So when we ask griddedInterpolant for 'linear' interpolation in multiple dimenssions, it does a bilinear or trilinear interpolation, as you wish.

x = 0:0.5:1;
y = linspace(0,1,4);

[xx,yy] = ndgrid(x,y)
xx =
            0            0            0            0
          0.5          0.5          0.5          0.5
            1            1            1            1
yy =
            0      0.33333      0.66667            1
            0      0.33333      0.66667            1
            0      0.33333      0.66667            1

Now create an array z, of the form z(x,y). So two independent variables, and thus a bilinear interpolant.

zz = xx.^2 + yy.^2
zz =
            0      0.11111      0.44444            1
         0.25      0.36111      0.69444         1.25
            1       1.1111       1.4444            2

z(x,y) is pretty boring, I'll admit. Sorry about that. It is a paraboloid of revolution, in two dimensions.

GI = griddedInterpolant(xx,yy,zz,'linear')
GI = 
  griddedInterpolant with properties:

            GridVectors: {[0 0.5 1]  [0 0.33333 0.66667 1]}
                 Values: [3×4 double]
                 Method: 'linear'
    ExtrapolationMethod: 'linear'

Can I convince you that GI is a bilienar interpolant? For example, fix x at some value.

yint = linspace(0,1,100);
xint = reshape(0.2,size(yint));
plot(yint,GI(xint,yint),'-r')
xlabel 'yint'
ylabel 'z(xint,yint)'
grid on

So 3 linear segments. This is a tensor product linear (bilinear) interpolant, as you asked about.

With 3 independent variables, so trilinear interpolation, you would create xx,yy,zz using ndgrid. Then create a 3-d array of values that defines the relation w(x,y,z). Now griddedInterpolant will do the trilinear interpolation you want.

If you prefer to use interp2 interp3, or interpn, you can also do so, but these tools require arrays that are ordered as meshgrid would create, instead of ndgrid.

Be careful though. One of the things people need to understand about these interpolants is they are not truly linear when viewed along any path through space. For example,

[xx,yy] = ndgrid([0 1],[0 1])
xx =
     0     0
     1     1
yy =
     0     1
     0     1
     
zz = [0 .5 ; .5 2]
zz =
            0          0.5
          0.5            2

Now, consider the path that moves acrosss the diagonal of this cell, running from the points (0,0) to (1,1), with x==y along that path.

GI = griddedInterpolant(xx,yy,zz,'linear');
xint = linspace(0,1,100);
yint = xint;
zint = GI(xint,yint);
plot(xint,GI(xint,yint),'-r')
xlabel 'xint == yint'
ylabel 'z(xint,yint)'
grid on

So even though this is a bilinear interpolant, when you vary two variables at once, the shape of the function can exhibit a quadratic polynomial shape.

In the case of a triinear interpolant, you will see cubic polynomials when traversed along a diagonal. This is completely as expected. Except of course, it seems to come as a surprise to many, probably because this is a "linear" interpolation. I will admit the first time I saw this behavior, it seemed surprising to me. Once I learned how to manipulate these tools, my next trick was to build a 3-d array, such that when interpolated along the diagonal using trilinear interpolation had the complete behavior of a user supplied cubic spline. Totally useless of course, unless you wanted exactly that behavior, which sometimes, we did. :)

[Jalal Hassan]

Thank you very much :)