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)

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