79 lines
2.1 KiB
Python
79 lines
2.1 KiB
Python
import math
|
||
|
||
# 列主元高斯消元法
|
||
def SovleRowMain(A,b):
|
||
ks = 0.00000001
|
||
n = len(A)
|
||
if len(A[0]) != n:
|
||
raise ValueError("A要为方阵")
|
||
if len(b) != n:
|
||
raise ValueError("b与A的行数不匹配")
|
||
p = list(range(n))
|
||
for i in range(n):
|
||
row_max = abs(A[i][i])
|
||
row_max_index = i
|
||
for j in range(i + 1, n):
|
||
if abs(A[j][i]) > row_max:
|
||
row_max = abs(A[j][i])
|
||
row_max_index = j
|
||
A[i], A[row_max_index] = A[row_max_index], A[i]
|
||
b[i], b[row_max_index] = b[row_max_index], b[i]
|
||
p[i], p[row_max_index] = p[row_max_index], p[i]
|
||
|
||
if abs(A[i][i]) < ks:
|
||
raise ValueError("A矩阵奇异，无法进行高斯消元")
|
||
for j in range(i + 1, n):
|
||
m = A[j][i] / A[i][i]
|
||
A[j][i] = m
|
||
for k in range(i + 1, n):
|
||
A[j][k] -= m * A[i][k]
|
||
b[j] -= m * b[i]
|
||
|
||
if abs(A[n - 1][n - 1]) < ks:
|
||
raise ValueError("A矩阵奇异，无法进行高斯消元")
|
||
|
||
# 回代求解
|
||
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] * b[j]
|
||
b[i] /= A[i][i]
|
||
return b
|
||
|
||
|
||
# 最小二乘法拟合
|
||
def LeastSquares(list_x,list_y,n):
|
||
m = len(list_x)
|
||
x_n = []
|
||
b = []
|
||
for i in range(2*n+1):
|
||
x_n.append(0)
|
||
b.append(0)
|
||
for j in range(m):
|
||
x_n[i]+=(list_x[j]**i)
|
||
b[i]+=(list_y[j]*list_x[j]**i)
|
||
b = b[:n+1]
|
||
A = []
|
||
for i in range(n+1):
|
||
tmp = []
|
||
for j in range(n+1):
|
||
tmp.append(x_n[i+j])
|
||
A.append(tmp)
|
||
print("A:", A)
|
||
print("b:", b)
|
||
result = SovleRowMain(A, b)
|
||
print("result:", result)
|
||
return result
|
||
|
||
#把x和y的值改为实际数据#######################
|
||
x = [1,2,4,8,16,32,64]
|
||
y = [4.22,4.02,3.85,3.59,3.44,3.02,2.59]
|
||
if __name__ == "__main__":
|
||
# 取对数 ln(W) = ln(C)+lamda*ln(t)
|
||
ln_W = [math.log(i) for i in y]
|
||
ln_t = [math.log(i) for i in x]
|
||
coeff = LeastSquares(ln_t,ln_W,1)
|
||
C = math.exp(coeff[0])
|
||
lamda = coeff[1]
|
||
print("C:", C)
|
||
print("lamda:", lamda) |