# Correct Interpolation on a pre-calculated grid that is a vector function of a vector?

6 views (last 30 days)
Brendan on 9 Jan 2020
Commented: Brendan on 10 Jan 2020
I have calculated vector data, such that the data can be written and represented as , where each component is a function of all three coordinates; i.e. . I need to interpolate to be able to represent the inputs of a different vector anywhere, and have it be a vector output. This different vector is . Written out fully, you can see that each component of the vector M is an array which depends on all three x, y, and z coordinates. As a result, I can't simply use Interp3 three times. This means that I need the G's as a grid, so that I can query at points (Gx_q, Gy_q, Gz_q) and find their (x,y,z) inputs that give G.
Interp3 repeatedly gives errors of "Input is not a Meshgrid."
My vector data arrays (, , and ) are sized 31x31x31. I'm so lost; I've tried so many tricks that have failed. Any tricks? Will scatteredInterpolant() work?
The situation is analogous to this user's question (https://www.mathworks.com/matlabcentral/answers/408712-interpolating-scattered-3-dimensional-data), but distinctly different. The answers provided in () also are not quite related – something there is wrong, as well.
Previously, I've poorly or incorrectly explained my problem, and so those questions (and, as a result, answers) were phrased incorrectly. This is fixed.
To avoid confusion, I'll present the case where the M vector depends only on one coordinate, each, such that , which lets us write a scalar G as , and just do the process three times. The code would look like:
x = 0:1000:30000;
y = 0:1000:30000;
z = 0:1000:30000;
Gx_q = 721.0085; %These are example query points
Gy_q = 15.1313;
Gz_q = 45.9302;
Gx = x - squeeze(Mx(1,1,:))'; %These should be 1D arrays; for simplicity, in this case I would only need one column of each M_i array.
Gy = y - squeeze(My(1,1,:))';
Gz = z - squeeze(Mz(1,1,:))';
Gx_out = interp1(Gx,x,Gx_q,'linear','extrap');
Gy_out = interp1(Gy,y,Gy_q,'linear','extrap');
Gz_out = interp1(Gz,z,Gz_q,'linear','extrap');

Matt J on 9 Jan 2020
Edited: Matt J on 9 Jan 2020
so that I can query at points (Gx_q, Gy_q, Gz_q) and find their (x,y,z) inputs that give G.
But there may exist no such point or many of them. Suppose for example that Gx=1,Gy=1,Gz=1 everywhere in x,y,z. Then if Gx_q=5, Gx_q=6,Gz_q=7 there will exist no (x,y,z) agreeing with this query point for any reasonable interpolator. If, on the other hand Gx_q=1, Gx_q=1,Gz_q=1 then all (x,y,z) agree with it.
The best I think you can do is solve for one coordinate (x,y,z) at which the interpolation of G is closest to (Gx_q, Gy_q, Gz_q) according to some error metric. That can be approached as an inverse problem as in the following example:
[x,y,z]=ndgrid(0:1000:30000);
Gx=x-Mx;
Gy=y-My;
Gz=z-Mz;
%Find best grid point (x,y,z) - later refine that with fminsearch
Err = abs(Gx(:)-Gx_q)+abs(Gy(:)-Gy_q)+abs(Gz(:)-Gz_q);
[~,i]=min(Err);
xyzInitial = [x(i), y(i), z(i)];
%Refine:
Fx=griddedInterpolant(x,y,z,Gx);
Fy=griddedInterpolant(x,y,z,Gy);
Fz=griddedInterpolant(x,y,z,Gz);
xyz=fminsearch(@(xyz)lookupError(xyz,Fx,Fy,Fz), xyzInitial);
x_out=xyz(1);
y_out=xyz(2);
z_out=xyz(3);
function Err=lookupError(xyz,Fx,Fy,Fz)
gx=Fx(xyz);
gy=Fy(xyz);
gz=Fz(xyz);
Err= abs(gx-Gx_q)+abs(gy-Gy_q)+abs(gz-Gz_q);
end

Show 1 older comment
Matt J on 9 Jan 2020
Try it now.
Brendan on 9 Jan 2020
It seems to work. Holy moly. This was several days' struggle. Thank you so much!
Brendan on 10 Jan 2020
Last question: To minimize the error function, would you think it smarter to package this entire function into its function to be passed to fminsearch, or would it be better to simply make an ndgrid with a smaller stepsize?
Thanks!