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

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


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)

  • $\begingroup$ Hi Simon. Thanks for the reply. I do use three points for the second derivative... For the perforated domain, my thinking was that having a neighbour (e.g. u[x+1][y]) as a hole was like having a Neumann boundary condition, such that u[x+1][y]-u[x][y] = 0... Also, any thoughts on why the last expression for the cross-derivative is wrong, even when all the neighbours are part of the domain? Many thanks! $\endgroup$ – Marta Sep 26 '18 at 12:28
  • $\begingroup$ Oh yeah. Didn't catch that error before. See the edits. $\endgroup$ – Simon Sep 26 '18 at 12:57
  • $\begingroup$ Thanks, @Simon. I am glad our initial formulae for the cross-derivative also match. My second formula was indeed asymmetric and included u_{xx} and u_{yy}, which I had never seen before and whose accuracy I was wondering about... $\endgroup$ – Marta Sep 26 '18 at 13:07
  • $\begingroup$ In the diffusion setting, my "holes" are like regions that the diffusing substance cannot enter, i.e. there can be no flux between them and the rest of the domain. My difficulty is figuring out how to implement this into the derivatives by always writing them as functions of u(neighbour)-u that I can effectively zero when the neighbour is a hole. Hence my clumsy attempts at rewriting the cross-derivative. $\endgroup$ – Marta Sep 26 '18 at 13:12
  • $\begingroup$ Your reply nevertheless gave me the idea of using use-sided schemes for the cross-derivatives. I will try that. Thank you again, @Simon! $\endgroup$ – Marta Sep 26 '18 at 13:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.