Python function to solve Ax = b by back substitution.

Question

Okay, for my numerical methods class I have the following question:

Write a Python function to solve Ax = b by back substitution, where A is an upper triangular nonsingular matrix. MATLAB code for this is on page 190 which you can use as a pseudocode guide if you wish. The function should take as input A and b and return x. Your function need not check that A is nonsingular. That is, assume that only nonsingular A will be passed to your function.

The MATLAB code that it refers to is:

x(n) = c(u)/U(n,n) 
for i = n-1 : -1 : 1
  x(i) = c(i);
  for j = i+1 : n 
    x(i) = x(i) - U(i,j)*x(j);
    end
  x(i) = x(i)/U(i,i);
end

My Python code, which I wrote using the MATLAB code snippet, is with an upper triangular test matrix(not sure if its nonsingular! How do I test for singularity?):

from scipy import mat
c=[3,2,1]
U=([[6,5,1],[0,1,7],[0,0,2]])
a=0
x=[]
while a<3:
    x.append(1)
    a=a+1

n=3
i=n-1
x[n-1]=c[n-1]/U[n-1][n-1]
while i>1: 
    x[i]=c[i]
    j=i+1
    while j<n-1:
        x[i]=x[i]-U[i][j]*x[j];
    x[i]=x[i]/U[i][i]
    i=i-1
print mat(x)

The answer I am getting is [[1 1 0]] for x. I not sure if I am doing this correctly. I assume it is wrong and can't figure out what to do next. Any clues?

Answer 1

j=i+1
while j<n-1:
    x[i]=x[i]-U[i][j]*x[j];

is infinite ... and never gets executed

your indexing is fubared:

for i in range(n-2,-1,-1):
....
    for j in range(i+1,n):

notice, range is half open unlike matlab

Answer 2

You ask how to test for singularity of an upper triangular matrix?

Please don't compute the determinant!

Simply look at the diagonal elements. Which ones are zero? Are any zero?

How about effective numerical singularity? Compare the smallest absolute value to the largest in absolute value. If that ratio is smaller than something on the order of eps, it is effectively singular.

Answer 3

One problem I see is that your input consists of integers, which means that Python is going to do integer division on them, which will turn 3/4 into 0 when what you want is floating point division. You can tell python to do floating point division by default by adding

from __future__ import division

To the top of your code. From the use of scipy, I'm assuming you're using Python 2.x here.