# finite difference method not working when going to two dimensions

I have two coupled ordinary differential equations in the steady state:

The following code solves, using the Jacobi finite difference method, in 1d using Dirichlet boundary conditions for function $$s$$ and Neumann boundary conditions for function $$n$$.

from pylab import imshow,gray,show, plot, xlim
import time

# Constants
M = 200         # Grid squares on a side

target = 1e-7   # Target accuracy

# Create arrays to hold  values
s = zeros([M+1],float)   # zeros everywhere
sprime = empty([M+1],float)   #  next  configuration
n = zeros([M+1],float)   # zeros everywhere
nprime = empty([M+1],float)   # s next  configuration

# Initialize grid
Lx = 20
x = linspace(-Lx, Lx, M+1)

a = x[2] - x[1]
print("a = ", a)

sigma =2.
source = exp(-(x**2) / (2 * sigma**2))

taus = 0.1
D = 0.5
b = -0.5

# Main loop
delta = 1.0
tic = time.time()

while delta>target:

# s has Dirichlet bcs
sprime[0] = 0.0
sprime[M] = 0.0

# n has Nuemann bcs
nprime[0] = n[1]
nprime[M] = n[M-1]

den = 2*D + a*a/taus
denn = 2*D
sprime[1:M] = (D*(s[0:M-1] + s[2:M+1]) + b * a*(n[2:M+1]-n[0:M-1])/2  + 2*source[1:M]*a*a)/den
nprime[1:M] = (D*(n[0:M-1] + n[2:M+1]) + b * a*(s[2:M+1]-s[0:M-1])/2 )/denn

# Calculate maximum difference from old values
deltas = max(abs(s-sprime))
deltan = max(abs(n-nprime))
delta = max([deltas, deltan])

# Swap the two arrays around
s,sprime = sprime,s
n,nprime = nprime,n

toc = time.time() - tic
print(toc)

# Make a plot
plot(x,s,'b')   #
plot(x,n,'r')  #
xlim(-Lx,Lx)

show()


The solutions looks great and I can confirm with analytic solution.

I now want to take the same equation but in two dimensions (allow for 2d diffusion and source function is now 2d):

This is where I'm stuck. My code attempt is

from pylab import imshow,gray,show, plot, xlim
import time

# Constants
M = 500         # Grid squares on a side
N = 500         # Grid squares on a side

target = 1e-14   # Target accuracy

# Create arrays to hold potential values
s = zeros([M+1,N+1],float)   # zeros everywhere
sprime = empty([M+1,N+1],float)   # phiprime is next phi configuration
n = zeros([M+1,N+1],float)   # zeros everywhere
nprime = empty([M+1,N+1],float)   # phiprime is next phi configuration

# Initialize grid
Lx = 50
Ly = 50
x = linspace(-Lx, Lx, M+1)
y = linspace(-Ly, Ly, N+1)
X, Y = meshgrid(x, y)

a = x[2] - x[1]

sigma = 0.2
source = exp(-(X**2 + Y**2) / (2 * sigma**2))

taus = 10000
D = 1.
b = -0.5

# Main loop
delta = 1.0
tic = time.time()

while delta>target:

# Dirichlet bc for s
sprime[0,:] = 0.0
sprime[M,:] = 0.0
sprime[:,0] = 0.0
sprime[:,N] = 0.0

#Neumann b.c. for n
nprime[0,:] = nprime[1,:]
nprime[M,:] = nprime[M-1,:]
nprime[:,0] = nprime[:,1]
nprime[:,N] = nprime[:,N-1]

den = 2*D + a*a/taus
denn = 2*D
sprime[1:M, 1:N] = (D*(s[0:M-1, 1:N] + s[2:M+1, 1:N] + s[1:M, 0:N-1] + s[1:M, 2:N+1]) + b * a*(n[2:M+1, 1:N]-n[0:M-1, 1:N])/2 + 2*source[1:M, 1:M]*a*a)/den
nprime[1:M, 1:N] = (D*(n[0:M-1, 1:N] + n[2:M+1, 1:N] + n[1:M, 0:N-1] + n[1:M, 2:N+1]) + b * a*(s[2:M+1, 1:N]-s[0:M-1, 1:N])/2)/denn

# Calculate maximum difference from old values
deltas = max(abs(s-sprime))
deltan = max(abs(n-nprime))
delta = max([deltas, deltan])

# Swap the two arrays around
s,sprime = sprime,s
n,nprime = nprime,n

toc = time.time() - tic
print(toc)

# Make a plot
plot(n[50,:],'b')


I've played a lot with the tolerance and the grid spacing but getting a lot of overflows. I'm hoping to just reproduce what I had in 1d with the addition of diffusion in 2d. The extension seems straightforward but I'm getting garbage.

• I think your diffusion terms should be divided by 4 in the 2D case. In the 1D case you divided by 2 but there is no division in your code for the 2D case Feb 28 at 17:01
• Thanks for catching that. I think it's still a factor of 2 to divide by and not 4. The other factors of 2 is in either den or denn Feb 28 at 19:27