# 2D lattice with fixed boundary condition

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)


## migrated from stackoverflow.comDec 3 '16 at 14:08

This question came from our site for professional and enthusiast programmers.

• 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 Dec 3 '16 at 18:53
• This is NOT working code! It does not work as posted. – hpaulj Dec 3 '16 at 18:56
• stackoverflow.com/questions/40946058/… is a similar problem - but with a 2d array with a wrapping boundary condition. – hpaulj Dec 3 '16 at 19:08
• 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. – hpaulj Dec 3 '16 at 21:04
• @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. – ajaydeep Dec 4 '16 at 4:20

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

• 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. – ajaydeep Dec 8 '16 at 10:15