2D lattice with fixed boundary condition

Question

I am writing code for the following equation with fixed boundary condition on a 2 dimensional lattice of \$L\times L\$ sites:

$$\begin{align} x_{i+1} =&\ (1-\varepsilon)r,円 x_i (1-x_i) + \\ &\ 0.25\varepsilon\left((r,円x_{i-1}(1-x_{i-1}) + r,円x_{i+1}(1-x_{i+1}) + r,円x_{i-L}(1-x_{i-L}) + r,円x_{i+L}(1-x_{i+L})\right) \end{align}$$

Fixed boundary condition means for end sites there are no neighboring sites beyond boundary.

Is there a simpler or more sophisticated way to write following code for the above equation with fixed boundary condition ?

def CML2d(x):
 eps = 0.3
 r = 3.9
 xn = np.zeros(N+1, float)
 for i in range(1, N+1):
 if i>L and i<=(L-1)*L:
 if i%L==1:
 xl, xr = 0., x[i+1]
 xu, xd = x[i-L], x[i+L]
 elif i%L==0:
 xl, xr = x[i-1], 0.
 xu, xd = x[i-L], x[i+L]
 else:
 xl, xr = x[i-1], x[i+1]
 xu, xd = x[i-L], x[i+L]
 elif i>1 and i<L:
 xl, xr = x[i-1], x[i+1]
 xu, xd = 0., x[i+L]
 elif i>(L-1)*L+1 and i<L*L:
 xl, xr = x[i-1], x[i+1]
 xu, xd = x[i-L], 0.
 elif i==1:
 xl, xr = 0., x[i+1]
 xu, xd = 0., x[i+L]
 elif i==L:
 xl, xr = x[i-1], 0.
 xu, xd = 0., x[i+L]
 elif i==(L-1)*L+1:
 xl, xr = 0., x[i+1]
 xu, xd = x[i-L], 0.
 elif i==L*L:
 xl, xr = x[i-1], 0.
 xu, xd = x[i-L], 0.
 xn[i] = (1-eps)*r*x[i]*(1-x[i]) + 0.25*eps*( r*xl*(1-xl) + r*xr*(1-xr) + r*xu*(1-xu) + r*xd*(1-xd) )
 return xn
L = 10 #side of 2d lattice
N = L*L #number of sites in 2d lattice
x0 = numpy.random.uniform(0.1, 0.9, N+1) #initial values for x
xf = [] # store iterate x
x = x0
for nt in np.arange(0.005, 50.005, 0.005):
 x = CML2d(x)
 xf.append(x)

• hpaulj: Other than creating the initial arrays, I don't see any use of numpy. It looks like the code would work just as well with lists. And the arrays are 1d, despite this being a 2d problem (shape (N,) rather than (L,L)).
• hpaulj: This is NOT working code! It does not work as posted.
• hpaulj: stackoverflow.com/questions/40946058/… is a similar problem - but with a 2d array with a wrapping boundary condition.
• hpaulj: Sorry for filling up the comments - but for a 10x10 problem, x and xn are 101 long - you ignore the first elements. It may make flat indexing a bit easier but it doesn't help with thinking in array terms.
• ajaydeep: @hpaulj I have not pasted the full code here but just the defined function. Code would work with list also, I know that. But I will be going for larger lattice size (say 1000*1000) in that case numpy array is better choice and also making code with 2 dimensional array would be computationally expensive for larger lattice size. Once I find the more optimistic solution I will change the indexing starting from zero.

Answer 1

Testing script from a few days ago:

import numpy as np

original:

def CML2d(x):
 eps = 0.3
 r = 3.9
 xn = np.zeros(N+1, float)
 for i in range(1, N+1):
 if i>L and i<=(L-1)*L:
 if i%L==1:
 xl, xr = 0., x[i+1]
 xu, xd = x[i-L], x[i+L]
 elif i%L==0:
 xl, xr = x[i-1], 0.
 xu, xd = x[i-L], x[i+L]
 else:
 xl, xr = x[i-1], x[i+1]
 xu, xd = x[i-L], x[i+L]
 elif i>1 and i<L:
 xl, xr = x[i-1], x[i+1]
 xu, xd = 0., x[i+L]
 elif i>(L-1)*L+1 and i<L*L:
 xl, xr = x[i-1], x[i+1]
 xu, xd = x[i-L], 0.
 elif i==1:
 xl, xr = 0., x[i+1]
 xu, xd = 0., x[i+L]
 elif i==L:
 xl, xr = x[i-1], 0.
 xu, xd = 0., x[i+L]
 elif i==(L-1)*L+1:
 xl, xr = 0., x[i+1]
 xu, xd = x[i-L], 0.
 elif i==L*L:
 xl, xr = x[i-1], 0.
 xu, xd = x[i-L], 0.
 value = (1-eps)*r*x[i]*(1-x[i]) + 0.25*eps*( r*xl*(1-xl) + r*xr*(1-xr) + r*xu*(1-xu) + r*xd*(1-xd) )
 xn[i] = value
 return xn

partial attempt to work with 2d x; the intention was to replace all the uses of L with 2d array indexing tests.

def CML2d_1(x):
 eps = 0.3
 r = 3.9
 L,_ = x.shape
 N = L*L
 # xn = np.zeros((L,L), float)
 xn = (1-eps)*r*x*(1-x)
 x = x.flat
 #for i in range(N):
 for j in range(L):
 for k in range(L):
 i = k+L*j
 i1 = i+1
 if i1>L and i1<=(L-1)*L:
 if i1%L==1:
 xl, xr = 0., x[i+1]
 xu, xd = x[i-L], x[i+L]
 ....
 #value = (1-eps)*r*x[i]*(1-x[i])
 value = 0.25*eps*( r*xl*(1-xl) + 
 r*xr*(1-xr) + 
 r*xu*(1-xu) + 
 r*xd*(1-xd) )
 #xn.flat[i] = value
 xn[j,k] += value
 return xn

But I then realized that I don't need to iterate over all points. Instead I could just sum subarrays (similar to the array of taking a 1d diff, x[1:]-x[:-1]:

# x[i+1] = (1-eps)* r*x[i]*(1-x[i]) + 
# 0.25*eps*( r*x[i-1]*(1-x[i-1]) + 
# r*x[i+1]*(1-x[i+1]) + 
# r*x[i-L]*(1-x[i-L]) + 
# r*x[i+L]*(1-x[i+L]) )
def CML2d_2(x):
 eps = 0.3
 r = 3.9
 x2 = r * x * (1-x)
 xn = (1-eps) * x2
 xn[1:,:] += 0.25 * eps * x2[:-1,:]
 xn[:-1,:] += 0.25 * eps * x2[1:,:]
 xn[:,1:] += 0.25 * eps * x2[:,:-1]
 xn[:,:-1] += 0.25 * eps * x2[:,1:]
 return xn

then making use of a covolve2d function in scipy.signal:

from scipy import signal
def CML2d_3(x):
 eps = 0.3
 r = 3.9
 in2 = np.zeros((3,3))
 in2[1,:] = 0.25 * eps
 in2[:,1] = 0.25 * eps
 in2[1,1] = (1-eps)
 print(in2)
 x2 = r * x * (1-x)
 res = signal.convolve2d(x2, in2, mode='same', boundary='fill', fillvalue=0)
 return res 

testing:

L = 10 #side of 2d lattice
#L = 4
N = L*L #number of sites in 2d lattice
x0 = np.random.uniform(0.1, 0.9, N+1) #initial values for x
x0[0]=np.nan # test if this is used
res=CML2d(x0.tolist())
print(res)
x2d = x0[1:].reshape(L,L)
res1=CML2d_1(x2d)
print(x0.shape)
print(res1.shape)
print(np.allclose(res1.flatten(), res[1:]))
print(np.allclose(res1, CML2d_2(x2d)))
print(np.allclose(res1, CML2d_3(x2d)))

• ajaydeep: Thank you for giving an attempt to problem. However, I liked CML2d_3() approach. Meanwhile I found an another solution using numpy.lib.pad which I will post soon.