"""
本函数通过创建增广矩阵,并调用高斯列主元消去法模块进行求解。
"""import numpy as np
import gaosi as gs
shape =int(input('请输入拟合函数的次数:'))
x = np.array([0.6,1.3,1.64,1.8,2.1,2.3,2.44])
y = np.array([7.05,12.2,14.4,15.2,17.4,19.6,20.2])
data =[]for i inrange(shape*2+1):if i !=0:
data.append(np.sum(x**i))else:
data.append(len(x))
b =[]for i inrange(shape+1):if i !=0:
b.append(np.sum(y*x**i))else:
b.append(np.sum(y))
b = np.array(b).reshape(shape+1,1)
n = np.zeros([shape+1,shape+1])for i inrange(shape+1):for j inrange(shape+1):
n[i][j]= data[i+j]
result = gs.Handle(n,b)ifnot result:print('增广矩阵求解失败!')
exit()
fun='f(x) = 'for i inrange(len(result)):iftype(result[i])==type(''):print('存在自由变量!')
fun = fun +str(result[i])elif i ==0:
fun = fun +'{:.3f}'.format(result[i])else:
fun = fun +'+{0:.3f}*x^{1}'.format(result[i],i)print('求得{0}次拟合函数为:'.format(shape))print(fun)
高斯模块
# 导入 numpy 模块import numpy as np
# 行交换defswap_row(matrix, i, j):
m, n = matrix.shape
if i >= m or j >= m:print('错误! : 行交换超出范围 ...')else:
matrix[i],matrix[j]= matrix[j].copy(),matrix[i].copy()return matrix
# 变成阶梯矩阵defmatrix_change(matrix):
m, n = matrix.shape
main_factor =[]
main_col = main_row =0while main_row < m and main_col < n:# 选择进行下一次主元查找的列
main_row =len(main_factor)# 寻找列中非零的元素
not_zeros = np.where(abs(matrix[main_row:,main_col])>0)[0]# 如果该列向下全部数据为零,则直接跳过列iflen(not_zeros)==0:
main_col +=1continueelse:# 将主元列号保存在列表中
main_factor.append(main_col)# 将第一个非零行交换至最前if not_zeros[0]!=[0]:
matrix = swap_row(matrix,main_row,main_row+not_zeros[0])# 将该列主元下方所有元素变为零if main_row < m-1:for k inrange(main_row+1,m):
a =float(matrix[k, main_col]/ matrix[main_row, main_col])
matrix[k]= matrix[k]- matrix[main_row]* matrix[k, main_col]/ matrix[main_row, main_col]
main_col +=1return matrix,main_factor
# 回代求解defback_solve(matrix, main_factor):# 判断是否有解iflen(main_factor)==0:print('主元错误,无主元! ...')returnNone
m, n = matrix.shape
if main_factor[-1]== n -1:print('无解! ...')returnNone# 把所有的主元元素上方的元素变成0for i inrange(len(main_factor)-1,-1,-1):
factor = matrix[i, main_factor[i]]
matrix[i]= matrix[i]/float(factor)for j inrange(i):
times = matrix[j, main_factor[i]]
matrix[j]= matrix[j]-float(times)* matrix[i]# 先看看结果对不对return matrix
# 结果打印defprint_result(matrix, main_factor):if matrix isNone:print('阶梯矩阵为空! ...')returnNone
m, n = matrix.shape
result =['']*(n -1)
main_factor =list(main_factor)for i inrange(n -1):# 如果不是主元列,则为自由变量if i notin main_factor:
result[i]='(free var)'# 否则是主元变量,从对应的行,将主元变量表示成非主元变量的线性组合else:# row_of_main表示该主元所在的行
row_of_main = main_factor.index(i)
result[i]= matrix[row_of_main,-1]return result
# 得到简化的阶梯矩阵和主元列defHandle(matrix_a, matrix_b):# 拼接成增广矩阵
matrix_01 = np.hstack([matrix_a, matrix_b])
matrix_01, main_factor = matrix_change(matrix_01)
matrix_01 = back_solve(matrix_01, main_factor)
result = print_result(matrix_01, main_factor)return result
if __name__ =='__main__':
a = np.array([[2,1,1],[3,1,2],[1,2,2]],dtype=float)
b = np.array([[4],[6],[5]],dtype=float)
a = Handle(a, b)
模块导入import numpy as npimport gaosi as gs代码"""本函数通过创建增广矩阵,并调用高斯列主元消去法模块进行求解。"""import numpy as npimport gaosi as gsshape = int(input('请输入拟合函数的次数:'))x = np.array([0.6,1.3,1.64,1.8,2.1,2.3,2.44])y = np.array([7.05,12.2,14.4,15.2,17.4,19.6,20.2]