Re: interp2 meshgrid problem

[Juan Pablo Carbajal]

On Sat, Apr 28, 2012 at 6:51 PM, Grant Stephens <address@hidden> wrote:
> Hi
>
> Running version 3.2.4 on Ubuntu 12.04 and have come across something that
> doesn't make sense to me.
> The following program gives the error that the data is not in meshgrid
> format but all the xx, yy and xxx, yyy data is generated from the meshgrid
> command. I suspect it is becuase it is a non-linear meshgrid? It there a
> workaround or it a limitation?
>
> Thank You
> Grant
>
> % Set up grids and tensor product Laplacian and solve for u:
>   N = 24; [D,x] = cheb(N); y = x;
>   [xx,yy] = meshgrid(x(2:N),y(2:N));
>   xx = xx(:); yy = yy(:);       % stretch 2D grids to 1D vectors
>   f = 10*sin(8*xx.*(yy-1));
>   D2 = D^2; D2 = D2(2:N,2:N); I = eye(N-1);
>   L = kron(I,D2) + kron(D2,I);                       % Laplacian
>   figure(1), clf, spy(L), drawnow
>   tic, u = L\f; toc          % solve problem and watch the clock
>
> % Reshape long 1D results onto 2D grid:
>   uu = zeros(N+1,N+1); uu(2:N,2:N) = reshape(u,N-1,N-1);
>   [xx,yy] = meshgrid(x,y);
>   value = uu(N/4+1,N/4+1);
>
> % Interpolate to finer grid and plot:
>   [xxx,yyy] = meshgrid(-1:.04:1,-1:.04:1);
>   uuu = interp2(xx,yy,uu,xxx,yyy,'cubic');
>   figure(2), mesh(xxx,yyy,uuu), colormap(1e-6*[1 1 1]);
>   xlabel x, ylabel y, zlabel u
>   text(.4,-.3,-.3,sprintf('u(2^{-1/2},2^{-1/2}) = %14.11f',value))
>
>  function [D,x] = cheb(N)
>   if N==0, D=0; x=1; return, end
>   x = cos(pi*(0:N)/N)';
>   c = [2; ones(N-1,1); 2].*(-1).^(0:N)';
>   X = repmat(x,1,N+1);
>   dX = X-X';
>   D  = (c*(1./c)')./(dX+(eye(N+1)));      % off-diagonal entries
>   D  = D - diag(sum(D'));
>
> _______________________________________________
> Help-octave mailing list
> address@hidden
> https://mailman.cae.wisc.edu/listinfo/help-octave
>

Hi,
I tried the code you sent and in Octave 3.6.1 I do not get the same
error but the interpolation fails. However interp2 works perfectly
x=linspace(-1,1,10); [
xx yy] = meshgrid(x,x);
x2=linspace(-1,1,100);
[xx2 yy2] = meshgrid(x,x);

zz = sin(xx*2*pi) + cos(yy*2*pi);
uu = interp2 (xx,yy,zz,xx2,yy2);
mesh(xx2,yy2,uu)

uu = interp2 (xx,yy,zz,xx2,yy2,'cubic');
mesh(xx2,yy2,uu)

uu = interp2 (xx,yy,zz,xx2,yy2,'pchip');
mesh(xx2,yy2,uu)

uu = interp2 (xx,yy,zz,xx2,yy2,'spline');
mesh(xx2,yy2,uu)

I think the problem is with your "uu" data and the cubic method (and
pchip) since this works fine (uu is your vector)

[xx2 yy2]=meshgrid(linspace(-1,1,100));
uu2 = interp2 (xx,yy,uu,xx2,yy2,'spline');
mesh(xx2,yy2,uu2)

It might be a bug in interp2 but I would first double check that uu is OK.

-- 
M. Sc. Juan Pablo Carbajal
-----
PhD Student
University of Zürich
http://ailab.ifi.uzh.ch/carbajal/