How can I implement cubic interpolation in Gradient Descent algorithm?

After computing the gradient of z along the entire computational grid, I want to calculate Z0=Z(X0,Y0) and gradient(Z0) using cubic interpolation. I'm unsure of the best way to do this but it is essential that X0 and Y0 are set outside of the function. They are coordinates of an arbitrary starting point to begin the Gradient Descent algorithm.

Below is the problem part of my code with some attempted statements:

%Function to implement the 2-Dimensional Gradient Descent method. The
%output of this function should be an [N x 3] matrix containing, on each
%row, the x, y and z co-ordinates of every position considered in the
%Gradient Descent algorithm (including start and end points)

function [xi, yi, Z] = gradient_descent2(Z,X0,Y0,gamma,tau)%Declare function
x = -9:0.2:9;% Set the input of comp_z function two be two vectors with a set stepsize
y = -8:0.2:8;
[Xg, Yg] = meshgrid(x,y);  %make mesh grid
[px,py] = gradient(Z); %compute gradient of z along entire computational
%mesh grid (z computed outside).

%---------------------------------2--------------------------------------%
%Use cubic interpolation to obtain starting point z(X0,Y0) and gradient at
%this point
Z0 = griddata(Xg,Yg,Z,X0,Y0,'cubic');
t = (x == X0) & (y == Y0);
indt = find(t);
Z0_grad = [px0(indt) py0(indt)];

This is supposed to be a function that implements a 2D Gradient Descent method. x, y & z are related as follows:

and is computed outside the function.

Answer 1

[Matt J]

Since you have an analytical formula for z(x,y), why don't you just use that to calculate Z and its gradient exactly? Why are you working with discrete approximations and interpolaton?

In other words, why don't you just do as follows?

       Z = @(x,y)  cos(x/2)*cos(y/2) - x/5 +y/10;
  gradZ = @(x,y) [-sin(x/2)*cos(y/2)/2 -1/5    ;   -cos(x/2)*sin(y/2)/2  + 1/10 ];

[OzzWal]

That would be a lot easier for me - I have written a function comp_z which computes z.

However, the above is the way that the rubrik seems to want me to do it (this is homework!) and is integral to the way the algorithm works. I'm clear on steps 3&4 just step 2 I'm really unsure of what query points z is being interpolated from and how to achieve it:

[Matt J]

OK. Well, I think it should look like this,

xgrid = -9:0.2:9;% Set the input of comp_z function two be two vectors with a set stepsize
ygrid = -8:0.2:8;
[Gx,Gy] = gradient(Z);

Z0=interp2(xgrid,ygrid,Z,X0,Y0,'cubic');
gradZ0_X=interp2(xgrid,ygrid,Gx,X0,Y0,'cubic');
gradZ0_Y=interp2(xgrid,ygrid,Gy,X0,Y0,'cubic');

[OzzWal]

Thank you!