Python function to solve Ax = b by inverse substitution

Question

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

Write a Python function to solve Ax = b by inverse substitution, where A is the upper triangular nonsingular matrix. The MATLAB code for this is on page 190, which you can use as a pseudo-code guide if you want. The function should take A and b as input and return x. Your function should not check that A is non-singular. That is, suppose that only non-singular A. will be passed to your function.

The MATLAB code to which it refers is:

x(n) = c(u)/U(n,n) for i = n-1 : -1 : 1 x(i) = c(i); for j = i+1 : nx(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, has an upper triangular test matrix (not sure if it is non-singular!). How can I check the 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 get is [[1 1 0]] for x. I'm not sure I'm doing it right. I guess this is wrong and cannot 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]; 

endless ... and never fulfilled.

your fubared indexing:

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

the range is half open, unlike matlab

Answer 2

One of the problems I see is that your input is made up of integers, which means Python will do integer division on them, which will turn 3/4 to 0 when you want floating point division. You can say that python does default floating point division by adding

 from __future__ import division 

To the top of your code. From using scipy, I assume you are using Python 2.x here.

Answer 3

You ask how to check the singularity of the upper triangular matrix?

Please do not calculate the qualifier!

Just look at the diagonal elements. Which of them are equal to zero? Any zero?

What about an effective numerical feature? Compare the smallest absolute value with the largest absolute value. If this ratio is less than something of the order of eps, it is effectively singular.