考虑n阶线性方程组:
高斯消去法python
import numpy as np
def gaussin(a, b):
m, n = a.shape
c = np.zeros(n)
for i in range(n):
#限制条件
if (a[i][i] == 0):
print("no answer")
#k表示第一层循环,(0,n-1)行
#i表示第二层循环,(k+1,n)行,计算该行消元的系数
#j表示列
for k in range(n - 1):
for i in range(k + 1, n):
c[i] = a[i][k] / a[k][k]
for j in range(m):
a[i][j] = a[i][j] - c[i] * a[k][j]
b[i] = b[i] - c[i] * b[k]
x = np.zeros(n)
x[n - 1] = b[n - 1] / a[n - 1][n - 1]
#回代求出方程解
for i in range(n - 2, -1, -1):
for j in range(i + 1, n):
b[i] -= a[i][j] * x[j]
x[i] = b[i] / a[i][i]
for i in range(n):
print("x" + str(i + 1) + " = ", x[i])
if __name__ == '__main__':
a = np.array([[2, 4, -6], [1, 5, 3], [1, 3, 2]])
b = np.array([-4.0, 10.0, 5.0])
gaussin(a, b)
print('*' * 20)
但是,在消元过程中,无法使主元素a(ii)≠0,但是很小时,用其做除数,会导致其他元素数量级的严重增长,舍入误差的扩展,最后导致计算结果不可靠。所以这次采用列主元素消去法来进行,思想就是将有小数的那行与该列中数最大的那行进行交换。
python代码如下:
import numpy as np
#高斯列主元消去法
#找到主元并交换
def swap(a,b,k,n):
ans = -1.
for i in range(k,n):
if ans < np.fabs(a[i][k]):
ans = a[i][k]
maxn = i
a[[k,maxn],:] = a[[maxn,k],:]
b[k],b[maxn] = b[maxn],b[k]
def gaussin(a, b):
m, n = a.shape
c = np.zeros(n)
for i in range(n):
#限制条件
if (a[i][i] == 0):
print("no answer")
#k表示第一层循环,(0,n-1)行
#i表示第二层循环,(k+1,n)行,计算该行消元的系数
#j表示列
for k in range(n - 1):
swap(a,b,k,n)
for i in range(k + 1, n):
c[i] = a[i][k] / a[k][k]
for j in range(m):
a[i][j] = a[i][j] - c[i] * a[k][j]
b[i] = b[i] - c[i] * b[k]
x = np.zeros(n)
x[n - 1] = b[n - 1] / a[n - 1][n - 1]
#回代求出方程解
for i in range(n - 2, -1, -1):
for j in range(i + 1, n):
b[i] -= a[i][j] * x[j]
x[i] = b[i] / a[i][i]
for i in range(n):
print("x" + str(i + 1) + " = ", x[i])
if __name__ == '__main__':
a = np.array([[2, 4, -6], [1, 5, 3], [1, 3, 2]])
b = np.array([-4.0, 10.0, 5.0])
gaussin(a, b)
print('*' * 20)
a = np.array([[0.01, 2, -0.5], [-1, -0.5, 2], [5, -4, 0.5]])
b = np.array([-5.0, 5.0, 9.0])
gaussin(a, b)
print('*' * 20)
版权声明:本文为qq_36711003原创文章,遵循 CC 4.0 BY-SA 版权协议,转载请附上原文出处链接和本声明。