Manually determining the gradient of a surface

Hi , I’m trying to write a function that can manually determine the gradient of a surface (represented by a multivariable function that I have created) at a particular point (x,y)

The code I've written is shown below.

function [pdx,pdy] = comp_gradient(x,y)

kernelx = [2.5,0,-2.5];

kernely = [2.5;0;-2.5];

z = cos(x/2).* cos(y) + (y/10) - (x/5);

pdx = conv2(kernelx,z,'same');

pdy = conv2(kernely,z,'same');

end

I'm not sure where to go from here, my output is clearly incorrect.

Any help would be greatly appreciated.

Answer 1

[John D'Errico]

Your attempt as you posted it is not actually that far off.

First, we need to know what the spacing is between elements, and in each direction. So what is an estimator of the first derivative of a function? Assume we have a function defined on a uniform spacing, with spacing dx and dy in the x and y directions.

First create a grid of points. Note that I used a different spacing in x and y to make things clear.

x = linspace(-1,1,20);
dx = x(2) - x(1);
y = linspace(-2,2,20);
dy = y(2) - y(1);
[X,Y] = meshgrid(x,y);
% and a function
fun = @(x,y) cos(x/2).* cos(y) + (y/10) - (x/5);
Z = fun(X,Y);

Now, you need to remember that meshgrid variex x as the SECOND index. I used meshgrid in this because we tend to think of X as the horizontal coordinate, Y as the vertical. But we will need to be careful, as the horizontal coordinate of an array is the second one. Anyway...

What is the first derivative in the x direction, or an estimate thereof? You used central differences, so lets stick with that. We would have something like this:

df/dx = (F(y, x + dx) - F(y , x - dx))/(2*dx)

df/dy = (F(y + dy, x) - F(y - dy , x))/(2*dy)

We can achieve that using a convolution, as you tried. The convolution kernel for the x and y gradients would be

kernelx = [-1 0 1]/(2*dx);
kernely = [-1;0;1]/(2*dy);

Now, how to use conv? Be careful, as edge effects will cause problems. conv2 will insert zeros for those terms that fall off the edge of the world. So I would strongly recommend the use of 'valid'.

pdx = conv2(Z,kernelx,'valid');
pdy = conv2(Z,kernely,'valid');

Note that this will cause the partial derivative estimate arrays to be 2 elements smaller in the indicated directions. We see that in the respective sizes of pdx and pdy.

size(pdx)
ans =
    20    18
size(pdy)
ans =
    18    20

If you wish, you could patch the edges, perhaps inserting a lower order approximation as new first and last columns for x, and new first and last rows for y.

pdx = [conv2(Z(:,1:2),[-1 1]/dx,'valid'), pdx, conv2(Z(:,[end-1,end]),[-1 1]/dx,'valid')];
size(pdx)
ans =
    20    20
surf(X,Y,pdx)

Then do something similar for the y partial.

So, as I said, you really were reasonably close in how you started.

[Carlos Pedro]

Thankyou John, this has helped so much. Happy new year my friend!

Answer 2

[Alexander Viktorov]

% Hi , I’m trying to write a function that can manually determine
% the gradient of a surface (represented by a multivariable 
% function that I have created) at a particular point (x,y) 
% The code I've written is shown below.
function [pdx, pdy] = comp_gradient(x, y, h)
% fun - char of function (ToDo)
% x - scalar numeric 
% y - scalar numeric coordinates of a particular point (x,y)
% h - step of each argument (the less h the greater precision)
% % % kernelx = [2.5,0,-2.5];
% hx = 2.5; % step of argument x
% kernelx = -2.5:hx:2.5;
		n = 5; % for example
		kernelx = (x - n*h):h:(x + n*h);
% % % kernely = [2.5;0;-2.5];
kernely = (y - n*h):h:(y + n*h);
[X, Y] = meshgrid(kernelx, kernely);% since (kernely == kernelx)

% % % z = cos(x/2).* cos(y) + (y/10) - (x/5);
z = cos(X/2).* cos(Y) + (Y/10) - (X/5);
% [pdx, pdy] = gradient(z, hx, hx);
[px, py] = gradient(z, h, h)
% returns a matrix (2n+1)by(2n+1) for each gradient's component
% take the result from center
pdx = px(n+1, n+1);
pdy = py(n+1, n+1);
% % % pdx = conv2(kernelx,z,'same');
% % % pdy = conv2(kernely,z,'same');
end
% I'm not sure where to go from here, 
% my output is clearly incorrect.
% Any help would be greatly appreciated.

%comp_gradientTest.m - Тестирование comp_gradient
% This script was generated by the AlVi Toolbox version 0.2.0.00
% development: 02-Jan-2019 23:48:05
clc; clear; close all;
format compact
name = 'comp_gradient' % имя тестируемой функции
disp('======  test of [pdx, pdy] = comp_gradient(x, y, h) ======')
x = 1
y = 3
h = 0.1
[pdx, pdy] = comp_gradient(x, y, h)
disp('== complete testd of comp_gradient ===')
help 'comp_gradient'
loc = which('comp_gradient')
% [fList, pList] = matlab.codetools.requiredFilesAndProducts('comp_gradient.m')
% shg %show graph window current
commandwindow

Answer 3

[Puro89]

Hi, how can I deal with the gradient of a image? I have this exercise that expressively ask to use the "conv2" function. Thanks in advance.