CODEWITHSUNDEEP
pythonlistnumpyfor-looplinear-algebra
Trevie3
Trevie3

Reputation: 15

Values in Python list won't change in for loop

I'm doing LU decomposition for a certain type of tridiagonal matrix and I'm having a problem with Python that I can't understand at all.

def trisolve2(b):
    length = len(b)

    # Solve Ly = b
    y = b.copy()
    for i in range(1, length):
        print('y[' + str(i) + '] = ' + str(b[i]) + ' - (-' + str(i) + ' / (' + str(i+1)
            + ')) * ' + str(y[i-1]))
        y[i] = b[i] - (-1 * i / (i+1)) * y[i-1]
    print("y = " + str(y))

    # Solve Ux = y
    x = y.copy()
    for i in range(length - 2, 0, -1):
        print("x = " + str(x))
        x[i] =  (y[i] + x[i+1]) / ((i+2) / (i+1))

    return x

test_b = np.array([[1], [0], [0], [0], [1]])
print("b = " + str(test_b))
test_x = trisolve2(test_b)
print("x = " + str(test_x))

I would expect when i = 1, for y[i] to be assigned to 1/2, and it even looks like it should with the print statements I put in there, but for some reason I can't wrap my head around it never changes. y remains b, and x remains y without changing. The output is:

b = [[1]
 [0]
 [0]
 [0]
 [1]]
y[1] = [0] - (-1 / (2)) * [1]
y[2] = [0] - (-2 / (3)) * [0]
y[3] = [0] - (-3 / (4)) * [0]
y[4] = [1] - (-4 / (5)) * [0]
y = [[1]
 [0]
 [0]
 [0]
 [1]]
x = [[1]
 [0]
 [0]
 [0]
 [1]]
x = [[1]
 [0]
 [0]
 [0]
 [1]]
x = [[1]
 [0]
 [0]
 [0]
 [1]]
x = [[1]
 [0]
 [0]
 [0]
 [1]]```

Upvotes: 1

Views: 79

Answers (1)

Adehad
Adehad

Reputation: 579

Your numpy array is defaulting to int32, 0.5 is getting rounded to the nearest integer: 1.

>> b = np.array([[1], [0], [0], [0], [1]])
>> b
array([[1],
       [0],
       [0],
       [0],
       [1]])
>> b.dtype
dtype('int32')

Try passing the dtype argument:

test_b = np.array([[1], [0], [0], [0], [1]], dtype=float)

Upvotes: 2

Related Questions