Matrix indexing and initialization in " for in" loop

[QuintenR]

Hi all, in order to make a MATLAB script more generally available, I'm rewriting it in Python.
Google couldn't provide me with conclusive answers, therefore I've set up a testfile that reproduces my problem (see below) and hope that someone here can help me out.

The results (matrix "x") seem to depent on whether the matrix "xnew" is (re-)initialized in line 20.
To my best knowledge, this should not be possible since all values of "xnew" will be overwritten in the for loops starting at lines 22 and 23.

Am I missing something? I hope someone can tell me more about this initialization problem.
Thank you in advance!

import numpy as np
nnx = int(33)
nny = int(33)
h = 0.03
volfrac = 0.4
dcdx = np.ones((nnx,nny))
x = np.tile(volfrac, (nnx, nny))
xnew = np.zeros((nnx,nny))
Am = h**2
loop = 1
while loop < 3:
    l1 = 0.
    l2 = 10000. 
    move = 0.2
    eta = 0.5

    while (l2-l1) >  1e-4:
            lmid = 0.5*(l2+l1)
            V = 0.00
            xnew = np.zeros((nnx,nny)) #Results are different when this line is commented. Why?
            
            for j in range(0,nny):
                for i in range(0,nnx):
              
                    Bc = dcdx[i,j]/(lmid*Am)
                    
                    if x[i,j]*(Bc**eta) <= max([(1.0-move)*x[i,j],0.001]):
                        xnew[i,j] = max([(1.0-move)*x[i,j],0.001])
                    elif x[i,j]*(Bc**eta) >= min([(1.0 + move)*x[i,j],1.0]):
                        xnew[i,j] = min([(1.0 + move)*x[i,j],1.0])
                    else:
                        xnew[i,j] = x[i,j]*(Bc**eta)

                    V = V + Am*xnew[i,j]
                    
            if (V - volfrac*(nnx-1.0)*(nny-1.0)*h**2.) > 0.0: 
                l1 = lmid
            else:
                l2 = lmid
                    
    x = xnew
    loop = loop + 1
print(x)

The script is part of a bi-sectioning algorithm that converges V to the value "volfrac*(nnx-1.0)*(nny-1.0)*h**2".

[deanhystad]

At the bottom of your while loop you set x = xnew and then you use x in the next iteration.

When you call "xnew = np.zeros((nnx,nny))" you create a new matrix filled with zeros. xnew references this new matrix and x references the old matrix.

If you don't reassign xnew to a new matrix, xnew and x both reference the same matrix. When you change values in xnew you are also changing values in x.

This is better demonstrated by the following short example.

x = [1, 2, 3, 4]
y = x
x.append(5)
print(y)

Output:
[1, 2, 3, 4, 5]
>>> id(x)
2714481334144
>>> id(y)
2714481334144

Even though 5 was appended to x, it gets printed when I print y. This is because x and y both reference the same list. For further proof I used id() to print the Python object referenced by x and y. Notice they are the same.

Now I will change my program slightly so x and y do not reference the same object. This is what you are doing when you reassign xnew to an array of zeroes.

x = [1, 2, 3, 4]
y = x.copy()
x.append(5)
print(y)

Output:
[1, 2, 3, 4]
>>> id(x)
1753300211136
>>> id(y)
1753288343232

This time printing y does not display the change made to x, and printing the id's shows that x and y reference different objects.

[QuintenR]

Thank you Dean, that was a very helpful reply!
Coming from a MATLAB environment and being relatively new to Python, I was not aware of this referencing system. It makes a lot of sense though, in terms of memory.

I will dive into this subject to get a better understanding of referencing and indexing in Python. Thank you again!