diff --git a/按方法整理/矩阵-平方根法.py b/按方法整理/矩阵-平方根法.py new file mode 100644 index 0000000..c0ce8ed --- /dev/null +++ b/按方法整理/矩阵-平方根法.py @@ -0,0 +1,49 @@ +import math + +# 获取下三角矩阵的索引 +# 列 行 +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)] + print(L) + 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 + +#把A,b换成题干的数值########################################### +if __name__ == "__main__": + ########################################################################## + # 储存下三角矩阵 a11, a21, a22, a31, a32, a33 ...,按照这个格式写 + A = [4,2,2,-2,-3,14] + + + b = [10,5,4] # b向量,常数项 + + print("x:") + print(SqrtSolve(A,b)) diff --git a/按方法整理/矩阵-模范数-点积.py b/按方法整理/矩阵-模范数-点积.py new file mode 100644 index 0000000..a8f4fb6 --- /dev/null +++ b/按方法整理/矩阵-模范数-点积.py @@ -0,0 +1,35 @@ +#模 范数 +def Norm(x,v): + if len(x[0]) == 1: + if v == 1: + return sum([abs(i[0]) for i in x]) + elif v == 2: + return (sum([i[0]**2 for i in x]))**0.5 + elif v == float("inf"): + return max([abs(i[0]) for i in x]) + else: + if v == 1: + return max([sum([abs(x[i][j]) for i in range(len(x))]) for j in range(len(x[0]))]) + elif v == float("inf"): + return max([sum([abs(i) for i in x[j]]) for j in range(len(x))]) + return None +# 计算矩阵的点积 +def Dot(A,B): + if len(A[0]) != len(B): + return None + return [[sum([A[i][j] * B[j][k] for j in range(len(A[0]))]) for k in range(len(B[0]))] for i in range(len(A))] + +#改成题干的矩阵########################## +if __name__ == "__main__": + A = [[1, 3], [-2, 4]] # 注意储存形式 + x = [[1], [-1]] # 注意储存形式,单独一列的相量,每个数字都要中括号 + + print("Norm of x with v=1:", Norm(x, 1)) + print("Norm of x with v=inf:", Norm(x, float("inf"))) + print("Norm of x with v=2:", Norm(x, 2)) + print("Dot product of A and x:", Dot(A, x)) + print("Norm of Ax with v=2:", Norm(Dot(A, x), 2)) + print("Norm of A with v=inf:", Norm(A, float("inf"))) + print("Norm of A with v=1:", Norm(A, 1)) + + \ No newline at end of file diff --git a/按方法整理/矩阵-追赶法.py b/按方法整理/矩阵-追赶法.py new file mode 100644 index 0000000..7d68815 --- /dev/null +++ b/按方法整理/矩阵-追赶法.py @@ -0,0 +1,42 @@ +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]) + print("beta:") + print(beta[:-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]) + print("y:") + print(b) + for i in range(n-2,-1,-1): + b[i] = b[i] - beta[i]*b[i+1] + + return b + + +if __name__ == "__main__": + ########################################################################## + #把A,b换成题干的数值########################################### + # 储存追赶法A矩阵 + A = [ + [0,4,-1], + [-1,4,-1], + [-1,4,-1], + [-1,4,-1], + [-1,4,0] + ] + # A第一行前面有个0,最后一行后面有个0 + + b = [100,200,200,200,100] + print("x:") + print(ZGsolve(A,b)) \ No newline at end of file diff --git a/按方法整理/非线性方程-抛物线法.py b/按方法整理/非线性方程-抛物线法.py index ae34b05..e5ac561 100644 --- a/按方法整理/非线性方程-抛物线法.py +++ b/按方法整理/非线性方程-抛物线法.py @@ -50,6 +50,7 @@ def MullerSolve(fx,x0,x1,x2,err1,err2,N): if __name__ == "__main__": + ############################################################################################################## err1 = 1e-5 # 精度要求 P152 err2 = 1e-5 # 精度要求 P152 N = 100 # 最大迭代次数