111
This commit is contained in:
92
69-5 - 副本.py
92
69-5 - 副本.py
@@ -1,92 +0,0 @@
|
|||||||
x = [0,1,2,3]
|
|
||||||
y = [0,0,0,0]
|
|
||||||
|
|
||||||
def ZGsolve(A,b):
|
|
||||||
n = len(b)
|
|
||||||
beta = [0]*n
|
|
||||||
for i in range(n):
|
|
||||||
if i == 0:
|
|
||||||
beta[i] = A[i][2] / A[i][1]
|
|
||||||
else:
|
|
||||||
beta[i] = A[i][2] / (A[i][1] - A[i][0]*beta[i-1])
|
|
||||||
|
|
||||||
for i in range(n):
|
|
||||||
if i == 0:
|
|
||||||
b[i] = b[i] / A[i][1]
|
|
||||||
else:
|
|
||||||
b[i] = (b[i] - A[i][0]*b[i-1]) / (A[i][1] - A[i][0]*beta[i-1])
|
|
||||||
|
|
||||||
for i in range(n-2,-1,-1):
|
|
||||||
b[i] = b[i] - beta[i]*b[i+1]
|
|
||||||
|
|
||||||
return b
|
|
||||||
|
|
||||||
|
|
||||||
def GetDList(list_r):
|
|
||||||
result = []
|
|
||||||
for i in range(1,len(list_r)):
|
|
||||||
result.append(list_r[i] - list_r[i-1])
|
|
||||||
return result
|
|
||||||
|
|
||||||
def GetDQList(list_x, list_y):
|
|
||||||
result = []
|
|
||||||
for i in range(1,len(list_y)):
|
|
||||||
result.append((list_y[i] - list_y[i-1]) / (list_x[i] - list_x[i-1]))
|
|
||||||
return result
|
|
||||||
|
|
||||||
def CubicSplineInterpolation(list_x,list_y,boundary_type,a1,a2):
|
|
||||||
list_h = GetDList(list_x)
|
|
||||||
list_dqxy = GetDQList(list_x, list_y)
|
|
||||||
list_mu = [list_h[i]/(list_h[i]+list_h[i+1]) for i in range(len(list_h)-1)]
|
|
||||||
list_lamda = [1-i for i in list_mu]
|
|
||||||
list_g = [6*(list_dqxy[i+1]-list_dqxy[i])/(list_h[i+1]+list_h[i]) for i in range(len(list_h)-1)]
|
|
||||||
|
|
||||||
A = []
|
|
||||||
b = []
|
|
||||||
if boundary_type == 0: # 自然边界条件
|
|
||||||
pass
|
|
||||||
elif boundary_type == 1: # 一阶导数边界条件
|
|
||||||
A.append([0,2,1])
|
|
||||||
b.append(6/list_h[0]*(list_dqxy[1]-a1))
|
|
||||||
for i in range(len(list_g)):
|
|
||||||
A.append([list_mu[i],2,list_lamda[i]])
|
|
||||||
b.append(list_g[i])
|
|
||||||
A.append([1,2,0])
|
|
||||||
b.append(6/list_h[-1]*(a2-list_dqxy[-1]))
|
|
||||||
M = ZGsolve(A,b)
|
|
||||||
|
|
||||||
elif boundary_type == 2: # 二阶导数边界条件
|
|
||||||
A.append([0,2,list_lamda[0]])
|
|
||||||
b.append(list_g[0]-list_mu[0]*a1)
|
|
||||||
for i in range(1,len(list_g)-1):
|
|
||||||
A.append([list_mu[i],2,list_lamda[i]])
|
|
||||||
b.append(list_g[i])
|
|
||||||
A.append([list_mu[-1],2,0])
|
|
||||||
b.append(list_g[-1]-list_lamda[-1]*a2)
|
|
||||||
M = ZGsolve(A,b)
|
|
||||||
M = [a1] + M + [a2]
|
|
||||||
|
|
||||||
return M,list_h
|
|
||||||
|
|
||||||
|
|
||||||
if __name__ == "__main__":
|
|
||||||
M,list_h = CubicSplineInterpolation(x,y,2,1,0)
|
|
||||||
|
|
||||||
print("list_h:",list_h)
|
|
||||||
print("M:",M)
|
|
||||||
for i in range(3):
|
|
||||||
k1 = (M[i+1]-M[i])/6/list_h[i]
|
|
||||||
k2 = (M[i]*x[i+1]-M[i+1]*x[i])/2/list_h[i]
|
|
||||||
k3 = (3*M[i+1]*x[i]**2-3*M[i]*x[i+1]**2-6*y[i]+M[i]*list_h[i]**2+6*y[i+1]-M[i+1]*list_h[i]**2)/6/list_h[i]
|
|
||||||
k4 = (M[i]*x[i+1]**3-M[i+1]*x[i]**3+6*y[i]*x[i+1]-M[i]*list_h[i]**2*x[i+1]-6*y[i+1]*x[i]+M[i+1]*list_h[i]**2*x[i])/6/list_h[i]
|
|
||||||
# print(k1,k2,k3,k4)
|
|
||||||
print("S(x)=%.10f*x^3+%.10f*x^2+%.10f*x+%.10f"%(k1,k2,k3,k4))
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
116
89-1 - 副本.py
116
89-1 - 副本.py
@@ -1,116 +0,0 @@
|
|||||||
x = [2,4,6,8]
|
|
||||||
y = [2,11,28,40]
|
|
||||||
|
|
||||||
|
|
||||||
# 列主元高斯消元法
|
|
||||||
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)
|
|
||||||
return SovleRowMain(A,b)
|
|
||||||
|
|
||||||
def CalculateY(list_x, coeff):
|
|
||||||
re = []
|
|
||||||
for i in range(len(list_x)):
|
|
||||||
re.append(0)
|
|
||||||
for j in range(len(coeff)):
|
|
||||||
re[i] += coeff[j]*list_x[i]**j
|
|
||||||
return re
|
|
||||||
|
|
||||||
def MeanSquareErr(list_y,list_y_approx):
|
|
||||||
m = len(list_y)
|
|
||||||
err = 0
|
|
||||||
for i in range(m):
|
|
||||||
err += (list_y[i] - list_y_approx[i])**2
|
|
||||||
return err**0.5
|
|
||||||
|
|
||||||
def MaxErr(list_y,list_y_approx):
|
|
||||||
m = len(list_y)
|
|
||||||
err = 0
|
|
||||||
for i in range(m):
|
|
||||||
if abs(list_y[i] - list_y_approx[i]) > err:
|
|
||||||
err = abs(list_y[i] - list_y_approx[i])
|
|
||||||
return err
|
|
||||||
|
|
||||||
def PrintEquation(coeff):
|
|
||||||
n = len(coeff)
|
|
||||||
str_ = str(coeff[0]) + "+"
|
|
||||||
for i in range(0,n-1):
|
|
||||||
str_ += str(coeff[i+1])
|
|
||||||
str_ += "*x^"+str(i+1)+"+"
|
|
||||||
str_ = str_[0:len(str_)-1]
|
|
||||||
print(str_)
|
|
||||||
|
|
||||||
if __name__ == "__main__":
|
|
||||||
print("一次拟合")
|
|
||||||
coeff = LeastSquares(x,y,1)
|
|
||||||
PrintEquation(coeff)
|
|
||||||
y_approx = CalculateY(x, coeff)
|
|
||||||
print("MeanSquareErr:")
|
|
||||||
print(MeanSquareErr(y,y_approx))
|
|
||||||
print("MaxErr:")
|
|
||||||
print(MaxErr(y,y_approx))
|
|
||||||
|
|
||||||
print("二次拟合")
|
|
||||||
coeff = LeastSquares(x,y,2)
|
|
||||||
PrintEquation(coeff)
|
|
||||||
y_approx = CalculateY(x, coeff)
|
|
||||||
print("MeanSquareErr:")
|
|
||||||
print(MeanSquareErr(y,y_approx))
|
|
||||||
print("MaxErr:")
|
|
||||||
print(MaxErr(y,y_approx))
|
|
||||||
|
|
||||||
76
89-2 - 副本.py
76
89-2 - 副本.py
@@ -1,76 +0,0 @@
|
|||||||
import math
|
|
||||||
|
|
||||||
x = [1,2,4,8,16,32,64]
|
|
||||||
y = [4.22,4.02,3.85,3.59,3.44,3.02,2.59]
|
|
||||||
|
|
||||||
|
|
||||||
# 列主元高斯消元法
|
|
||||||
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)
|
|
||||||
return SovleRowMain(A,b)
|
|
||||||
|
|
||||||
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)
|
|
||||||
69
89-3 - 副本.py
69
89-3 - 副本.py
@@ -1,69 +0,0 @@
|
|||||||
|
|
||||||
x = [19,25,31,38,44]
|
|
||||||
y = [19.0,32.3,49.0,73.3,97.8]
|
|
||||||
|
|
||||||
# 列主元高斯消元法
|
|
||||||
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)
|
|
||||||
return SovleRowMain(A,b)
|
|
||||||
|
|
||||||
if __name__ == "__main__":
|
|
||||||
x_square = [i**2 for i in x]
|
|
||||||
coeff = LeastSquares(x_square, y, 1)
|
|
||||||
print("拟合方程: y = %.6f + %.6f*x^2" % (coeff[0], coeff[1]))
|
|
||||||
Reference in New Issue
Block a user