Report 2
为顺序消元法和列主元消元法我分别写了两个函数,即gauss_se
(gauss sequntial elimination) 与 gauss_cpe
(gauss column-pivoting elimination) 。它们的区别只在于,在 gauss_se
中,我每次循环开始检查一下 a[i, i]
是否为 0 ,如果为 0 ,则将其下不为 0 的行与当前行交换,如果下面的对应元素也都为 0 ,那就跳出。而在 gauss_cpe
里,每次循环开始不会检查,而是直接将 a[i, i]
下绝对值最大的那行与当前行交换,当然,如果这个绝对值最大也才是 0 ,那么就跳出。
import numpy as np
from numpy import dot
def gauss_se(a,b):
a = np.copy(a)
b = np.copy(b)
n = len(b)
for i in range(n-1):
if a[i, i] == 0:
ind = aux_func(a[i:, i+1])
if ind is None:
break
a[i,:], a[i+ind+1, :] = a[i+ind+1, :], a[i,:]
b[i], b[i+ind+1] = b[i+ind+1], b[i]
for j in range(i+1,n):
if a[j,i] != 0:
lam = a[j,i] / a[i,i]
a[j,:] = a[j,:] - lam*a[i,:]
b[j] = b[j] - lam*b[i]
for k in range(n-1,-1,-1):
b[k] = (b[k] - dot(a[k,(k+1):n],b[(k+1):n]))/a[k,k]
result = b
return result
def gauss_cpe(a, b):
a = np.copy(a)
b = np.copy(b)
n = len(b)
for i in range(n-1):
ind = np.argmax(np.abs(a[i:, i]))
if a[i, i+ind]==0:
break
a[[i, i+ind],:]= a[[i+ind, i],:]
b[[i, i+ind]] = b[[i+ind, i]]
for j in range(i+1,n):
if a[j,i] != 0:
lam = a[j,i] / a[i,i]
a[j,:] = a[j,:] - lam*a[i,:]
b[j] = b[j] - lam*b[i]
for k in range(n-1,-1,-1):
b[k] = (b[k] - dot(a[k,(k+1):n],b[(k+1):n]))/a[k,k]
result = b
return result
def aux_func_1(x):
for i,t in enumerate(x):
if t!=0:
return i
return None
def main():
A = [[1.1348, 3.8326, 1.1651, 3.4017], [0.5301, 1.7875, 2.5330, 1.5435],
[3.4129, 4.9317, 8.7643, 1.3142], [1.2371, 4.9998, 10.6721, 0.0147]]
b = [9.5342, 6.3941, 18.4231, 16.9237]
A = np.array(A)
b = np.array(b).T
print(gauss_se(A, b))
print(gauss_cpe(A, b))
if __name__ == '__main__':
main()
输出如下:
[ 1. 1. 1. 1.]
[ 1. 1. 1. 1.]