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);
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);
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