From 6baaac12c577e92ea5ac90a9949f726eac2eb1cb Mon Sep 17 00:00:00 2001 From: lhye200 Date: Mon, 16 Jun 2025 20:44:29 +0800 Subject: [PATCH] aaa --- 按方法整理/常微分方程-阿达姆斯方法.py | 54 +++++++ .../常微分方程-阿达姆斯预测-校正方法.py | 65 +++++++++ 按方法整理/拟合-最小二乘法.py | 133 ++++++++++++++++++ 按方法整理/插值-三次样条.py | 122 ++++++++++++++++ 按方法整理/插值-牛顿前插-牛顿基本插值.py | 93 ++++++++++++ 按方法整理/插值-线性-抛物线-拉格朗日多项式.py | 79 +++++++++++ 按方法整理/数值积分-复合梯形-复合辛普生法.py | 38 +++++ 按方法整理/数值积分-逐次分半梯形递推公式.py | 38 +++++ 按方法整理/数值积分-龙贝格算法.py | 71 ++++++++++ 按方法整理/矩阵-列主元高斯消元法.py | 89 ++++++++++++ 按方法整理/非线性方程-二分法.py | 29 ++++ 按方法整理/非线性方程-抛物线法.py | 65 +++++++++ 按方法整理/非线性方程-正割法.py | 38 +++++ 按方法整理/非线性方程-牛顿下山法.py | 50 +++++++ 按方法整理/非线性方程-牛顿法.py | 50 +++++++ 按方法整理/非线性方程-画图找零点.py | 40 ++++++ 按方法整理/非线性方程-迭代法.py | 46 ++++++ 17 files changed, 1100 insertions(+) create mode 100644 按方法整理/常微分方程-阿达姆斯方法.py create mode 100644 按方法整理/常微分方程-阿达姆斯预测-校正方法.py create mode 100644 按方法整理/拟合-最小二乘法.py create mode 100644 按方法整理/插值-三次样条.py create mode 100644 按方法整理/插值-牛顿前插-牛顿基本插值.py create mode 100644 按方法整理/插值-线性-抛物线-拉格朗日多项式.py create mode 100644 按方法整理/数值积分-复合梯形-复合辛普生法.py create mode 100644 按方法整理/数值积分-逐次分半梯形递推公式.py create mode 100644 按方法整理/数值积分-龙贝格算法.py create mode 100644 按方法整理/矩阵-列主元高斯消元法.py create mode 100644 按方法整理/非线性方程-二分法.py create mode 100644 按方法整理/非线性方程-抛物线法.py create mode 100644 按方法整理/非线性方程-正割法.py create mode 100644 按方法整理/非线性方程-牛顿下山法.py create mode 100644 按方法整理/非线性方程-牛顿法.py create mode 100644 按方法整理/非线性方程-画图找零点.py create mode 100644 按方法整理/非线性方程-迭代法.py diff --git a/按方法整理/常微分方程-阿达姆斯方法.py b/按方法整理/常微分方程-阿达姆斯方法.py new file mode 100644 index 0000000..e4ecc0a --- /dev/null +++ b/按方法整理/常微分方程-阿达姆斯方法.py @@ -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}") \ No newline at end of file diff --git a/按方法整理/常微分方程-阿达姆斯预测-校正方法.py b/按方法整理/常微分方程-阿达姆斯预测-校正方法.py new file mode 100644 index 0000000..d3dc496 --- /dev/null +++ b/按方法整理/常微分方程-阿达姆斯预测-校正方法.py @@ -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}") \ No newline at end of file diff --git a/按方法整理/拟合-最小二乘法.py b/按方法整理/拟合-最小二乘法.py new file mode 100644 index 0000000..5d70403 --- /dev/null +++ b/按方法整理/拟合-最小二乘法.py @@ -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)) + \ No newline at end of file diff --git a/按方法整理/插值-三次样条.py b/按方法整理/插值-三次样条.py new file mode 100644 index 0000000..77a546e --- /dev/null +++ b/按方法整理/插值-三次样条.py @@ -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) + + + + diff --git a/按方法整理/插值-牛顿前插-牛顿基本插值.py b/按方法整理/插值-牛顿前插-牛顿基本插值.py new file mode 100644 index 0000000..66d4cc1 --- /dev/null +++ b/按方法整理/插值-牛顿前插-牛顿基本插值.py @@ -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)) \ No newline at end of file diff --git a/按方法整理/插值-线性-抛物线-拉格朗日多项式.py b/按方法整理/插值-线性-抛物线-拉格朗日多项式.py new file mode 100644 index 0000000..92782cc --- /dev/null +++ b/按方法整理/插值-线性-抛物线-拉格朗日多项式.py @@ -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)) diff --git a/按方法整理/数值积分-复合梯形-复合辛普生法.py b/按方法整理/数值积分-复合梯形-复合辛普生法.py new file mode 100644 index 0000000..e0bc6ec --- /dev/null +++ b/按方法整理/数值积分-复合梯形-复合辛普生法.py @@ -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代表是辛普生公式############### + + diff --git a/按方法整理/数值积分-逐次分半梯形递推公式.py b/按方法整理/数值积分-逐次分半梯形递推公式.py new file mode 100644 index 0000000..a1454c8 --- /dev/null +++ b/按方法整理/数值积分-逐次分半梯形递推公式.py @@ -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}") diff --git a/按方法整理/数值积分-龙贝格算法.py b/按方法整理/数值积分-龙贝格算法.py new file mode 100644 index 0000000..987ca24 --- /dev/null +++ b/按方法整理/数值积分-龙贝格算法.py @@ -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}") diff --git a/按方法整理/矩阵-列主元高斯消元法.py b/按方法整理/矩阵-列主元高斯消元法.py new file mode 100644 index 0000000..80fc456 --- /dev/null +++ b/按方法整理/矩阵-列主元高斯消元法.py @@ -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) + + \ No newline at end of file diff --git a/按方法整理/非线性方程-二分法.py b/按方法整理/非线性方程-二分法.py new file mode 100644 index 0000000..9441a49 --- /dev/null +++ b/按方法整理/非线性方程-二分法.py @@ -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}") \ No newline at end of file diff --git a/按方法整理/非线性方程-抛物线法.py b/按方法整理/非线性方程-抛物线法.py new file mode 100644 index 0000000..ae34b05 --- /dev/null +++ b/按方法整理/非线性方程-抛物线法.py @@ -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不收敛") diff --git a/按方法整理/非线性方程-正割法.py b/按方法整理/非线性方程-正割法.py new file mode 100644 index 0000000..92eeb79 --- /dev/null +++ b/按方法整理/非线性方程-正割法.py @@ -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，无法收敛") \ No newline at end of file diff --git a/按方法整理/非线性方程-牛顿下山法.py b/按方法整理/非线性方程-牛顿下山法.py new file mode 100644 index 0000000..7170889 --- /dev/null +++ b/按方法整理/非线性方程-牛顿下山法.py @@ -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，无法收敛") \ No newline at end of file diff --git a/按方法整理/非线性方程-牛顿法.py b/按方法整理/非线性方程-牛顿法.py new file mode 100644 index 0000000..78d53c7 --- /dev/null +++ b/按方法整理/非线性方程-牛顿法.py @@ -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，无法收敛") diff --git a/按方法整理/非线性方程-画图找零点.py b/按方法整理/非线性方程-画图找零点.py new file mode 100644 index 0000000..ec0671d --- /dev/null +++ b/按方法整理/非线性方程-画图找零点.py @@ -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) \ No newline at end of file diff --git a/按方法整理/非线性方程-迭代法.py b/按方法整理/非线性方程-迭代法.py new file mode 100644 index 0000000..3f172d3 --- /dev/null +++ b/按方法整理/非线性方程-迭代法.py @@ -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("不收敛") +