1
\$\begingroup\$

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)
\$\endgroup\$

migrated from stackoverflow.com Dec 3 '16 at 14:08

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

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

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.