Cross derivatives using finite differences I am numerically calculating the derivatives of a scalar function u(x,y) in a domain defined in a 2D-Cartesian grid (x,y) implementing finite differences. I have been using a centred scheme so far, but would now like to account for the possibility that my domain may have holes. For example, u(x+1,y) may not be assigned a value and thus cannot be used to compute u(x,y) or u(x+2,y).
For the first derivatives, I replace the centered scheme for a one-sided scheme as needed (and set it to 0 if both relevant neighbours are missing). For the second derivative, I can rewrite: 
dudx2 = (u[x - 1][y] + u[x + 1][y] - 2 * u[x][y]) / (dx*dx);
as
dudx2 = (isdomain[x-1][y]*(u[x-1][y]-u[x][y]) + isdomain[x+1][y]*(u[x+1][y]-u[x][y])) / (dx*dx);
I am struggling to write the cross-derivative in an appropriate format. Before, I had:
dudxdy = (u[x + 1][y + 1] + u[x - 1][y - 1] - u[x + 1][y - 1] - u[x - 1][y + 1])/(4*dx*dx);
By expanding u[x+1][y+1] and u[x-1][y-1] as a Taylor series, it seems to be I should be able to write:
dudxdy = (isdomain[x-1][y-1]*(u[x-1][y-1]-u[x][y]) + isdomain[x+1][y+1]*(u[x+1][y+1]-u[x][y])) / (2*dx2) - 0.5*(dudx2+dudy2);
But some quick numerical tests tell me this is not right. (Why not?)
Do you have any suggestions on how I can get round this problem? The context is the numerical solution of an anisotropic reaction-diffusion equation in this perforated domain.
Many thanks!
Marta
 A: For the first derivative, you are right. You can either use
$$\frac{\partial u(x)}{\partial x} = \frac{1}{2dx}(u[x+1] - u[x-1]) + O(dx^2)$$
or
$$\frac{\partial u(x)}{\partial x} = \frac{1}{dx}(u[x+1] - u[x]) + O(dx)$$
(The first version is generally preferred due to the smaller error term). Note that in both cases you need to evaluate $u$ at exactly two points (though more could be used to improve accuracy, but thats not the point here).
Now for the second derivative, you will always need do evalue (at least) three points. The best numerical accuracy will be achieved with the symmetric formula
$$ \frac{\partial^2u(x)}{\partial x^2} = \frac{1}{dx^2}(u[x+1]-2u[x]+u[x-1]) + O(dx^2)$$
But in principle any three points can be used (if they are too far away from $x$, the accuracy will be very bad). For example at a border you might use
$$ \frac{\partial^2u(x)}{\partial x^2} = \frac{1}{dx^2}(u[x+2]-2u[x+1]+u[x]) + O(dx)$$
This means that you formula for the second derivative is wrong (even in the 1-dimensional case). It is impossible to estimate the second derivative using only two points.
EDIT 1 (cross-derivative):
Your cross-derivative formula is incorrect even without any holes. The commenly used symmetric one is
$$\frac{\partial^2 u(x,y)}{\partial x\partial y} = \frac{\partial}{\partial y}\frac{\partial u(x,y)}{\partial x}\\=\frac{\partial}{\partial y}\left(\frac{1}{2dx}(u[x+1][y] - u[x-1][y])\right) = \frac{1}{4dx dy}\left(u[x+1][y+1] - u[x+1][y-1] - u[x-1][y+1] + u[x-1][y-1]\right) + \text{error-terms}$$
(You might note that these are four terms. Theoretically, a version with three terms is possible, but on a rectangular grid this would always be asymmetric and thus have bad accuracy)
EDIT 2 (your intuition about Neumann-Boundary):
There is another basic approach to dealing with holes: removing them by giving them a value and then use the nice symmetric formulas for derivatives. There are multiple approaches of doing so, for example:


*

*just set $u[x] = 0$ (Dirichlet boundary)

*set it equal to a neighbour $u[x] = u[x+1]$ (Neumann boundary, which you effectively did)

*interpolation (if you have values on both sides of the hole): $u[x] = \frac{1}{2} (u[x+1]+u[x-1])$
If you know which one of these is sensible for your problem, go ahead and use it. But this depends very much on the physical problem you are trying to model. I.e. what does a hole in your system mean? Should the physical sytem actually have a hole there, or is just the simulation missing a point? (I dont know your problem well enough to answer this)
