Gradient Vector for partial derivative seems backward for dfdx

clc; clear;
xo = 0.8;
yo = 0.2;
%analytic solution
syms f(x,y)       %tells matlab to treat f(x) symbolically
f(x,y) = x*exp(-(x^2+y^2));
%Calculate Derivatives
Dfx = diff(f,x);      %partial of f w.r.t x
Dfy = diff(f,y);      %partial of f w.r.t y
fprintf('Analytic Solutions:\n')

Analytic Solutions:

fprintf('Dfx = %9.5f\n',Dfx(xo,yo))

Dfx =  -0.14185

fprintf('Dfy = %9.5f\n',Dfy(xo,yo))

Dfy =  -0.16212

%numerical solution
h = 0.01;        %grid size deltax = deltay = 0.01 
[nx,ny] = meshgrid(0:h:2);   %x and y range from 0 to 2, step size h
z = nx.*exp(-(nx.^2+ny.^2));
%locate index values for xo = 0.8 and yo = 0.2.
nxo = find(nx(1,:)==xo);
nyo = find(ny(:,1)==yo);

Calculate Derivatives

[nDfx, nDfy] = gradient(z,h);   %calculate the partials w.r.t x and y
fprintf('Numerical Solutions:\n')

Numerical Solutions:

fprintf('Dfx = %9.5f\n',nDfx(nxo,nyo))   %incorrect value

Dfx =   0.46604

fprintf('Dfx = %9.5f\n',nDfx(nyo,nxo))   %correct value, but why?

Dfx =  -0.14180

fprintf('Dfy = %9.5f\n',nDfy(nxo,nyo))

Dfy =  -0.16211

These are the outputs:

Analytic Solutions:

Dfx = -0.14185

Dfy = -0.16212

Numerical Solutions:

Dfx = 0.46604

Dfx = -0.14180

Dfy = -0.16211

Why must I switch the order for dfdx in order to get the correct result?

Answer 1

[Paul]

Hi Kevin,

It turns out that nDfy only gives the correct result because of the symmetry.

xo = 0.8;
yo = 0.2;
h = 0.01;        %grid size deltax = deltay = 0.01 
[nx,ny] = meshgrid(0:h:2);   %x and y range from 0 to 2, step size h
z = nx.*exp(-(nx.^2+ny.^2));
%locate index values for xo = 0.8 and yo = 0.2.
nxo = find(nx(1,:)==xo);
nyo = find(ny(:,1)==yo);
[nDfx, nDfy] = gradient(z,h);   %calculate the partials w.r.t x and y
nDfx(nxo,nyo)   %incorrect value

ans = 0.4660

nDfx(nyo,nxo)   %correct value, but why?

ans = -0.1418

nDfy(nxo,nyo) % happens to be the correct value

ans = -0.1621

nDfy(nyo,nxo)  % corret value

ans = -0.1621

When we use meshgrid, the x-values vary across the 2nd dimension (horizontal) and the y-values vary down the 1st dimension (vertical). Sounds counterintutive, but that's how it works. Similarly, gradient assumes the x-direction is horizontal across the array and the y-direction is vertical down the array, which is consistent with the meshgrid convention. But that means that indexing into the results also has to reverse as well, which is what we're seeing.

Consider a simpler problem

x = 1:3;
y = 101:103;
[X,Y] = meshgrid(x,y);
Z = 5*X + Y

Z = 3×3

106 111 116 107 112 117 108 113 118

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

If we want the value of Z that corresponds to the second element of x and the first element y (z = 5*2 + 101 = 111), we'd need Z(1,2)

[Kevin]

So what it appears I should have done is reverse how the values are labeled from meshgrid/gradient, namely:

[ny,nx] = meshgrid(0:h:2);

and

[nDfy, nDfx] = gradient(z,h);

Then, it works as expected (at least to me):

%locate index values for xo = 0.2 and yo = 0.3.

nyo = find(ny(1,:)==xo);

nxo = find(nx(:,1)==yo);

and

fprintf('Dfx = %9.5f\n',nDfx(nyo,nxo))

fprintf('Dfy = %9.5f\n',nDfy(nyo,nxo))

[Paul]

If it works and is less confusing for you, that's all that counts.

But this code looks a lot more confusing to me and I suspect would also to many others. The solution I showed was exactly what you had originally, except explained the need for the reversed indexing into nDfx and nDfy. Here, we have a bunch of other changes, but still have the reversed indexing into nDfx and nDfy.

How do you want the code to work? Perhaps have x vary vertically and y vary horizontally so that we'd index as nDfx(nxo,nyo)?

[Kevin]

I will play around with this further, but it seems odd to switch the order of the arrays. To me, it seems like you use f(xo,yo) to evaluate the function at (xo,yo), but dfdx(yo,xo) for the derivative. Why would you want to do that?

[Paul]

When using meshgrid, indexing into the function array is the same as indexing into the gradient arrays:

xo = 0.8;
yo = 0.2;
h = 0.01;        %grid size deltax = deltay = 0.01 
v = 0:h:2;
nxo = find(v==xo);
nyo = find(v==yo);
[nx,ny] = meshgrid(v);   %x and y range from 0 to 2, step size h
z = nx.*exp(-(nx.^2+ny.^2));
[z(nyo,nxo) xo.*exp(-(xo.^2 + yo.^2))]  % z(nyo,nxo)

ans = 1×2

0.4053 0.4053

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

If all we want is to index into the function array as z(nxo,nyo), then we can use ndgrid

[nx,ny] = ndgrid(v);
z = nx.*exp(-(nx.^2+ny.^2));
[z(nxo,nyo) xo.*exp(-(xo.^2 + yo.^2))]  % z(nxo,nyo)

ans = 1×2

0.4053 0.4053

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

However, gradient expects the input array to be in meshgrid format, so if z is in ndgrid format we'd have to go through some additional gymnastics.