3
$\begingroup$

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

enter image description here

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

enter image description here

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

enter image description here

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.

$\endgroup$
2
  • 1
    $\begingroup$ 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 $\endgroup$
    – whpowell96
    Feb 28 at 17:01
  • 1
    $\begingroup$ 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 $\endgroup$
    – BeauGeste
    Feb 28 at 19:27

0

Your Answer

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