《计算方法》---列主元高斯消元法（python实现）

计算方法的第三单元的第一个算法，列主元高斯消元法
我将它分为3个部分
第一部分，输入所要求的的方程式组
第二部分，消元，简而言之进行排序
第三部分，回代求解，解决问题

这里通过这个例子来协助编程

1. 第一部分：

import numpy as np                  #本次仅需要用到这一个库
#a = input('需要计算的方程组为：')
a = np.array([[10,-2,-1,3],[-2,10,-1,15],[-1,-2,5,10]])      #双括号   a即要计算的

i = a.shape[0]                       #得到a的行数
j = a.shape[1]                       #得到a的列数
print(i)        
print(j)
X = []                               #建立一个X空列表，便于放入所得解
print(X)

2.第二部分—消元部分（出的问题最多）：
分两小部分
第一小部分排序 ，这里提供了一种思路，先找到最大值，然后分别与最大值比较

for z in range(0,j-2):                   #1.当z从第1列到j-1列循环，每一次循环中，都是针对所在列，先排序，后消元（以下对应这两部分） 2.z还代表循环次数     
    for x in range(z,i):                 #x从z行到第最后一行循环，每次循环扔掉上次循环的第一行，因为每次的第一行都用来存储每次的最大
        n = np.argmax(abs(a[x:i,z]))      #n为剩下的几行中，绝对值最大的值所在行的索引
        print(n+x)                        #***索引再加x，因为n为剩下的几行中的索引，而我们接下来的交换针对所有行
        a[[x,n+x],:] = a[[n+x,x],:]      #数组特殊的交换数值方式，吃亏
        print(a)
    
    for y in range(1+z,i):              #之所以从1+z，开始，最开始的z循环不仅仅是列，也是循环次数，每循环一次，向下移动一列                
        k = a[z,z]/a[y,z]               #充分理解这个式子
        print(k)
        a[y,] = a[z,]-k*a[y,]         
print(a)

output:
0
[[10 -2 -1 3]
[-2 10 -1 15]
[-1 -2 5 10]]
1
[[10 -2 -1 3]
[-2 10 -1 15]
[-1 -2 5 10]]
2
[[10 -2 -1 3]
[-2 10 -1 15]
[-1 -2 5 10]]
-5.0
-10.0
1
[[ 10 -2 -1 3]
[ 0 48 -6 78]
[ 0 -22 49 103]]
2
[[ 10 -2 -1 3]
[ 0 48 -6 78]
[ 0 -22 49 103]]
-2.18181818182
[[ 10 -2 -1 3]
[ 0 48 -6 78]
[ 0 0 100 302]]

3.第三部分—回带求解
回带求解的思路：
q从最后一行开始，循环到第一行，每一行的时候分别计算X
w仅在分母a[q,w]中用到，这里需要随着行数的上升，同时所在列向左移动，所以列在每次循环时候都-1
因为求的解X[n]一开始是空的，在不停增加，所以不能每次n循环，循环的次数一样，否则报错。所以通过h来限制改变每一次的循环

#for w in range(j-2,-1,-1):           #这里从倒数第二列开始向前循环   列！  注意：步长默认为1，所以改成-1
w = j-1
g = 0
for q in range(i-1,-1,-1):  #i为矩阵行数  这里是从最后一行开始向上循环   行！
    w = w - 1    
    h = 2 - q              #注意！这里的2，在不同的方程组是不一样，这是3元，所以2
    
    if h == 0:
        X.append(a[q,j-1]/a[q,w])
        print(X)
    else:
        for n in range(0,h):               #e+1决定让n循环几次
            g = g + a[q,2 - n] * X[n]
            print(g)
        X.append((a[q,j-1] - g)/a[q,w])   #公式
        print(X)
        g = 0                            #很重要，容易遗漏！推一推就好

output:
[3.02]
-18.12
[3.02, 2.0024999999999999]
-3.02
-7.025
[3.02, 2.0024999999999999, 1.0024999999999999]

反思总结提升：
1.a = np.array([[10,-2,-1,3],[-2,10,-1,15],[-1,-2,5,10]]) #双括号,格式还是要熟记得
2.i = a.shape[0] #得到a的行数
j = a.shape[1] #得到a的列数
3.np.argmax(abs(a[x:i,z])) numpy库里的argmax函数，找到最大值的索引
！！！需要注意的是，它返回的索引值是！argmax()括号里的范围，从0开始计数
！4.a[[x,n+x],:] = a[[n+x,x],:] #数组特殊的交换数值方式，吃亏（与单一数值交换不同）。如果用学C或者PYTHON特有的数值交换方式，总是会存在最后一行不变的异样
例：n = np.argmax(abs(a[x:i,y]))
A = a[x,]
a[x,] = a[n,]
a[n,] = A #在这里对换了数值
5.IndexError: list index out of range #如果取了X[n]没有的位置，会报错的