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

aaa

Browse Source
This commit is contained in:
lhye200
2025-06-16 20:44:29 +08:00
parent 3a41897ba2
commit 6baaac12c5
17 changed files with 1100 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

54
按方法整理/常微分方程-阿达姆斯方法.py Normal file
View File

@@ -0,0 +1,54 @@
import math
def AdamusExplicitly(k,first_ys,x0,y0,h,xk,fxy,fx_real):
betas = [
[1],
[3/2, -1/2],
[23/12, -16/12, 5/12],
[55/24, -59/24, 37/24, -9/24],
[1901/720, -2774/720, 2626/720, -1274/720, 251/720],
[4277/1440, -7923/1600, 9982/1440, -7298/1440, 2877/1440, -476/1440]
]
# Bks = [1/2,5/12,3/8,251/720,95/288,10987/60480]
result = []
if k<0 or k >= len(betas):
print("k超出范围")
return None
if len(first_ys) != k + 1:
print("first_ys(前几个y的值)与k长度不匹配")
return None
fxys = [fxy(x0+i*h, first_ys[i]) for i in range(k + 1)]
for i in range(k + 1):
result.append((x0 + i * h, first_ys[i]))
x1 = x0 + h*k
y0 = first_ys[k]
while x1 < xk:
x1 += h
delta_y = h * sum(betas[k][i] * fxys[k-i] for i in range(k + 1))
y1 = y0 + delta_y
print(f"x={x1}, y={y1}, delta_y={delta_y}, real_y={fx_real(x1)}, abs(real_y-y)={abs(fx_real(x1) - y1)}")
result.append((x1, y1))
fxys.append(fxy(x1, y1))
fxys.pop(0)
y0 = y1
return result
if __name__ == "__main__":
# 定义初始值和参数#####################################################################################
k = 3 # 阿达姆斯k+1步显式方法，精度为k+1阶，最常用k=3，其他阶数我没试过 P253
first_ys = [1, 0.904837418036, 0.818730753078, 0.7408182206817] # 前几个y的值，可用龙格-库塔计算或者知道精确解自己算出来
x0 = 0.0 # 初始x值
y0 = 1.0 # 初始y值
h = 0.1 # 步长
xk = 1.0 # 最终x值
fxy = lambda x, y: -y # 定义f(x,y)函数，dx/dy = f(x,y)，导函数
fx_real = lambda x: math.exp(-x) # 实际解函数，用于验证结果，如果不知道或者不用算误差，可以直接写个 lambda x: 0
# 调用阿达姆斯显式方法
result = AdamusExplicitly(k, first_ys, x0, y0, h, xk, fxy, fx_real)
print("计算结果看到有几.99999或者几.00000就自己四舍五入一下，有可能会多算一点，自己比较一下")
if result:
for xy in result:
print(f"(x, y): {xy}")

65
按方法整理/常微分方程-阿达姆斯预测-校正方法.py Normal file
View File

@@ -0,0 +1,65 @@
import math
def AdamusFix(k,first_ys,x0,y0,h,xk,fxy,fx_real):
betas_1 = [
[1],
[3/2, -1/2],
[23/12, -16/12, 5/12],
[55/24, -59/24, 37/24, -9/24],
[1901/720, -2774/720, 2626/720, -1274/720, 251/720],
[4277/1440, -7923/1600, 9982/1440, -7298/1440, 2877/1440, -476/1440]
]
# Bks = [1/2,5/12,3/8,251/720,95/288,10987/60480]
betas_2 = [
[1],
[1/2, 1/2],
[5/12, 8/12, -1/12],
[9/24, 19/24, -5/24, 1/24],
[251/720, 646/720, -261/720, 106/720, -19/720],
[475/1440, 1427/1440, -798/1440, 482/1440, -173/1440, 27/1440]
]
result = []
if k<0 or k >= len(betas_1):
print("k超出范围")
return None
if len(first_ys) != k + 1:
print("first_ys(前几个y的值)与k长度不匹配")
return None
fxys = [fxy(x0+i*h, first_ys[i]) for i in range(k + 1)]
for i in range(k + 1):
result.append((x0 + i * h, first_ys[i]))
x1 = x0 + h*k
y0 = first_ys[k]
while x1 < xk:
x1 += h
y1_guess = y0 + h * sum(betas_1[k][i] * fxys[k-i] for i in range(k + 1))
fxy_guess = fxy(x1, y1_guess)
fxys.append(fxy_guess)
fxys.pop(0)
delta_y = h * sum(betas_2[k][i] * fxys[k-i] for i in range(k + 1))
y1 = y0 + delta_y
print(f"x={x1}, y={y1}, y_guess={y1_guess}, real_y={fx_real(x1)}, abs(real_y-y)={abs(fx_real(x1) - y1)}")
result.append((x1, y1))
fxys[k] = fxy(x1, y1)
y0 = y1
return result
if __name__ == "__main__":
# 定义初始值和参数#####################################################################################
k = 3 # 阿达姆斯k+1步显式方法，精度为k+1阶，最常用k=3，其他阶数我没试过 P253
first_ys = [1, 1.0954, 1.1832, 1.2649] # 前几个y的值，可用龙格-库塔计算或者知道精确解自己算出来
x0 = 0.0 # 初始x值
y0 = 1.0 # 初始y值
h = 0.1 # 步长
xk = 1.0 # 最终x值
fxy = lambda x, y: y - 2*x/y # 定义f(x,y)函数，dx/dy = f(x,y)，导函数
fx_real = lambda x: 0 # 实际解函数，用于验证结果，如果不知道或者不用算误差，可以直接写个 lambda x: 0
# 调用阿达姆斯显式方法
result = AdamusFix(k, first_ys, x0, y0, h, xk, fxy, fx_real)
print("计算结果看到有几.99999或者几.00000就自己四舍五入一下，有可能会多算一点，自己比较一下")
if result:
for xy in result:
print(f"(x, y): {xy}")

133
按方法整理/拟合-最小二乘法.py Normal file
View File

@@ -0,0 +1,133 @@
import math
# 列主元高斯消元法
def SovleRowMain(A,b):
ks = 0.00000001
n = len(A)
if len(A[0]) != n:
print("A要为方阵")
return None,None,None,None
if len(b) != n:
print("b与A的行数不匹配")
return None,None,None,None
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:
print("A矩阵奇异，无法进行高斯消元")
return None,None,None,None
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:
print("A矩阵奇异，无法进行高斯消元")
return None,None,None,None
# 回代求解
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值上的值
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__":
##########################################################################
# 下面的数值根据题意改成适合的，比如要取对数就先取对数，下面的x和y是直接用于拟合的
list_x = [2,4,6,8] # 已给出的x数值，与y数值对应
list_y = [2,11,28,40] # 已给出的y数值，与x数值对应
# 记得改下面最后一个参数，为拟合阶数
print("一次拟合")
coeff = LeastSquares(list_x,list_y,1)
PrintEquation(coeff)
y_approx = CalculateY(list_x, coeff)
print("均方根误差:")
print(MeanSquareErr(list_y,y_approx))
print("最大误差:")
print(MaxErr(list_y,y_approx))
print("二次拟合")
coeff = LeastSquares(list_x,list_y,2)
PrintEquation(coeff)
y_approx = CalculateY(list_x, coeff)
print("均方根误差:")
print(MeanSquareErr(list_y,y_approx))
print("最大误差:")
print(MaxErr(list_y,y_approx))

122
按方法整理/插值-三次样条.py Normal file
View File

@@ -0,0 +1,122 @@
import math
# 追赶法
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)
print("h:", list_h)
list_dqxy = GetDQList(list_x, list_y)
print("f[xi,xi+1]:", list_dqxy)
list_mu = [list_h[i]/(list_h[i]+list_h[i+1]) for i in range(len(list_h)-1)]
print("miu:", list_mu)
list_lamda = [1-i for i in list_mu]
print("lambda:", list_lamda)
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 = []
M = []
copy_b = []
if boundary_type == 0: # 自然边界条件
a1 = 0
a2 = 0
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)
copy_b = b.copy()
print("g1~gn-1:", list_g)
M = ZGsolve(A,b)
M = [a1] + M + [a2]
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]))
copy_b = b.copy()
print("g0~gn:", copy_b)
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)
copy_b = b.copy()
print("g1~gn-1:", list_g)
M = ZGsolve(A,b)
M = [a1] + M + [a2]
print("A:", A)
print("b:", copy_b)
print("M:", M)
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__":
##############################################################################################################
list_x = [0,1,2,3] # 已给出的x数值，与y数值对应
list_y = [0,0,0,0] # 已给出的y数值，与x数值对应
# 记得改后三个参数，分别为边界条件类型(0=自然边界条件，1=一阶导数边界条件，2=二阶导数边界条件)，边界条件a1和a2
M,list_h = CubicSplineInterpolation(list_x,list_y,2,1,0)
print("二阶导数边界条件:")
PrintResult(M,list_h,list_x,list_y)
M,list_h = CubicSplineInterpolation(list_x,list_y,1,1,0)
print("一阶导数边界条件:")
PrintResult(M,list_h,list_x,list_y)

93
按方法整理/插值-牛顿前插-牛顿基本插值.py Normal file
View File

@@ -0,0 +1,93 @@
import math
# 求差分表
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,FxDiff_n):
h = list_x[1] - list_x[0]
t = (x - list_x[0]) / h
ks = max([abs(FxDiff_n(list_x[i],n)) for i in range(n)])
result = 1
for i in range(n):
result *= (t-i)*h/(i+1)
return result * ks
# 牛顿前插
def NewtonForwardInterpolation(x,n,list_x,list_y,FxDiff_n):
list_dyk = GetDyList(list_y)
print("\n差分表")
print("--------------------------------------------------")
for i in range(len(list_dyk)):
print(list_dyk[i])
print("--------------------------------------------------")
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,FxDiff_n))
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)
print("\n差商表")
print("--------------------------------------------------")
for i in range(len(list_dqy)):
print(list_dqy[i])
print("--------------------------------------------------")
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
# 定义原函数和其导函数计算结果，用于计算牛顿前插的截断误差，如果不需要则不用管
def FxDiff_n1(x,n):
result = 0
if n == 0:
# 下面改成原函数 ############################################################
result = math.exp(x)
else:
# 下面改成n阶导数 ##############################################################################
result = math.exp(x)
return result
if __name__ == '__main__':
##############################################################################################################
list_x = [0.0,0.2,0.4,0.6,0.8] # 已给出的x数值，与y数值对应
list_y = [1.0,1.2214,1.4918,1.8221,2.2255] # 已给出的y数值，与x数值对应
x_to_predict = 0.12 # 要预测的x值
# 记得改下面的点数(第二个参数)，n是点数，阶数为n-1 ################################
print("三点牛顿前插结果为%f, 截断误差%f" % NewtonForwardInterpolation(x_to_predict, 3, list_x, list_y,FxDiff_n1))
print("四点牛顿前插结果为%f, 截断误差%f" % NewtonForwardInterpolation(x_to_predict, 4, list_x, list_y,FxDiff_n1))
print("三点牛顿基本插值结果为%f" % NewtonBaseInterpolation(x_to_predict, 3, list_x, list_y))
print("四点牛顿基本插值结果为%f" % NewtonBaseInterpolation(x_to_predict, 4, list_x, list_y))

79
按方法整理/插值-线性-抛物线-拉格朗日多项式.py Normal file
View File

@@ -0,0 +1,79 @@
import math
# 获取与待求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,FxDiff_n):
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,FxDiff_n):
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,FxDiff_n):
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,FxDiff_n)
return (result,r)
# 抛物线插值
def ParabolaInterpolation(x,list_x,list_y,FxDiff_n):
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,FxDiff_n)
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
# 定义原函数和其导函数计算结果，用于计算插值的截断误差，如果不需要则不用管
def FxDiff_n1(x,n):
result = 0
if n == 0:
# 下面改成原函数 ############################################################
result = math.log(x)
else:
# 下面改成n阶导数 ##############################################################################
result = (-1)**(n+1) * math.factorial(n-1) / (x**n)
return result
if __name__ == "__main__":
##############################################################################
list_x = [10,11,12,13] # 已给出的x数值，与y数值对应
list_y = [2.3026,2.3979,2.4849,2.5649] # 已给出的y数值，与x数值对应
x_to_predict = 11.75 # 要预测的x值
print("线性插值结果为%f, 截断误差%f" % LinearInterpolation(x_to_predict,list_x,list_y,FxDiff_n1))
print("抛物线插值结果为%f, 截断误差%f" % ParabolaInterpolation(x_to_predict,list_x,list_y,FxDiff_n1))
# print(LagrangeInterpolation(x_to_predict,list_x,list_y))

38
按方法整理/数值积分-复合梯形-复合辛普生法.py Normal file
View File

@@ -0,0 +1,38 @@
# 复合Newton-Cotes公式 复合牛顿-柯特斯公式
# n等分参数，x1到x2的区间，type=1表示梯形法，type=2表示辛普生法
def CompositeNewtonCotes(x_start,x_end,fx,n, type):
if type == 1:
h = (x_end - x_start) / n
result = 0
for i in range(n):
result += (fx(x_start + i * h) + fx(x_start + (i + 1) * h))
result *= (h / 2)
return result
elif type == 2:
h = (x_end - x_start) / n
result = -fx(x_start) + fx(x_end)
for i in range(n):
result += (4 * fx(x_start + (i + 0.5) * h) + 2 * fx(x_start + i * h))
result *= (h / 6)
return result
# 积分原函数 ##############################################################
def fx(x):
if x == 0:
x = 1e-10 # Avoid division by zero #如果x能为0，注释掉这行##############
pass
return x/(4+x**2) #把函数改成题干的形式###################
if __name__ == "__main__":
##############################################################################################################
x_start = 3.0 # 积分下限
x_end = 6.0 # 积分上限
# 复合梯形公式，点数为n+1
print("复合梯形公式\n", CompositeNewtonCotes(x_start,x_end,fx,8, 1)) #8等分，1代表是梯形公式####################
# 复合辛普生公式，点数为2n+1
print("复合辛普生公式\n", CompositeNewtonCotes(x_start,x_end,fx,4, 2)) #4等分，2代表是辛普生公式###############

38
按方法整理/数值积分-逐次分半梯形递推公式.py Normal file
View File

@@ -0,0 +1,38 @@
import math
# 逐次分半梯形递推公式
def SplitTrapezoidal(a,b,fx,err):
count = 1
t1 = (b-a)*(fx(a)+fx(b))/2
print(f"t{count}={t1}")
k = 1
while True:
tmp = 0
for i in range(1, 2**(k-1)+1):
tmp += fx(a + (b-a)*(2*i-1)/(2**k))
t2 = t1/2+(b-a)*tmp/(2**k)
count *= 2
print(f"t{count}={t2}")
if abs(t2-t1) < err:
break
t1 = t2
k += 1
return t2,k
# 积分原函数 ##############################################################
def fx(x):
if x == 0:
x = 1e-10 # Avoid division by zero #如果x能为0，注释掉这行##############
pass
return 1/x #把函数改成题干的形式###################
if __name__ == "__main__":
##############################################################################################################
x_start = 1 # 积分下限
x_end = 3 # 积分上限
err = 1e-2 # 精度要求 P106
result,k = SplitTrapezoidal(x_start, x_end, fx, err)
print(f"Result: {result},k={k}")

71
按方法整理/数值积分-龙贝格算法.py Normal file
View File

@@ -0,0 +1,71 @@
import math
# 龙贝格方法 积分
def Romberg(a, b, fx, err):
table = [[],[],[],[]]
t00 = (b-a)*(fx(a)+fx(b))/2
table[0].append(t00)
print(f"t0(0)={t00}")
t01 = t10 = t11 = t20 = t21 = t30 = t31 = 0
k = 1
while True:
tmp = 0
for i in range(1, 2**(k-1)+1):
tmp += fx(a + (b-a)*(2*i-1)/(2**k))
t01 = t00/2+(b-a)*tmp/(2**k)
print(f"t0({k})={t01}")
table[0].append(t01)
if k>1:
t11 = (4*t01-t00)/3
print(f"t1({k-1})={t11}")
table[1].append(t11)
if k>2:
t21 = (16*t11-t10)/15
print(f"t2({k-2})={t21}")
table[2].append(t21)
if k>3:
t31 = (64*t21-t20)/63
print(f"t3({k-3})={t31}")
table[3].append(t31)
if abs(t31-t30) < err:
break
t30 = t31
else:
t30 = (64*t21-t20)/63
print(f"t3(0)={t30}")
table[3].append(t30)
t20 = t21
else:
t20 = (16*t11-t10)/15
print(f"t2(0)={t20}")
table[2].append(t20)
t10 = t11
else:
t10 = (4*t01-t00)/3
print(f"t1(0)={t10}")
table[1].append(t10)
t00 = t01
k += 1
print("Romberg table:")
for i in range(len(table)):
print(f"t{i}: {table[i]}")
return t31, k
# 积分原函数 ##############################################################
def fx(x):
if x == 0:
x = 1e-10 # Avoid division by zero #如果x能为0，注释掉这行##############
pass
return math.sin(x)/x #把函数改成题干的形式###################
if __name__ == "__main__":
##############################################################################################################
x_start = 0 # 积分下限
x_end = 1 # 积分上限
err = 0.5e-6 # 精度要求 P113
result, k = Romberg(x_start, x_end, fx, err)
print(f"Result: {result}, k={k}")

89
按方法整理/矩阵-列主元高斯消元法.py Normal file
View File

@@ -0,0 +1,89 @@
import math
# 列主元高斯消元法
def SovleRowMain(A,b):
ks = 0.00000001
n = len(A)
if len(A[0]) != n:
print("A要为方阵")
return None,None,None,None
if len(b) != n:
print("b与A的行数不匹配")
return None,None,None,None
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:
print("A矩阵奇异，无法进行高斯消元")
return None,None,None,None
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:
print("A矩阵奇异，无法进行高斯消元")
return None,None,None,None
# 回代求解
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和b改成题干要求的#####################################
A = [
[0, 3, 4],
[1, -1, 1],
[2, 1, 2]
]
b = [1, 2, 3]
P,L,U,x = SovleRowMain(A, b)
print("P:")
prettyPrintMatrix(P)
print("L:")
prettyPrintMatrix(L)
print("U:")
prettyPrintMatrix(U)
print("x:")
print(x)

29
按方法整理/非线性方程-二分法.py Normal file
View File

@@ -0,0 +1,29 @@
import math
# 二分法求解方程的根
def SolveByDivTwo(x1,x2,fx,err):
count = 1
while abs(x2-x1) >= err:
x = (x1+x2)/2
print(f"k={count},a{count}={x1},b{count}={x2},x{count}={x},fx(x)={fx(x)}")
if fx(x) * fx(x1) < 0:
x2 = x
else:
x1 = x
count += 1
return (x1+x2)/2
# 原函数换成题干的 #############################
def fx(x):
return x**4-3*x+1
if __name__ == "__main__":
##############################################################################################################
x1 = 0.3 # x下限范围
x2 = 0.4 # x上限范围
err = 0.5e-5 # 精度要求 P125
root = SolveByDivTwo(x1, x2,fx,err)
print(f"Root: {root}")

65
按方法整理/非线性方程-抛物线法.py Normal file
View File

@@ -0,0 +1,65 @@
import math
#抛物线法解方程
def MullerSolve(fx,x0,x1,x2,err1,err2,N):
count = 0
f0 = fx(x0)
f1 = fx(x1)
f2 = fx(x2)
q = (x2 - x1) / (x1 - x0)
p = 0
a = 0
b = 0
c = 0
while True:
p = (x2 - x0) / (x1 - x0)
a = q**2 * f0 - q*p*f1 + q*f2
b = q**2 *f0 - p**2 *f1 + (p + q)*f2
c = p*f2
h1 = 0
if b < 0:
h1 = -2 * c / (b - (b**2 - 4*a*c)**0.5)
else:
h1 = -2 * c / (b + (b**2 - 4*a*c)**0.5)
x3 = x2 + h1 * (x2 - x1)
f3 = fx(x3)
print(f"k={count}: x{count}={x0:.5}, f(x{count})={f0:.5}; x{count+1}={x1:.5}, f(x{count+1})={f1:.5}; x{count+2}={x2:.5}, f(x{count+2})={f2:.5}; x{count+3}={x3:.5}, f(x{count+3})={f3:.5}")
print(f"p={p:.5}, q={q:.5}, a={a:.5}, b={b:.5}, c={c:.5}, h={h1:.5}")
k = err1 + 1
if abs(f3) < 1:
k = abs(x3 - x2)
else:
k = abs(x3 - x2) / abs(f3)
if abs(f3) < err2 or k < err1:
return x3, 1
count += 1
if count > N:
return None, 0
x0 = x1
x1 = x2
x2 = x3
f0 = f1
f1 = f2
f2 = f3
q = h1
if __name__ == "__main__":
err1 = 1e-5 # 精度要求 P152
err2 = 1e-5 # 精度要求 P152
N = 100 # 最大迭代次数
x0 = 0.3 # 初始值1
x1 = 0.5 # 初始值2
x2 = 0.4 # 初始值3
fx = lambda x: 8*x**4 - 8*x**2 + 1 #原函数
result, status = MullerSolve(fx, x0, x1, x2, err1, err2, N)
if status == 1:
print(f"fx收敛 解为: {result}")
else:
print("fx不收敛")

38
按方法整理/非线性方程-正割法.py Normal file
View File

@@ -0,0 +1,38 @@
import math
#正割法计算方程
def SecantSolve(fx, x0, x1, err=1e-10, N0=100):
count = 0
print(f"k={count}: x{count}={x0}, f(x{count})={fx(x0)}")
count += 1
print(f"k={count}: x{count}={x1}, f(x{count})={fx(x1)}")
count += 1
while abs(x1 - x0) > err or abs(fx(x1)) > err:
if fx(x1) == fx(x0):
return None,0
x2 = x1 - fx(x1) * (x1 - x0) / (fx(x1) - fx(x0))
print(f"k={count}: x{count}={x2}, f(x{count})={fx(x2)}")
count += 1
if count > N0:
return None,-1
x0 = x1
x1 = x2
return x2,1
if __name__ == "__main__":
##############################################################################################################
err = 1e-5 # 根的误差限 P148
N0 = 100 # 最大迭代次数
x0 = 0.3 # 初始值1
x1 = 0.4 # 初始值2
fx = lambda x: x**4 - 3*x + 1 #原函数
result, status = SecantSolve(fx, x0, x1, err, N0)
if status == 1:
print(f"fx收敛 解为: {result}")
elif status == -1:
print("fx不收敛")
else:
print("分母为0，无法收敛")

50
按方法整理/非线性方程-牛顿下山法.py Normal file
View File

@@ -0,0 +1,50 @@
import math
def NewtonDownHillSolve(fx, dfx, x0, err1,err2, N0,min_t):
count = 0
print(f"k={count}, x0={x0}")
x1 = x0 + 1 + err1
while abs(x1 - x0) > err1 or abs(fx(x1)) > err2:
t = 1
if abs(dfx(x1)) < 1e-10:
print("导数为0，无法下山")
return None, 0
print(f"当前点: x0={x0}")
while t >= min_t:
x1 = x0 - t * fx(x0) / dfx(x0)
print(f"下山: t={t}, x1={x1}, abs(fx(x1))={abs(fx(x1))}, abs(fx(x0))={abs(fx(x0))}")
if abs(fx(x1)) < abs(fx(x0)):
break
t *= 0.5
if t < min_t:
print("达到最小t，下山失败")
return None, -2
# x1 = x0 - fx(x0) / dfx(x0)
count += 1
print(f"k={count}, x{count}={x1},x1-x0={abs(x1-x0)}")
if count > N0:
return None, -1
x0 = x1
print(f"收敛: x1={x1}, fx(x1)={fx(x1)}")
return x1, 1
if __name__ == "__main__":
##############################################################################################################
err1 = 1e-5 # 根的误差限 见P147
err2 = 1e-5 # 残量精度 见P147
N0 = 100 # 最大迭代次数
min_t = 1e-10 # 最小t值
x0 = 0.6 # 初始值
fx = lambda x: x**3 - x - 1 # 原函数
dfx = lambda x: 3*x**2 - 1 # 导函数
result, status = NewtonDownHillSolve(fx, dfx, x0, err1, err2, N0, min_t)
if status == 1:
print(f"收敛 解为: {result}")
elif status == -1:
print("不收敛")
elif status == -2:
print("下山失败")
else:
print("导数为0，无法收敛")

50
按方法整理/非线性方程-牛顿法.py Normal file
View File

@@ -0,0 +1,50 @@
import math
# 牛顿方法求解方程
def NewtonSolve(fx, dfx, x0, err, N0):
count = 0
print(f"k={count}, x0={x0}")
x1 = x0 + 1 + err
while abs(x1 - x0) > err or abs(fx(x1)) > err: # 添加条件修正根误差太大的问题
if abs(dfx(x1)) < 1e-10:
return None, 0
x1 = x0 - fx(x0) / dfx(x0)
count += 1
print(f"k={count}, x{count}={x1},x1-x0={abs(x1-x0)}")
if count > N0:
return None, -1
x0 = x1
return x1, 1
# 查找根区间
def FindRootZone(fx,start,stop,step):
x = start
while x < stop:
if fx(x) * fx(x+step) < 0:
return x
x += step
return None
# 定义原函数 ###############################
def fx(x):
return x**2 + 10*math.cos(x)
# 定义其导函数 ###############################
def dfx(x):
return 2*x - 10*math.sin(x)
if __name__ == "__main__":
##############################################################################################################
err = 1e-5 # 根的误差限 见P141
N0 = 100 # 最大迭代次数
x0 = 1.6 # 可以指定初始值
# x0 = FindRootZone(fx, 0, 2, 0.1) # 也可以用找根函数给定范围找根的缩小范围
result,status = NewtonSolve(fx, dfx, x0, err,N0)
if status == 1:
print(f"fx收敛 解为: {result}")
elif status == -1:
print("fx不收敛")
else:
print("fx导数为0，无法收敛")

40
按方法整理/非线性方程-画图找零点.py Normal file
View File

@@ -0,0 +1,40 @@
import math
import matplotlib.pyplot as plt
# 绘制函数图像 并标记可能的零点
def DrawGraph(a, b, stepper):
x = [a + (b-a)*i*stepper for i in range(int(1/stepper+1))]
y = [fx(i) for i in x]
plt.xlim(a-abs(b-a)*0.1, b+abs(b-a)*0.1)
plt.plot(x, y)
plt.axhline(0, color='black', lw=0.5, ls='-')
plt.axvline(a, color='red', lw=0.5, ls='--')
plt.axvline(b, color='red', lw=0.5, ls='--')
plt.text(a, y[0], f'({a:.5f},{y[0]:.5f})', fontsize=8, ha='left')
plt.text(b, y[-1], f'({b:.5f},{y[-1]:.5f})', fontsize=8, ha='right')
for i in range(len(x)-1):
if y[i] * y[i+1] < 0:
print(f"可能存在零点: ({x[i]:.5f},{y[i]:.5f})和({x[i+1]:.5f},{y[i+1]:.5f})之间")
plt.plot((x[i]+x[i+1])/2, (y[i]+y[i+1])/2, 'ro', markersize=3)
plt.title("Graph of f(x)")
plt.xlabel("x")
plt.ylabel("f(x)")
plt.show()
return x, y
# 把原函数换成题干的形式，自己注意定义域########################
def fx(x):
return math.exp(x)-math.sin(x)
if __name__ == "__main__":
##############################################################################################################
a = -2*math.pi # 左边界
b = math.pi # 右边界
step = 0.00001 # 步长
print(f"边界点: {a}, {b}")
print(f"步长: {step}")
x, y = DrawGraph(a, b, step)

46
按方法整理/非线性方程-迭代法.py Normal file
View File

@@ -0,0 +1,46 @@
def fd1(x):
return 1+1/x**2
def fd2(x):
return 1/(x-1)**0.5
# 迭代法求解方程
def Renew(x,fd,err):
count=0
i = 0
try:
while True:
xk = fd(x)
print(f"当前迭代值: {xk}, 上一次迭代值: {x}, 误差: {abs(xk-x)}")
if abs(i) < (xk-x):
count += 1
else:
count = 0
if count > 10:
print("不收敛")
break
if abs(xk-x) < err:
return xk
i = xk-x
x = xk
except Exception as e:
print(f"发生错误: {e}")
return None
return None
# 定义迭代公式 x = fd(x) ###############################
def fd(x):
return 1+1/x**2
if __name__ == "__main__":
##############################################################################################################
x0 = 1.5 # 初始值
err = 1e-5 # 精度要求 P136
result = Renew(x0, fd, err)
if result is not None:
print(f"收敛 解为: {result:.5f}")
else:
print("不收敛")