Back substitution in Python

Question

import numpy as np


A = np.matrix([[1,1,-1],[0,1,3],[0,0,-6]])
b = np.array([9,3,8])


def back_sub(A,b):
    n = len(A)
    print('n is:', n)
    x = [0]*n
    for i in range(n-1,-1,-1): #this refers to the rows; i goes 2,1,0
        for j in range(i+1,n): #j goes 1,2 @ i = 0
                               #j goes 2   @ i = 1
            b[i] = b[i] - A[i,j]*x[j]
        x[i] = b[i]/A[i,i]

    return x

I am trying to back substitute an upper-triangular matrix A, and come up with an answer: $x_1 = 4/3, x_2 = 7, x_3 = -4/3$.

When I run my code, it gives correct answer for $x_2$ and $x_3$, but $x_1$ comes out as 0, which is confounding. I have tried to debug my code, and it seems to run fine until $i = 0$. Why isn't this code giving me the correct value of $x_1$?

Answer 1

I added a few lines of code:

        print(b[i])
        b[i] = b[i] - A[i,j]*x[j]
        print(i,A[i,j],x[j])
        print(b[i])
        print('--')

the output of the final iteration is:

--
2
0 -1 -1.33333333333
0
--

Apparently, it thinks $2 - (-1\cdot -4/3) = 0$. This is due to intermediate rounding. With normal python 3 code this does not happen, but apparently with numpy you still have to be careful. Adding "b = b.astype(float)" on top resolves the issue.