在python中旋转的高斯消除
声明：本页面是StackOverFlow热门问题的中英对照翻译，遵循CC BY-SA 4.0协议，如果您需要使用它，必须同样遵循CC BY-SA许可，注明原文地址和作者信息，同时你必须将它归于原作者（不是我）：StackOverFlow
原文地址: http://stackoverflow.com/questions/31957096/
Warning: these are provided under cc-by-sa 4.0 license. You are free to use/share it, But you must attribute it to the original authors (not me):
StackOverFlow
Gaussian elimination with pivoting in python
提问by Olli
I am trying to write a function that will solve a linear system using gaussian elimination with pivoting. I am not allowed to use any modules either. Can someone help me out here? I don't know what I'm doing wrong. Thanks in advance :) :) :)
我正在尝试编写一个函数，该函数将使用带旋转的高斯消元来求解线性系统。我也不允许使用任何模块。有人可以帮我吗？我不知道我做错了什么。提前致谢 ：） ：） ：）
First i set up augmented matrix M, then i do the pivoting and row operations and finally i do the back substitution
首先我设置了增广矩阵 M，然后我进行了旋转和行操作，最后我进行了反向替换
def linearsolver(A,b):
n = len(A)
M = A
i = 0
for x in M:
x.append(b[i])
i += 1
for k in range(n):
for i in range(k,n):
if abs(M[i][k]) > abs(M[k][k]):
M[k], M[i] = M[i],M[k]
else:
pass
for j in range(k+1,n):
q = M[j][k] / M[k][k]
for m in range(k, n+1):
M[j][m] += q * M[k][m]
x = [0 for i in range(n)]
x[n] =float(M[n][n+1])/M[n][n]
for i in range (n-1,-1,-1):
z = 0
for j in range(i+1,n):
z = z + float(M[i][j])*x[j]
x[i] = float(M[i][n+1] - z)/M[i][i]
print x
回答by izidor
Let's make it educational as it is clearly your homework and explain the mistakes in the code, shall we? Maybe you can learn something in the process.
让我们让它具有教育意义，因为它显然是您的作业并解释代码中的错误，好吗？也许你可以在这个过程中学到一些东西。
Assume there is the following matrix to solve:
假设有以下矩阵需要求解：
A = [[3, 2, -4], [2, 3, 3], [5, -3, 1]]
b = [3, 15, 14]
Your initial code:
您的初始代码：
def linearsolver(A,b):
n = len(A)
M = A
i = 0
for x in M:
x.append(b[i])
i += 1
for k in range(n):
for i in range(k,n):
if abs(M[i][k]) > abs(M[k][k]):
M[k], M[i] = M[i],M[k]
else:
pass
for j in range(k+1,n):
q = M[j][k] / M[k][k]
for m in range(k, n+1):
M[j][m] += q * M[k][m]
x = [0 for i in range(n)]
x[n] =float(M[n][n+1])/M[n][n]
for i in range (n-1,-1,-1):
z = 0
for j in range(i+1,n):
z = z + float(M[i][j])*x[j]
x[i] = float(M[i][n+1] - z)/M[i][i]
print x
As first you get an error, something about index out of range or something
首先你会得到一个错误，关于索引超出范围之类的东西
Traceback (most recent call last):
File "solver.py", line 32, in <module>
print linearsolver([[3, 2, -4], [2, 3, 3], [5, -3, 1]], [3, 15, 14])
File "solver.py", line 24, in linearsolver
x[n] =float(M[n][n+1])/M[n][n]
IndexError: list index out of range
How you can debug this? Print the current state of the calculation:
你怎么能调试这个？打印计算的当前状态：
x = [0 for i in range(n)]
print "n = ", n
print "x = ", x
for row in M:
print row
x[n] =float(M[n][n+1])/M[n][n]
You will get output:
你会得到输出：
n = 3
x = [0, 0, 0]
[5, -3, 1, 14]
[2, 3, 3, 15]
[3, 2, -4, 3]
Traceback (most recent call last):
File "solver.py", line 37, in <module>
print linearsolver([[3, 2, -4], [2, 3, 3], [5, -3, 1]], [3, 15, 14])
File "solver.py", line 29, in linearsolver
x[n] =float(M[n][n+1])/M[n][n]
IndexError: list index out of range
If you count it correctly, you try to write to x[3]and access M[3][4]and M[3][3]. However, python counts from 0, meaning that the last element is -1 smaller than expected. Last element are x[2]and M[2][3]. This can be fixed by applying -1 to all indexes:
如果计算正确，则尝试写入x[3]和访问M[3][4]和M[3][3]。但是，python 从 0 开始计数，这意味着最后一个元素比预期小 -1。最后一个元素是x[2]和M[2][3]。这可以通过对所有索引应用 -1 来解决：
x[n-1] =float(M[n-1][n])/M[n-1][n-1]
Run it and....
运行它并......
n = 3
x = [0, 0, 0]
[5, -3, 1, 14]
[2, 3, 3, 15]
[3, 2, -4, 3]
Traceback (most recent call last):
File "solver.py", line 37, in <module>
print linearsolver([[3, 2, -4], [2, 3, 3], [5, -3, 1]], [3, 15, 14])
File "solver.py", line 34, in linearsolver
x[i] = float(M[i][n+1] - z)/M[i][i]
IndexError: list index out of range
Damn, another traceback! However, we're making progress, its on higher line now! When you look at it closer, you use n+1again. Originally Ahad nrows and ncolumns. Because you appended bto the matrix, there are n+1columns. Because of indexing from 0, the last element is -1 of the total columns: n+1-1 = nWe fix it and there is code:
妈的，又一个回溯！但是，我们正在取得进展，现在处于更高的水平！当你仔细观察它时，你会n+1再次使用它。最初A有n行和n列。因为您附加b到矩阵，所以有n+1列。由于从 0 开始索引，最后一个元素是总列数的 -1：n+1-1 = n我们修复它并且有代码：
def linearsolver(A,b):
n = len(A)
M = A
i = 0
for x in M:
x.append(b[i])
i += 1
for k in range(n):
for i in range(k,n):
if abs(M[i][k]) > abs(M[k][k]):
M[k], M[i] = M[i],M[k]
else:
pass
for j in range(k+1,n):
q = M[j][k] / M[k][k]·
for m in range(k, n+1):
M[j][m] += q * M[k][m]
x = [0 for i in range(n)]
print "n = ", n
print "x = ", x
for row in M:
print row
x[n-1] =float(M[n-1][n])/M[n-1][n-1]
for i in range (n-1,-1,-1):·
z = 0·
for j in range(i+1,n):
z = z + float(M[i][j])*x[j]·
x[i] = float(M[i][n] - z)/M[i][i]
print x
And when we run it, we get result!!!1!
当我们运行它时，我们得到了结果！！！1！
n = 3
x = [0, 0, 0]
[5, -3, 1, 14]
[2, 3, 3, 15]
[3, 2, -4, 3]
[6.4, 5.75, -0.75]
Just proof checking that it works:
只是证明检查它有效：
3 * 6.4 + 2 * 5.75 - 4 * -0.75 = 33.7 (not 3)
2 * 6.4 + 3 * 5.75 + 3 * -0.75 = 27.8 (not 15)
5 * 6.4 - 3 * 5.75 + 1 * -0.75 = 14.0 (correct)
Hmm... Our matrix wasn't solved but there at least we have one row solved properly. The algorithm requires for the final step to have matrix in certain format, where most rows starts with 0. But that's not the case as you can see. Let's add additional prints to show the matrix as we compute it:
嗯......我们的矩阵没有解决，但至少我们已经正确解决了一行。该算法要求最后一步具有特定格式的矩阵，其中大多数行以 0 开头。但正如您所见，情况并非如此。让我们添加额外的打印来显示我们计算的矩阵：
def linearsolver(A,b):
n = len(A)
M = A
i = 0
for x in M:
x.append(b[i])
i += 1
for k in range(n):
print "Iteration ", k
for i in range(k,n):
if abs(M[i][k]) > abs(M[k][k]):
M[k], M[i] = M[i],M[k]
else:
pass
# Show the matrix after swapping rows
for row in M:
print row
print
for j in range(k+1,n):
q = M[j][k] / M[k][k]
for m in range(k, n+1):
M[j][m] += q * M[k][m]
# Show matrix after multiplying rows
for row in M:
print row
print
x = [0 for i in range(n)]
print "n = ", n
print "x = ", x
for row in M:
print row
x[n-1] =float(M[n-1][n])/M[n-1][n-1]
for i in range (n-1,-1,-1):
z = 0
for j in range(i+1,n):
z = z + float(M[i][j])*x[j]
x[i] = float(M[i][n] - z)/M[i][i]
print x
And run it:
并运行它：
Iteration 0
[5, -3, 1, 14]
[2, 3, 3, 15]
[3, 2, -4, 3]
[5, -3, 1, 14]
[2, 3, 3, 15]
[3, 2, -4, 3]
Iteration 1
[5, -3, 1, 14]
[2, 3, 3, 15]
[3, 2, -4, 3]
[5, -3, 1, 14]
[2, 3, 3, 15]
[3, 2, -4, 3]
Iteration 2
[5, -3, 1, 14]
[2, 3, 3, 15]
[3, 2, -4, 3]
[5, -3, 1, 14]
[2, 3, 3, 15]
[3, 2, -4, 3]
n = 3
x = [0, 0, 0]
[5, -3, 1, 14]
[2, 3, 3, 15]
[3, 2, -4, 3]
[6.4, 5.75, -0.75]
As you can see, the row multiplication didn't happen at all! Let's inspect the state:
如您所见，行乘法根本没有发生！让我们检查一下状态：
for j in range(k+1,n):
q = M[j][k] / M[k][k]·
print "j=", j, "q=", q
print "M[j][k]=", M[j][k], "M[k][k]=", M[k][k]
for m in range(k, n+1):
M[j][m] += q * M[k][m]
For the iteration 0, we get:
对于迭代 0，我们得到：
Iteration 0
[5, -3, 1, 14]
[2, 3, 3, 15]
[3, 2, -4, 3]
j= 1 q= 0
M[j][k]= 2 M[k][k]= 5
j= 2 q= 0
M[j][k]= 3 M[k][k]= 5
[5, -3, 1, 14]
[2, 3, 3, 15]
[3, 2, -4, 3]
Wait a moment, qshouldn't be 0! It should be 0.4 and 0.6! This happens because both numbers are integers and python provided result as an integer instead of a float. To fix this, change it to float():
等一下，q应该不是0！应该是 0.4 和 0.6！发生这种情况是因为两个数字都是整数，而 python 提供的结果是整数而不是浮点数。要解决此问题，请将其更改为float()：
for j in range(k+1,n):
q = float(M[j][k]) / M[k][k]·
print "j=", j, "q=", q
print "M[j][k]=", M[j][k], "M[k][k]=", M[k][k]
for m in range(k, n+1):
M[j][m] += q * M[k][m]
Now we get result:
现在我们得到结果：
Iteration 0
[5, -3, 1, 14]
[2, 3, 3, 15]
[3, 2, -4, 3]
j= 1 q= 0.4
M[j][k]= 2 M[k][k]= 5
j= 2 q= 0.6
M[j][k]= 3 M[k][k]= 5
[5, -3, 1, 14]
[4.0, 1.7999999999999998, 3.4, 20.6]
[6.0, 0.20000000000000018, -3.4, 11.4]
Although it has a different output, we are missing zeros in the first column. Why? In place of 2, there is 4 instead of 0, In place of 3, there is 6 instead of 0. I know! You should subtract instead of adding multiplied row:
尽管它有不同的输出，但我们在第一列中缺少零。为什么？2 的位置是 4 而不是 0，3 的位置是 6 而不是 0。我知道！您应该减去而不是添加相乘的行：
Iteration 0
[5, -3, 1, 14]
[2, 3, 3, 15]
[3, 2, -4, 3]
j= 1 q= 0.4
M[j][k]= 2 M[k][k]= 5
j= 2 q= 0.6
M[j][k]= 3 M[k][k]= 5
[5, -3, 1, 14]
[0.0, 4.2, 2.6, 9.399999999999999]
[0.0, 3.8, -4.6, -5.4]
Iteration 1
[5, -3, 1, 14]
[0.0, 4.2, 2.6, 9.399999999999999]
[0.0, 3.8, -4.6, -5.4]
j= 2 q= 0.904761904762
M[j][k]= 3.8 M[k][k]= 4.2
[5, -3, 1, 14]
[0.0, 4.2, 2.6, 9.399999999999999]
[0.0, 0.0, -6.952380952380952, -13.904761904761903]
Iteration 2
[5, -3, 1, 14]
[0.0, 4.2, 2.6, 9.399999999999999]
[0.0, 0.0, -6.952380952380952, -13.904761904761903]
[5, -3, 1, 14]
[0.0, 4.2, 2.6, 9.399999999999999]
[0.0, 0.0, -6.952380952380952, -13.904761904761903]
n = 3
x = [0, 0, 0]
[5, -3, 1, 14]
[0.0, 4.2, 2.6, 9.399999999999999]
[0.0, 0.0, -6.952380952380952, -13.904761904761903]
[2.9999999999999996, 0.9999999999999996, 2.0]
Can you see those 0.0 at beginning of most rows? Let's proof check it now:
你能在大多数行的开头看到那些 0.0 吗？现在让我们证明检查一下：
3 * 2.9999999999999996 + 2 * 0.9999999999999996 - 4 * 2.0 = 2.9999999999999964 (almost 3)
2 * 2.9999999999999996 + 3 * 0.9999999999999996 + 3 * 2.0 = 14.999999999999998 (almost 15)
5 * 2.9999999999999996 - 3 * 0.9999999999999996 + 1 * 2.0 = 14.0 (correct)
Computers have issues with floating numbers but the algorithm got us to the correct solution [3, 1, 2]
计算机有浮点数问题，但算法让我们找到了正确的解决方案 [3, 1, 2]
After removing the helper prints, there is the code:
去除辅助打印后，有代码：
def linearsolver(A,b):
n = len(A)
M = A
i = 0
for x in M:
x.append(b[i])
i += 1
for k in range(n):
for i in range(k,n):
if abs(M[i][k]) > abs(M[k][k]):
M[k], M[i] = M[i],M[k]
else:
pass
for j in range(k+1,n):
q = float(M[j][k]) / M[k][k]
for m in range(k, n+1):
M[j][m] -= q * M[k][m]
x = [0 for i in range(n)]
x[n-1] =float(M[n-1][n])/M[n-1][n-1]
for i in range (n-1,-1,-1):
z = 0
for j in range(i+1,n):
z = z + float(M[i][j])*x[j]
x[i] = float(M[i][n] - z)/M[i][i]
print x
Things to take away:
带走的东西：
- Python counts indexes from 0, meaning that the last index is n-1
- Be aware of difference between integer and float division
- Print statements can help you debug the issues.
- Python从0开始计数索引，意思是最后一个索引是n-1
- 注意整数和浮点除法之间的区别
- 打印语句可以帮助您调试问题。
