反替换求解Ax=b的Python函数

Question

好的，对于我的数值方法课，我有以下问题：

编写一个Python函数，通过反代换求解Ax=b，其中a是上三角非奇异矩阵。这方面的MATLAB代码在第190页，如果您愿意，可以将其用作伪代码指南。函数应该接受输入A和b并返回x。您的函数不需要检查A是否非奇异。也就是说，假设只有非奇异的A才会传递给函数。

它引用的MATLAB代码是：

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

我的Python代码是用MATLAB代码片段编写的，上面有一个三角形测试矩阵（不确定它是否是非奇异的！如何测试奇点？）以下内容：

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)

我得到的答案是x的[[110]]，我不确定我做的是否正确。我想这是错的，不知道下一步该怎么办。有线索吗？

Answer 1

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

是无限的。。。永远不会被处决

你的索引功能很强大：

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

注意，范围是半开的，不像matlab

Answer 2

你问如何检验上三角矩阵的奇异性？

请不要计算行列式！

只需看看对角线元素。哪些是零？有零吗？

有效的数值奇异性如何？比较最小绝对值和最大绝对值。如果这个比率小于eps阶上的某个值，那么它实际上是奇异的。

Answer 3

我看到的一个问题是，您的输入由整数组成，这意味着Python将对它们进行整数除法，这将把3/4变成0，而您需要的是浮点除法。默认情况下，可以通过添加

from __future__ import division

在你的代码顶部。从使用scipy来看，我假设您在这里使用的是Python 2.x。