3
$\begingroup$

I would like to model the 2D diffusion equation with Neumann BC's inside the following egg-shaped domain:

IMG

I would like to use the finite difference method with the discretization implied by the image above, but while I can handle the interior points, I don't know how to handle the boundary values. I know how to search for and identify the boundary points, but I'm not sure what finite-difference approximation to the no-flux condition I should use. Currently I'm thinking of estimating the direction of the normal at each boundary point by looking at it's 8 nearest neighbors, and then, if the normal is approximately

$$\textbf{n} = [n_1, n_2]$$

I would add an equation like:

$$\mathbf{\nabla} c \cdot \textbf{n} \approx \frac{c_{i,j}-c_{i-1,j}}{\Delta x} \times n_1 + \frac{c_{i,j}-c_{i,j-1}}{\Delta y} \times n_2 = 0$$

Though of course I'll need to be careful which side I take the derivative from.

Does anyone have any suggestions or resources for this type of problem? Examples in 3D would be particuarly useful, as that's where I'm working towards.

$\endgroup$
9
  • $\begingroup$ What is the domain? The red egg? $\endgroup$ Oct 1, 2014 at 10:58
  • $\begingroup$ Yes, sorry, I'll change that. $\endgroup$
    – tom
    Oct 2, 2014 at 0:27
  • $\begingroup$ There are plenty of methods to deal with boundary conditions on oddly shaped domains using finite difference methods. What literature have you already looked for? $\endgroup$ Oct 2, 2014 at 2:59
  • $\begingroup$ I can't say I've found a lot :) My experience with FDM's doesn't extend much beyond Numerical Recipes. The problem I actually want to solve involves modelling anisotropic diffusion inside a 3D brain, and the boundary there is already discretized for me and is terribly irregular. Most of what I've found so far relates to solving PDE's on a domain with a known function (e.g. on a ellipse) and you make clever discretizations. But I can't do that: my discretization is given. I would very much appreciate recommendations of textbooks/papers dealing with this issue :) $\endgroup$
    – tom
    Oct 3, 2014 at 0:44
  • $\begingroup$ Try google.com/… It shows 179,000 hits. $\endgroup$ Oct 3, 2014 at 1:51

1 Answer 1

2
$\begingroup$

It seems like you'll get the boundary conditions you want automatically if you construct your discrete Laplacian in terms of an adjacency matrix, minus the sum of the number of neighbors along the diagonal.

Consider the 1D case:

  • In the interior, the discrete Laplacian is $\{1,-2,1\}$. That is, the adjacent cells $\{1,0,1\}$ minus their sum ($-2$) along the diagonal.
  • At i.g. a right boundary, the cell to the right is missing so you get $\{1,-1,0\}$ for the local operation

The analogous construction in 2D is similar. Following the guidance here, you can include the four corners in your Laplacian in a way that maximizes rotational symmetry within the constraints of a square grid, which will help with the "corners"

Here's some scipy code exploring this. Note how pixels outside the boundary are masked out of the adjacency matrix before calculating the value of the diagonal for the discrete Laplacian operator.

from pylab import *
import scipy

# Define square grid and region mask
L       = 64
grid    = linspace(-1,1,L)[:,None] + 1j*linspace(-1.5,1.5,L)[None,:]
mask    = abs(grid)**2<.9
nanmask = float32([NaN,1])[int32(mask)]
imshow(mask,cmap='gray');
title('Region to simulate'); axis('off'); 
savefig('1.png')
show()

# Build a discrete Laplacian operator
A1d = eye(L,k=-1) + eye(L,k=1) + eye(L,k=L-1) + eye(L,k=1-L) 
A2d_cross  = scipy.sparse.kronsum(A1d,A1d).toarray()
A2d_corner = scipy.sparse.kron(A1d,A1d).toarray()
Ad = A2d_cross*2/3+A2d_corner*1/3
Ad[~mask.ravel(),:] = Ad[:,~mask.ravel()] = 0
Lp = Ad - diag(sum(Ad,0))

# Simplify adjacency matrix for speed
M = float32(Lp)
keep = any(M>0,0) | any(M>0,1)
M = M[keep][:,keep]
NPOINTS = sum(keep)

# Get eigenmodes
NMODES = 49
η,V  = scipy.sparse.linalg.eigsh(scipy.sparse.bsr_matrix(M),
                                 k=NMODES,which='LA')
# Plot eigenmodes
figure(figsize=(5,5))
subplots_adjust(0,0,1,1,0,0)
k = int(ceil(sqrt(NMODES)))
for i in range(NMODES):
    subplot(k,k,i+1)
    the_mode = zeros((L,L))
    the_mode[keep.reshape(L,L)] = V[:,i]
    scale = np.max(abs(the_mode))
    imshow(the_mode*nanmask,vmin=-scale,vmax=scale,cmap='coolwarm')
    axis('off')
suptitle('Eigenmodes—Neumann boundary zero')
savefig('2.png')
show()

And its figure outputs:

enter image description here

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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