From bded38c54f79de182e4c77a28ee21d04ba5010f3 Mon Sep 17 00:00:00 2001 From: lhye200 Date: Wed, 16 Apr 2025 21:40:31 +0800 Subject: [PATCH] all --- .gitignore | 1 + 227-3.py | 82 +++++++++++++++++++ 227-4.py | 45 +++++++++++ 227-7.py | 37 +++++++++ 68-1.py | 77 ++++++++++++++++++ 69-3.py | 71 +++++++++++++++++ 69-5 - 副本.py | 92 +++++++++++++++++++++ 69-5.py | 91 +++++++++++++++++++++ 89-1 - 副本.py | 116 +++++++++++++++++++++++++++ 89-1.py | 116 +++++++++++++++++++++++++++ 89-2 - 副本.py | 76 ++++++++++++++++++ 89-2.py | 76 ++++++++++++++++++ 89-3 - 副本.py | 69 ++++++++++++++++ 89-3.py | 69 ++++++++++++++++ test.py | 7 ++ 注释完整版/68-1-zhushi.py | 163 ++++++++++++++++++++++++++++++++++++++ 注释完整版/69-3-zhushi.py | 132 ++++++++++++++++++++++++++++++ 17 files changed, 1320 insertions(+) create mode 100644 .gitignore create mode 100644 227-3.py create mode 100644 227-4.py create mode 100644 227-7.py create mode 100644 68-1.py create mode 100644 69-3.py create mode 100644 69-5 - 副本.py create mode 100644 69-5.py create mode 100644 89-1 - 副本.py create mode 100644 89-1.py create mode 100644 89-2 - 副本.py create mode 100644 89-2.py create mode 100644 89-3 - 副本.py create mode 100644 89-3.py create mode 100644 test.py create mode 100644 注释完整版/68-1-zhushi.py create mode 100644 注释完整版/69-3-zhushi.py diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..e1f17b1 --- /dev/null +++ b/.gitignore @@ -0,0 +1 @@ +png/* \ No newline at end of file diff --git a/227-3.py b/227-3.py new file mode 100644 index 0000000..9f18459 --- /dev/null +++ b/227-3.py @@ -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) + \ No newline at end of file diff --git a/227-4.py b/227-4.py new file mode 100644 index 0000000..94390e1 --- /dev/null +++ b/227-4.py @@ -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)) diff --git a/227-7.py b/227-7.py new file mode 100644 index 0000000..a28d919 --- /dev/null +++ b/227-7.py @@ -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)) \ No newline at end of file diff --git a/68-1.py b/68-1.py new file mode 100644 index 0000000..0671d53 --- /dev/null +++ b/68-1.py @@ -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)) diff --git a/69-3.py b/69-3.py new file mode 100644 index 0000000..7e5b2ae --- /dev/null +++ b/69-3.py @@ -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)) diff --git a/69-5 - 副本.py b/69-5 - 副本.py new file mode 100644 index 0000000..66f221a --- /dev/null +++ b/69-5 - 副本.py @@ -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)) + + + + + + + + + diff --git a/69-5.py b/69-5.py new file mode 100644 index 0000000..e2589f1 --- /dev/null +++ b/69-5.py @@ -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) + + + + diff --git a/89-1 - 副本.py b/89-1 - 副本.py new file mode 100644 index 0000000..21c228f --- /dev/null +++ b/89-1 - 副本.py @@ -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)) + \ No newline at end of file diff --git a/89-1.py b/89-1.py new file mode 100644 index 0000000..21c228f --- /dev/null +++ b/89-1.py @@ -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)) + \ No newline at end of file diff --git a/89-2 - 副本.py b/89-2 - 副本.py new file mode 100644 index 0000000..bbedc1a --- /dev/null +++ b/89-2 - 副本.py @@ -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) \ No newline at end of file diff --git a/89-2.py b/89-2.py new file mode 100644 index 0000000..bbedc1a --- /dev/null +++ b/89-2.py @@ -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) \ No newline at end of file diff --git a/89-3 - 副本.py b/89-3 - 副本.py new file mode 100644 index 0000000..1f5ceb5 --- /dev/null +++ b/89-3 - 副本.py @@ -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])) \ No newline at end of file diff --git a/89-3.py b/89-3.py new file mode 100644 index 0000000..1f5ceb5 --- /dev/null +++ b/89-3.py @@ -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])) \ No newline at end of file diff --git a/test.py b/test.py new file mode 100644 index 0000000..6c2802f --- /dev/null +++ b/test.py @@ -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)) \ No newline at end of file diff --git a/注释完整版/68-1-zhushi.py b/注释完整版/68-1-zhushi.py new file mode 100644 index 0000000..e481201 --- /dev/null +++ b/注释完整版/68-1-zhushi.py @@ -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)) diff --git a/注释完整版/69-3-zhushi.py b/注释完整版/69-3-zhushi.py new file mode 100644 index 0000000..f536fee --- /dev/null +++ b/注释完整版/69-3-zhushi.py @@ -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))