lhye200/CalWay_Python
1
1
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

all

Browse Source
This commit is contained in:
lhye200
2025-04-16 21:40:31 +08:00
commit bded38c54f
17 changed files with 1320 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
.gitignore vendored Normal file
View File

@@ -0,0 +1 @@
png/*

82
227-3.py Normal file
View File

@@ -0,0 +1,82 @@
A = [
[0, 3, 4],
[1, -1, 1],
[2, 1, 2]
]
b = [1, 2, 3]
# 列主元高斯消元法
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]
# 得到L,U和P矩阵
L = [[0 for i in range(n)] for j in range(n)]
U = [[0 for i in range(n)] for j in range(n)]
P = [[0 for i in range(n)] for j in range(n)]
for i in range(n):
for j in range(n):
if i == j:
L[i][j] = 1
U[i][j] = A[i][j]
elif i < j:
U[i][j] = A[i][j]
else:
L[i][j] = A[i][j]
P[i][p[i]] = 1
return P,L,U,b
def prettyPrintMatrix(matrix):
for row in matrix:
print(row)
if __name__ == "__main__":
# A = np.array(A, dtype=float)
# b = np.array(b, dtype=float)
P,L,U,x = SovleRowMain(A, b)
print("P:")
prettyPrintMatrix(P)
print("L:")
prettyPrintMatrix(L)
print("U:")
prettyPrintMatrix(U)
print("x:")
print(x)

45
227-4.py Normal file
View File

@@ -0,0 +1,45 @@
# 储存下三角矩阵 a11, a21, a22, a31, a32, a33 ...
A = [4,2,2,-2,-3,14]
b = [10,5,4]
#列 行
def getIndexFromDownMatrix(col, row):
if row > col:
row, col = col, row
return (1+col)*col//2+row
# 平方根法求解
def SqrtSolve(A,b):
n = len(b)
L = [0]*(n*(n+1)//2)
for j in range(n):
for i in range(j,n):
if i == j:
L[getIndexFromDownMatrix(i,j)] = A[getIndexFromDownMatrix(i,j)]
for k in range(j):
L[getIndexFromDownMatrix(i,j)] -= L[getIndexFromDownMatrix(j,k)]**2
L[getIndexFromDownMatrix(i,j)] = L[getIndexFromDownMatrix(i,j)]**0.5
else:
L[getIndexFromDownMatrix(i,j)] = A[getIndexFromDownMatrix(i,j)]
for k in range(j):
L[getIndexFromDownMatrix(i,j)] -= L[getIndexFromDownMatrix(i,k)]*L[getIndexFromDownMatrix(j,k)]
L[getIndexFromDownMatrix(i,j)] /= L[getIndexFromDownMatrix(j,j)]
for i in range(n):
for k in range(i):
b[i] -= L[getIndexFromDownMatrix(i,k)]*b[k]
b[i] /= L[getIndexFromDownMatrix(i,i)]
for i in range(n-1,-1,-1):
for k in range(i+1,n):
b[i] -= L[getIndexFromDownMatrix(k,i)]*b[k]
b[i] /= L[getIndexFromDownMatrix(i,i)]
return b
if __name__ == "__main__":
print("x:")
print(SqrtSolve(A,b))

37
227-7.py Normal file
View File

@@ -0,0 +1,37 @@
# 储存追赶法A矩阵
A = [
[0,4,-1],
[-1,4,-1],
[-1,4,-1],
[-1,4,-1],
[-1,4,0]
]
b = [100,200,200,200,100]
# 追赶法求解
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
if __name__ == "__main__":
print("x:")
print(ZGsolve(A,b))

77
68-1.py Normal file
View File

@@ -0,0 +1,77 @@
import math
list_x = [10,11,12,13]
list_y = [2.3026,2.3979,2.4849,2.5649]
# 定义原函数和其导函数计算结果
def FxDiff_n(x,n):
result = 0
if n == 0:
result = math.log(x)
else:
result = (-1)**(n+1) * math.factorial(n-1) / (x**n)
return result
# 获取与待求x最接近的两个点
def GetClosestTwo(x,list_x):
for i in range(0, len(list_x)):
if x < list_x[i]:
return i-1, i
return len(list_x)-2, len(list_x)-1
# 获取与待求x最接近的三个点
def GetClosestThree(x,list_x):
if x < list_x[1]:
return 0, 1, 2
for i in range(3, len(list_x)):
if x < list_x[i]:
return i-2, i-1, i
return len(list_x)-3, len(list_x)-2, len(list_x)-1
# 线性插值余项计算
def LinearRegression(x,list_x):
i, j = GetClosestTwo(x,list_x)
ks = max([abs(FxDiff_n(list_x[i],2)), abs(FxDiff_n(list_x[j],2))])
omg = (x-list_x[i])*(x-list_x[j])
return abs(ks*omg/2)
# 抛物线插值余项计算
def ParabolaRegression(x,list_x):
i, j, k = GetClosestThree(x,list_x)
ks = max([abs(FxDiff_n(list_x[i],3)), abs(FxDiff_n(list_x[j],3)), abs(FxDiff_n(list_x[k],3))])
omg = (x-list_x[i])*(x-list_x[j])*(x-list_x[k])
return abs(ks*omg/6)
# 线性插值
def LinearInterpolation(x,list_x,list_y):
i, j = GetClosestTwo(x,list_x)
result = list_y[i] + (x - list_x[i]) * (list_y[j] - list_y[i]) / (list_x[j] - list_x[i])
r = LinearRegression(x,list_x)
return (result,r)
# 抛物线插值
def ParabolaInterpolation(x,list_x,list_y):
i, j, k = GetClosestThree(x,list_x)
result = list_y[i] * (x - list_x[j]) * (x - list_x[k]) / (list_x[i] - list_x[j]) / (list_x[i] - list_x[k])
result += list_y[j] * (x - list_x[i]) * (x - list_x[k]) / (list_x[j] - list_x[i]) / (list_x[j] - list_x[k])
result += list_y[k] * (x - list_x[i]) * (x - list_x[j]) / (list_x[k] - list_x[i]) / (list_x[k] - list_x[j])
r = ParabolaRegression(x,list_x)
return (result,r)
# 拉格朗日插值
def LagrangeInterpolation(x,list_x,list_y):
result = 0
for i in range(0, len(list_x)):
temp = 1
for j in range(0, len(list_x)):
if i != j:
temp *= (x - list_x[j]) / (list_x[i] - list_x[j])
result += temp * list_y[i]
return result
if __name__ == "__main__":
print("线性插值 ln11.75 结果为%f, 截断误差%f" % LinearInterpolation(11.75,list_x,list_y))
print("抛物线插值 ln11.75 结果为%f, 截断误差%f" % ParabolaInterpolation(11.75,list_x,list_y))
# print(LagrangeInterpolation(11.75,list_x,list_y))

71
69-3.py Normal file
View File

@@ -0,0 +1,71 @@
import math
list_x = [0.0,0.2,0.4,0.6,0.8]
list_y = [1.0,1.2214,1.4918,1.8221,2.2255]
# 定义原函数和其导函数计算结果
def FxDiff_n(x,n):
return math.exp(x)
# 求差分表
def GetDyList(list_y):
result = []
result.append(list_y)
for i in range(len(list_y)-1,0,-1):
tmp = []
for j in range(i):
tmp.append(result[len(list_y)-1-i][j+1] - result[len(list_y)-1-i][j])
result.append(tmp)
return result
# 牛顿前插余项计算
def NewtonForwardRegression(x,n,list_x):
h = list_x[1] - list_x[0]
t = (x - list_x[0]) / h
ks = max([abs(FxDiff_n(list_x[i],n+1)) for i in range(n+1)])
result = 1
for i in range(n+1):
result *= (t-i)*h/(i+1)
return result * ks
# 牛顿前插
def NewtonForwardInterpolation(x,n,list_x,list_y):
list_dyk = GetDyList(list_y)
h = list_x[1] - list_x[0]
t = (x - list_x[0]) / h
result = list_y[0]
mul = 1
result = list_y[0]
for i in range(1, n):
mul *= ((t-i+1)/i)
result += list_dyk[i][0] * mul
r = abs(NewtonForwardRegression(x,n,list_x))
return (result,r)
# 求差商表
def GetDQyList(list_x,list_y):
result = []
result.append(list_y)
for i in range(len(list_y)-1,0,-1):
tmp = []
for j in range(i):
tmp.append((result[len(list_y)-1-i][j+1] - result[len(list_y)-1-i][j])/(list_x[j+len(list_y)-i] - list_x[j]))
result.append(tmp)
return result
# 牛顿基本插值
def NewtonBaseInterpolation(x,n,list_x,list_y):
list_dqy = GetDQyList(list_x,list_y)
result = list_dqy[0][0]
mul = 1
for i in range(1,n):
mul *= (x - list_x[i-1])
result += list_dqy[i][0] * mul
return result
if __name__ == '__main__':
print("三点牛顿前插 e^0.12 结果为%f, 截断误差%f" % NewtonForwardInterpolation(0.12, 2, list_x, list_y))
print("四点牛顿前插 e^0.12 结果为%f, 截断误差%f" % NewtonForwardInterpolation(0.12, 3, list_x, list_y))
print("三点牛顿基本插值 e^0.12 结果为%f" % NewtonBaseInterpolation(0.12, 2, list_x, list_y))
print("四点牛顿基本插值 e^0.12 结果为%f" % NewtonBaseInterpolation(0.12, 3, list_x, list_y))

92
69-5 - 副本.py Normal file
View File

@@ -0,0 +1,92 @@
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))

91
69-5.py Normal file
View File

@@ -0,0 +1,91 @@
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
def PrintResult(M,list_h,list_x,list_y):
for i in range(len(list_h)):
k1 = (M[i+1]-M[i])/6/list_h[i]
k2 = (M[i]*list_x[i+1]-M[i+1]*list_x[i])/2/list_h[i]
k3 = (3*M[i+1]*list_x[i]**2-3*M[i]*list_x[i+1]**2-6*list_y[i]+M[i]*list_h[i]**2+6*list_y[i+1]-M[i+1]*list_h[i]**2)/6/list_h[i]
k4 = (M[i]*list_x[i+1]**3-M[i+1]*list_x[i]**3+6*list_y[i]*list_x[i+1]-M[i]*list_h[i]**2*list_x[i+1]-6*list_y[i+1]*list_x[i]+M[i+1]*list_h[i]**2*list_x[i])/6/list_h[i]
print("S(x)=%.6f*x^3+%.6f*x^2+%6f*x+%6f"%(k1,k2,k3,k4),"x=[%.6f,%.6f]"%(list_x[i],list_x[i+1]))
if __name__ == "__main__":
M,list_h = CubicSplineInterpolation(x,y,2,1,0)
print("二阶导数边界条件:")
PrintResult(M,list_h,x,y)
M,list_h = CubicSplineInterpolation(x,y,1,1,0)
print("一阶导数边界条件:")
PrintResult(M,list_h,x,y)

116
89-1 - 副本.py Normal file
View File

@@ -0,0 +1,116 @@
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))

116
89-1.py Normal file
View File

@@ -0,0 +1,116 @@
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 Normal file
View File

@@ -0,0 +1,76 @@
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)

76
89-2.py Normal file
View File

@@ -0,0 +1,76 @@
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 Normal file
View File

@@ -0,0 +1,69 @@
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]))

69
89-3.py Normal file
View File

@@ -0,0 +1,69 @@
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]))

7
test.py Normal file
View File

@@ -0,0 +1,7 @@
from sympy import diff
from sympy import symbols
def func(x):
return x**4
x = symbols("x")
print(diff(func(x),x))

163
注释完整版/68-1-zhushi.py Normal file
View File

@@ -0,0 +1,163 @@
import math
# 定义插值点
list_x = [10, 11, 12, 13] # 自变量列表
list_y = [2.3026, 2.3979, 2.4849, 2.5649] # 因变量列表
# 定义原函数和其导函数计算结果
def FxDiff_n(x, n):
"""
计算自然对数函数 ln(x) 的 n 阶导数值。
参数:
x (float): 自变量值，要求 x > 0。
n (int): 导数的阶数，n >= 0。
返回:
float: ln(x) 的 n 阶导数值。
"""
result = 0
if n == 0:
result = math.log(x) # 原函数 ln(x)
else:
result = (-1)**(n+1) * math.factorial(n-1) / (x**n) # n 阶导数公式
return result
# 获取与待求 x 最接近的两个点
def GetClosestTwo(x, list_x):
"""
获取与待求点 x 最接近的两个插值点的索引。
参数:
x (float): 待插值的自变量值。
list_x (list of float): 自变量列表，要求已按升序排列。
返回:
tuple of int: 最接近的两个点的索引 (i, j)，满足 list_x[i] <= x < list_x[j]。
"""
for i in range(0, len(list_x)):
if x < list_x[i]:
return i-1, i
return len(list_x)-2, len(list_x)-1
# 获取与待求 x 最接近的三个点
def GetClosestThree(x, list_x):
"""
获取与待求点 x 最接近的三个插值点的索引。
参数:
x (float): 待插值的自变量值。
list_x (list of float): 自变量列表，要求已按升序排列。
返回:
tuple of int: 最接近的三个点的索引 (i, j, k)，满足 list_x[i] <= x < list_x[k]。
"""
if x < list_x[1]:
return 0, 1, 2
for i in range(3, len(list_x)):
if x < list_x[i]:
return i-2, i-1, i
return len(list_x)-3, len(list_x)-2, len(list_x)-1
# 线性插值余项计算
def LinearRegression(x, list_x):
"""
计算线性插值的余项估计值。
参数:
x (float): 待插值的自变量值。
list_x (list of float): 自变量列表，要求已按升序排列。
返回:
float: 线性插值的余项估计值。
"""
i, j = GetClosestTwo(x, list_x)
ks = max([abs(FxDiff_n(list_x[i], 2)), abs(FxDiff_n(list_x[j], 2))]) # 二阶导数的最大值
omg = (x - list_x[i]) * (x - list_x[j]) # 乘积项
return abs(ks * omg / 2)
# 抛物线插值余项计算
def ParabolaRegression(x, list_x):
"""
计算抛物线插值的余项估计值。
参数:
x (float): 待插值的自变量值。
list_x (list of float): 自变量列表，要求已按升序排列。
返回:
float: 抛物线插值的余项估计值。
"""
i, j, k = GetClosestThree(x, list_x)
ks = max([abs(FxDiff_n(list_x[i], 3)), abs(FxDiff_n(list_x[j], 3)), abs(FxDiff_n(list_x[k], 3))]) # 三阶导数的最大值
omg = (x - list_x[i]) * (x - list_x[j]) * (x - list_x[k]) # 乘积项
return abs(ks * omg / 6)
# 线性插值
def LinearInterpolation(x, list_x, list_y):
"""
使用线性插值法计算插值结果及其余项估计值。
参数:
x (float): 待插值的自变量值。
list_x (list of float): 自变量列表，要求已按升序排列。
list_y (list of float): 因变量列表，与 list_x 一一对应。
返回:
tuple: 插值结果和余项估计值 (result, r)。
"""
i, j = GetClosestTwo(x, list_x)
result = list_y[i] + (x - list_x[i]) * (list_y[j] - list_y[i]) / (list_x[j] - list_x[i]) # 插值公式
r = LinearRegression(x, list_x) # 余项估计
return (result, r)
# 抛物线插值
def ParabolaInterpolation(x, list_x, list_y):
"""
使用抛物线插值法计算插值结果及其余项估计值。
参数:
x (float): 待插值的自变量值。
list_x (list of float): 自变量列表，要求已按升序排列。
list_y (list of float): 因变量列表，与 list_x 一一对应。
返回:
tuple: 插值结果和余项估计值 (result, r)。
"""
i, j, k = GetClosestThree(x, list_x)
result = list_y[i] * (x - list_x[j]) * (x - list_x[k]) / (list_x[i] - list_x[j]) / (list_x[i] - list_x[k])
result += list_y[j] * (x - list_x[i]) * (x - list_x[k]) / (list_x[j] - list_x[i]) / (list_x[j] - list_x[k])
result += list_y[k] * (x - list_x[i]) * (x - list_x[j]) / (list_x[k] - list_x[i]) / (list_x[k] - list_x[j])
r = ParabolaRegression(x, list_x) # 余项估计
return (result, r)
# 拉格朗日插值
def LagrangeInterpolation(x, list_x, list_y):
"""
使用拉格朗日插值法计算插值结果。
参数:
x (float): 待插值的自变量值。
list_x (list of float): 自变量列表，要求已按升序排列。
list_y (list of float): 因变量列表，与 list_x 一一对应。
返回:
float: 插值结果。
"""
result = 0
for i in range(0, len(list_x)):
temp = 1
for j in range(0, len(list_x)):
if i != j:
temp *= (x - list_x[j]) / (list_x[i] - list_x[j]) # 拉格朗日基函数
result += temp * list_y[i]
return result
if __name__ == "__main__":
# 测试线性插值
print(LinearInterpolation(11.75, list_x, list_y))
# 测试抛物线插值
print(ParabolaInterpolation(11.75, list_x, list_y))
# 测试拉格朗日插值
print(LagrangeInterpolation(11.75, list_x, list_y))

132
注释完整版/69-3-zhushi.py Normal file
View File

@@ -0,0 +1,132 @@
import math
import numpy as np
# 定义插值点
list_x = [0.0, 0.2, 0.4, 0.6, 0.8] # 自变量列表
list_y = [1.0, 1.2214, 1.4918, 1.8221, 2.2255] # 因变量列表
def GetDyList(list_y):
"""
计算前向差分表。
参数:
list_y (list of float): 因变量列表。
返回:
list of list of float: 前向差分表，每一行是对应阶数的差分值。
"""
result = []
result.append(list_y)
for i in range(len(list_y) - 1, 0, -1):
tmp = []
for j in range(i):
tmp.append(result[len(list_y) - 1 - i][j + 1] - result[len(list_y) - 1 - i][j])
result.append(tmp)
return result
def FxDiff_n(x, n):
"""
计算函数 e^x 的 n 阶导数值。
参数:
x (float): 自变量值。
n (int): 导数的阶数，n >= 0。
返回:
float: e^x 的 n 阶导数值。
"""
return math.exp(x)
def NewtonForwardRegression(x, n, list_x):
"""
计算牛顿前向插值的余项估计值。
参数:
x (float): 待插值的自变量值。
n (int): 插值多项式的阶数。
list_x (list of float): 自变量列表，要求等间距。
返回:
float: 牛顿前向插值的余项估计值。
"""
h = list_x[1] - list_x[0] # 等间距
t = (x - list_x[0]) / h # 归一化变量
ks = max([abs(FxDiff_n(list_x[i], n + 1)) for i in range(n + 1)]) # n+1 阶导数的最大值
result = 1
for i in range(n + 1):
result *= (t - i) * h / (i + 1) # 余项公式中的乘积项
return result * ks
def NewtonForwardInterpolation(x, n, list_x, list_y):
"""
使用牛顿前向插值法计算插值结果及其余项估计值。
参数:
x (float): 待插值的自变量值。
n (int): 插值多项式的阶数。
list_x (list of float): 自变量列表，要求等间距。
list_y (list of float): 因变量列表，与 list_x 一一对应。
返回:
tuple: 插值结果和余项估计值 (result, r)。
"""
list_dyk = GetDyList(list_y) # 差分表
h = list_x[1] - list_x[0] # 等间距
t = (x - list_x[0]) / h # 归一化变量
result = list_y[0] # 插值结果初值
mul = 1 # 乘积项
for i in range(1, n):
mul *= ((t - i + 1) / i) # 计算乘积项
result += list_dyk[i][0] * mul # 累加差分项
r = abs(NewtonForwardRegression(x, n, list_x)) # 余项估计
return (result, r)
def GetDQyList(list_x, list_y):
"""
计算牛顿基函数插值的差商表。
参数:
list_x (list of float): 自变量列表。
list_y (list of float): 因变量列表。
返回:
list of list of float: 差商表，每一行是对应阶数的差商值。
"""
result = []
result.append(list_y)
for i in range(len(list_y) - 1, 0, -1):
tmp = []
for j in range(i):
tmp.append((result[len(list_y) - 1 - i][j + 1] - result[len(list_y) - 1 - i][j]) /
(list_x[j + len(list_y) - i] - list_x[j]))
result.append(tmp)
return result
def NewtonBaseInterpolation(x, n, list_x, list_y):
"""
使用牛顿基函数插值法计算插值结果。
参数:
x (float): 待插值的自变量值。
n (int): 插值多项式的阶数。
list_x (list of float): 自变量列表。
list_y (list of float): 因变量列表，与 list_x 一一对应。
返回:
float: 插值结果。
"""
list_dqy = GetDQyList(list_x, list_y) # 差商表
result = list_dqy[0][0] # 插值结果初值
mul = 1 # 乘积项
for i in range(1, n):
mul *= (x - list_x[i - 1]) # 计算乘积项
result += list_dqy[i][0] * mul # 累加差商项
return result
if __name__ == '__main__':
# 测试牛顿前向插值
print(NewtonForwardInterpolation(0.12, 2, list_x, list_y))
print(NewtonForwardInterpolation(0.12, 3, list_x, list_y))
# 测试牛顿基函数插值
print(NewtonBaseInterpolation(0.12, 2, list_x, list_y))
print(NewtonBaseInterpolation(0.12, 3, list_x, list_y))