$\begingroup$

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$?

$\endgroup$ 1

1 Answer

$\begingroup$

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.

$\endgroup$ 1