diff --git a/69-5 - 副本.py b/69-5 - 副本.py deleted file mode 100644 index 66f221a..0000000 --- a/69-5 - 副本.py +++ /dev/null @@ -1,92 +0,0 @@ -x = [0,1,2,3] -y = [0,0,0,0] - -def ZGsolve(A,b): - n = len(b) - beta = [0]*n - for i in range(n): - if i == 0: - beta[i] = A[i][2] / A[i][1] - else: - beta[i] = A[i][2] / (A[i][1] - A[i][0]*beta[i-1]) - - for i in range(n): - if i == 0: - b[i] = b[i] / A[i][1] - else: - b[i] = (b[i] - A[i][0]*b[i-1]) / (A[i][1] - A[i][0]*beta[i-1]) - - for i in range(n-2,-1,-1): - b[i] = b[i] - beta[i]*b[i+1] - - return b - - -def GetDList(list_r): - result = [] - for i in range(1,len(list_r)): - result.append(list_r[i] - list_r[i-1]) - return result - -def GetDQList(list_x, list_y): - result = [] - for i in range(1,len(list_y)): - result.append((list_y[i] - list_y[i-1]) / (list_x[i] - list_x[i-1])) - return result - -def CubicSplineInterpolation(list_x,list_y,boundary_type,a1,a2): - list_h = GetDList(list_x) - list_dqxy = GetDQList(list_x, list_y) - list_mu = [list_h[i]/(list_h[i]+list_h[i+1]) for i in range(len(list_h)-1)] - list_lamda = [1-i for i in list_mu] - list_g = [6*(list_dqxy[i+1]-list_dqxy[i])/(list_h[i+1]+list_h[i]) for i in range(len(list_h)-1)] - - A = [] - b = [] - if boundary_type == 0: # 自然边界条件 - pass - elif boundary_type == 1: # 一阶导数边界条件 - A.append([0,2,1]) - b.append(6/list_h[0]*(list_dqxy[1]-a1)) - for i in range(len(list_g)): - A.append([list_mu[i],2,list_lamda[i]]) - b.append(list_g[i]) - A.append([1,2,0]) - b.append(6/list_h[-1]*(a2-list_dqxy[-1])) - M = ZGsolve(A,b) - - elif boundary_type == 2: # 二阶导数边界条件 - A.append([0,2,list_lamda[0]]) - b.append(list_g[0]-list_mu[0]*a1) - for i in range(1,len(list_g)-1): - A.append([list_mu[i],2,list_lamda[i]]) - b.append(list_g[i]) - A.append([list_mu[-1],2,0]) - b.append(list_g[-1]-list_lamda[-1]*a2) - M = ZGsolve(A,b) - M = [a1] + M + [a2] - - return M,list_h - - -if __name__ == "__main__": - M,list_h = CubicSplineInterpolation(x,y,2,1,0) - - print("list_h:",list_h) - print("M:",M) - for i in range(3): - k1 = (M[i+1]-M[i])/6/list_h[i] - k2 = (M[i]*x[i+1]-M[i+1]*x[i])/2/list_h[i] - k3 = (3*M[i+1]*x[i]**2-3*M[i]*x[i+1]**2-6*y[i]+M[i]*list_h[i]**2+6*y[i+1]-M[i+1]*list_h[i]**2)/6/list_h[i] - k4 = (M[i]*x[i+1]**3-M[i+1]*x[i]**3+6*y[i]*x[i+1]-M[i]*list_h[i]**2*x[i+1]-6*y[i+1]*x[i]+M[i+1]*list_h[i]**2*x[i])/6/list_h[i] - # print(k1,k2,k3,k4) - print("S(x)=%.10f*x^3+%.10f*x^2+%.10f*x+%.10f"%(k1,k2,k3,k4)) - - - - - - - - - diff --git a/89-1 - 副本.py b/89-1 - 副本.py deleted file mode 100644 index 21c228f..0000000 --- a/89-1 - 副本.py +++ /dev/null @@ -1,116 +0,0 @@ -x = [2,4,6,8] -y = [2,11,28,40] - - -# 列主元高斯消元法 -def SovleRowMain(A,b): - ks = 0.00000001 - n = len(A) - if len(A[0]) != n: - raise ValueError("A要为方阵") - if len(b) != n: - raise ValueError("b与A的行数不匹配") - p = list(range(n)) - for i in range(n): - row_max = abs(A[i][i]) - row_max_index = i - for j in range(i + 1, n): - if abs(A[j][i]) > row_max: - row_max = abs(A[j][i]) - row_max_index = j - A[i], A[row_max_index] = A[row_max_index], A[i] - b[i], b[row_max_index] = b[row_max_index], b[i] - p[i], p[row_max_index] = p[row_max_index], p[i] - - if abs(A[i][i]) < ks: - raise ValueError("A矩阵奇异，无法进行高斯消元") - for j in range(i + 1, n): - m = A[j][i] / A[i][i] - A[j][i] = m - for k in range(i + 1, n): - A[j][k] -= m * A[i][k] - b[j] -= m * b[i] - - if abs(A[n - 1][n - 1]) < ks: - raise ValueError("A矩阵奇异，无法进行高斯消元") - - # 回代求解 - b[n - 1] /= A[n - 1][n - 1] - for i in range(n - 2, -1, -1): - for j in range(i + 1, n): - b[i] -= A[i][j] * b[j] - b[i] /= A[i][i] - return b - - - -def LeastSquares(list_x,list_y,n): - m = len(list_x) - x_n = [] - b = [] - for i in range(2*n+1): - x_n.append(0) - b.append(0) - for j in range(m): - x_n[i]+=(list_x[j]**i) - b[i]+=(list_y[j]*list_x[j]**i) - b = b[:n+1] - A = [] - for i in range(n+1): - tmp = [] - for j in range(n+1): - tmp.append(x_n[i+j]) - A.append(tmp) - return SovleRowMain(A,b) - -def CalculateY(list_x, coeff): - re = [] - for i in range(len(list_x)): - re.append(0) - for j in range(len(coeff)): - re[i] += coeff[j]*list_x[i]**j - return re - -def MeanSquareErr(list_y,list_y_approx): - m = len(list_y) - err = 0 - for i in range(m): - err += (list_y[i] - list_y_approx[i])**2 - return err**0.5 - -def MaxErr(list_y,list_y_approx): - m = len(list_y) - err = 0 - for i in range(m): - if abs(list_y[i] - list_y_approx[i]) > err: - err = abs(list_y[i] - list_y_approx[i]) - return err - -def PrintEquation(coeff): - n = len(coeff) - str_ = str(coeff[0]) + "+" - for i in range(0,n-1): - str_ += str(coeff[i+1]) - str_ += "*x^"+str(i+1)+"+" - str_ = str_[0:len(str_)-1] - print(str_) - -if __name__ == "__main__": - print("一次拟合") - coeff = LeastSquares(x,y,1) - PrintEquation(coeff) - y_approx = CalculateY(x, coeff) - print("MeanSquareErr:") - print(MeanSquareErr(y,y_approx)) - print("MaxErr:") - print(MaxErr(y,y_approx)) - - print("二次拟合") - coeff = LeastSquares(x,y,2) - PrintEquation(coeff) - y_approx = CalculateY(x, coeff) - print("MeanSquareErr:") - print(MeanSquareErr(y,y_approx)) - print("MaxErr:") - print(MaxErr(y,y_approx)) - \ No newline at end of file diff --git a/89-2 - 副本.py b/89-2 - 副本.py deleted file mode 100644 index bbedc1a..0000000 --- a/89-2 - 副本.py +++ /dev/null @@ -1,76 +0,0 @@ -import math - -x = [1,2,4,8,16,32,64] -y = [4.22,4.02,3.85,3.59,3.44,3.02,2.59] - - -# 列主元高斯消元法 -def SovleRowMain(A,b): - ks = 0.00000001 - n = len(A) - if len(A[0]) != n: - raise ValueError("A要为方阵") - if len(b) != n: - raise ValueError("b与A的行数不匹配") - p = list(range(n)) - for i in range(n): - row_max = abs(A[i][i]) - row_max_index = i - for j in range(i + 1, n): - if abs(A[j][i]) > row_max: - row_max = abs(A[j][i]) - row_max_index = j - A[i], A[row_max_index] = A[row_max_index], A[i] - b[i], b[row_max_index] = b[row_max_index], b[i] - p[i], p[row_max_index] = p[row_max_index], p[i] - - if abs(A[i][i]) < ks: - raise ValueError("A矩阵奇异，无法进行高斯消元") - for j in range(i + 1, n): - m = A[j][i] / A[i][i] - A[j][i] = m - for k in range(i + 1, n): - A[j][k] -= m * A[i][k] - b[j] -= m * b[i] - - if abs(A[n - 1][n - 1]) < ks: - raise ValueError("A矩阵奇异，无法进行高斯消元") - - # 回代求解 - b[n - 1] /= A[n - 1][n - 1] - for i in range(n - 2, -1, -1): - for j in range(i + 1, n): - b[i] -= A[i][j] * b[j] - b[i] /= A[i][i] - return b - - - -def LeastSquares(list_x,list_y,n): - m = len(list_x) - x_n = [] - b = [] - for i in range(2*n+1): - x_n.append(0) - b.append(0) - for j in range(m): - x_n[i]+=(list_x[j]**i) - b[i]+=(list_y[j]*list_x[j]**i) - b = b[:n+1] - A = [] - for i in range(n+1): - tmp = [] - for j in range(n+1): - tmp.append(x_n[i+j]) - A.append(tmp) - return SovleRowMain(A,b) - -if __name__ == "__main__": - # 取对数 ln(W) = ln(C)+lamda*ln(t) - ln_W = [math.log(i) for i in y] - ln_t = [math.log(i) for i in x] - coeff = LeastSquares(ln_t,ln_W,1) - C = math.exp(coeff[0]) - lamda = coeff[1] - print("C:", C) - print("lamda:", lamda) \ No newline at end of file diff --git a/89-3 - 副本.py b/89-3 - 副本.py deleted file mode 100644 index 1f5ceb5..0000000 --- a/89-3 - 副本.py +++ /dev/null @@ -1,69 +0,0 @@ - -x = [19,25,31,38,44] -y = [19.0,32.3,49.0,73.3,97.8] - -# 列主元高斯消元法 -def SovleRowMain(A,b): - ks = 0.00000001 - n = len(A) - if len(A[0]) != n: - raise ValueError("A要为方阵") - if len(b) != n: - raise ValueError("b与A的行数不匹配") - p = list(range(n)) - for i in range(n): - row_max = abs(A[i][i]) - row_max_index = i - for j in range(i + 1, n): - if abs(A[j][i]) > row_max: - row_max = abs(A[j][i]) - row_max_index = j - A[i], A[row_max_index] = A[row_max_index], A[i] - b[i], b[row_max_index] = b[row_max_index], b[i] - p[i], p[row_max_index] = p[row_max_index], p[i] - - if abs(A[i][i]) < ks: - raise ValueError("A矩阵奇异，无法进行高斯消元") - for j in range(i + 1, n): - m = A[j][i] / A[i][i] - A[j][i] = m - for k in range(i + 1, n): - A[j][k] -= m * A[i][k] - b[j] -= m * b[i] - - if abs(A[n - 1][n - 1]) < ks: - raise ValueError("A矩阵奇异，无法进行高斯消元") - - # 回代求解 - b[n - 1] /= A[n - 1][n - 1] - for i in range(n - 2, -1, -1): - for j in range(i + 1, n): - b[i] -= A[i][j] * b[j] - b[i] /= A[i][i] - return b - - - -def LeastSquares(list_x,list_y,n): - m = len(list_x) - x_n = [] - b = [] - for i in range(2*n+1): - x_n.append(0) - b.append(0) - for j in range(m): - x_n[i]+=(list_x[j]**i) - b[i]+=(list_y[j]*list_x[j]**i) - b = b[:n+1] - A = [] - for i in range(n+1): - tmp = [] - for j in range(n+1): - tmp.append(x_n[i+j]) - A.append(tmp) - return SovleRowMain(A,b) - -if __name__ == "__main__": - x_square = [i**2 for i in x] - coeff = LeastSquares(x_square, y, 1) - print("拟合方程: y = %.6f + %.6f*x^2" % (coeff[0], coeff[1])) \ No newline at end of file