I want to compute partial derivatives of complex functions via finite difference approximation on two dimensional grid using NDsolve`FiniteDifferenceDerivative function (I'm on V10 on OSX). Let construct a discrete, rectangular grid $(x_i, y_j)$, $i=1,\ldots,n_x$, $j=1,\ldots,n_y$
xg = Range[-2, 2, 0.1];
yg = Range[-2, 2, 0.025];
{nx, ny} = Length/@{xg, yg};
xyg = Flatten[Outer[List, xg, yg], 1];
for which we create FiniteDifferenceDerivativeFunction's
Do[FiniteDifferenceDerivative[i, j] =
NDSolve`FiniteDifferenceDerivative[Derivative[i, j], {xg, yg},
DifferenceOrder -> 4], {i, 0, 2}, {j, 0, 2}];
Let's take simple nontrivial example
f = Function[{x, y}, Cos[x*y] + I*Sin[x*y]];
Calculating $\partial_{x}f$ gives reasonable result. Problematic is $\partial_{y}f$
ListPlot3D[
Flatten /@
Transpose@{xyg,
Re@Flatten[
FiniteDifferenceDerivative[0, 1][Partition[f @@@ xyg, ny]]]},
PlotRange -> All]
For comparison we can plot analytic expression for $\partial_{y}f$
Plot3D[Re@Derivative[0, 1][f][x, y], {x, -2, 2}, {y, -2, 2}]
I've checked that this gives incorrect result for complex functions only. Another interesting thing is that FiniteDifferenceDerivative[0, 1][Partition[f @@@ xyg, ny]] seems to generate random results, compare
with previous plot. The quesiton is this a bug or I'm doing something wrong? Is this OS or Mathematica version (and so NDSolve version) dependent issue? One possible resolution is to take out differentiation matrices from FiniteDifferenceDerivativeFunction object and apply them by hand (I thought this is what is done when FiniteDifferenceDerivativeFunction is called, please correct me if I'm wrong), eg.
Transpose @ Dot[
FiniteDifferenceDerivative[0, 1][[3, 2, 1]],
Transpose @ Partition[f @@@ xyg, ny]]
to get $\partial_yf$ (the result can be plotted as before with ListPlot3D[Flatten /@ Transpose@{xyg, Re @ Flatten @ %}, PlotRange -> All]), but is this the only way to get correct result?
