修改龙格库塔bug
This commit is contained in:
29
按方法整理/常微分方程-经典龙格-库塔格式.py
Normal file
29
按方法整理/常微分方程-经典龙格-库塔格式.py
Normal file
@@ -0,0 +1,29 @@
|
||||
import math
|
||||
|
||||
#经典R-K，龙格-库塔法
|
||||
def ClassicRK(x0,y0,h,xk,fxy):
|
||||
k1=k2=k3=k4=0
|
||||
result = [(x0,y0)]
|
||||
while x0<=xk:
|
||||
k1 = fxy(x0,y0)
|
||||
k2 = fxy(x0+h/2,y0+h*k1/2)
|
||||
k3 = fxy(x0+h/2,y0+h*k2/2)
|
||||
k4 = fxy(x0+h,y0+h*k3)
|
||||
y0 += h*(k1+2*k2+2*k3+k4)/6
|
||||
x0 += h
|
||||
result.append((x0,y0))
|
||||
return result
|
||||
|
||||
|
||||
if __name__=="__main__":
|
||||
##########################################################################
|
||||
x0 = 0 # x的初始值（左边界）换成题干里面的#################
|
||||
y0 = -1 # y的初始值换成题干里面的#################
|
||||
fxy = lambda x,y: x + y #f(x,y)换成题干里面的#################
|
||||
real_fx = lambda x: -x-1 #真实函数换成题干里面的#################
|
||||
h = 0.1 #步长换成题干里面的#################
|
||||
xk = 2 #x的右边界换成题干里面的#################
|
||||
result = ClassicRK(x0, y0, h, xk, fxy)
|
||||
print("x\ty\t\treal_y\t\t\t误差")
|
||||
for x, y in result:
|
||||
print(f"{x:.2f}\t{y:.10f}\t\t{real_fx(x):.10f}\t\t\t{abs(y - real_fx(x)):.5f}")
|
||||
62
按方法整理/矩阵-SOR逐次超松弛迭代法.py
Normal file
62
按方法整理/矩阵-SOR逐次超松弛迭代法.py
Normal file
@@ -0,0 +1,62 @@
|
||||
import math
|
||||
|
||||
#模 范数
|
||||
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
|
||||
|
||||
#SOR方法 逐次超松弛迭代
|
||||
def SOR(A,b,x,w,err,N):
|
||||
count = 0
|
||||
n = len(A)
|
||||
while True:
|
||||
count += 1
|
||||
for i in range(n):
|
||||
sum1 = sum(A[i][j] * x[j] for j in range(i))
|
||||
sum2 = sum(A[i][j] * x[j] for j in range(i, n))
|
||||
x[i] += w*(b[i] - sum1 - sum2) / A[i][i]
|
||||
r = [[b[i] - sum(A[i][j] * x[j] for j in range(len(A[0])))] for i in range(n)]
|
||||
err_now = Norm(r, float("inf"))
|
||||
x_t = [round(i,5) for i in x]
|
||||
print(f"第{count}次迭代, 误差 = {err_now:.5}, x = {x_t}")
|
||||
if err_now < err:
|
||||
return x, count, 1
|
||||
|
||||
if count > N:
|
||||
return None,count, 0
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
##########################################################################
|
||||
#把矩阵改成题干的矩阵，b改成题干结果，err精度要求修改##########################
|
||||
A = [
|
||||
[4,-1,0,-1,0,0],
|
||||
[-1,4,-1,0,-1,0],
|
||||
[0,-1,4,0,0,-1],
|
||||
[-1,0,0,4,-1,0],
|
||||
[0,-1,0,-1,4,-1],
|
||||
[0,0,-1,0,-1,4]
|
||||
]
|
||||
b = [2,3,2,2,1,2]
|
||||
|
||||
x = [0,0, 0, 0, 0, 0] # 初始解
|
||||
err = 1e-5 # 精度要求
|
||||
w = 1 # 松弛因子,题干要求 P201
|
||||
x1,k,sta = SOR(A, b, x, w, err, 100)
|
||||
print(f"w = {w}, 解为: {x1}, 迭代次数: {k}, 状态: {'收敛' if sta == 1 else '未收敛'}")
|
||||
|
||||
w = 1.1
|
||||
x = [0, 0, 0, 0, 0, 0]
|
||||
x2,k,sta = SOR(A, b, x, w, err, 100)
|
||||
print(f"w = {w}, 解为: {x2}, 迭代次数: {k}, 状态: {'收敛' if sta == 1 else '未收敛'}")
|
||||
94
按方法整理/矩阵-特征值-谱半径.py
Normal file
94
按方法整理/矩阵-特征值-谱半径.py
Normal file
@@ -0,0 +1,94 @@
|
||||
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.real < 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)
|
||||
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
|
||||
#计算矩阵的行列式
|
||||
def Det(A):
|
||||
if len(A) == 2:
|
||||
return A[0][0] * A[1][1] - A[0][1] * A[1][0]
|
||||
det = 0
|
||||
for c in range(len(A)):
|
||||
sub_matrix = [row[:c] + row[c+1:] for row in A[1:]]
|
||||
det += ((-1) ** c) * A[0][c] * Det(sub_matrix)
|
||||
return det
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
##########################################################################
|
||||
#把矩阵换成题干的矩阵，通常建议自己再验算一遍，这个能不能算不一定#########################
|
||||
A =[
|
||||
[1,0,1],
|
||||
[2,2,1],
|
||||
[-1,0,0]
|
||||
]
|
||||
|
||||
lam = []
|
||||
|
||||
count = 0
|
||||
k = -100
|
||||
fx = lambda x: Det([[A[i][j] - x * (1 if i == j else 0) for j in range(len(A))] for i in range(len(A))])
|
||||
k = 10
|
||||
while len(lam) < len(A):
|
||||
for i in range(-k,k):
|
||||
re,sta = MullerSolve(fx, i, i + 1, i + 2, 1e-10, 1e-10, 100)
|
||||
if sta == 1:
|
||||
a = round(re.real, 9)
|
||||
b = round(re.imag, 9)
|
||||
re_t = complex(a, b)
|
||||
if re_t not in lam:
|
||||
if re_t.imag != 0:
|
||||
lam.append(re_t)
|
||||
lam.append(re_t.conjugate())
|
||||
else:
|
||||
lam.append(a)
|
||||
if len(lam) == len(A):
|
||||
break
|
||||
k *= 10
|
||||
p = abs(lam[0])
|
||||
for i in range(len(lam)):
|
||||
print(f"λ{i+1} = {lam[i]}")
|
||||
if abs(lam[i]) > p:
|
||||
p = abs(lam[i])
|
||||
|
||||
print(f"谱半径 = {p:.3f}")
|
||||
102
按方法整理/矩阵-迭代改善法.py
Normal file
102
按方法整理/矩阵-迭代改善法.py
Normal file
@@ -0,0 +1,102 @@
|
||||
import math
|
||||
|
||||
#列主元高斯消元法
|
||||
def SovleRowMain(A,b,round_num=15):
|
||||
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 = round(A[j][i] / A[i][i],round_num)
|
||||
A[j][i] = m
|
||||
for k in range(i + 1, n):
|
||||
A[j][k] -= round(m * A[i][k],round_num)
|
||||
b[j] -= round(m * b[i],round_num)
|
||||
|
||||
if abs(A[n - 1][n - 1]) < ks:
|
||||
print("A矩阵奇异，无法进行高斯消元")
|
||||
return None, None, None, None
|
||||
|
||||
# 回代求解
|
||||
b[n - 1] = round(b[n - 1]/A[n - 1][n - 1],round_num)
|
||||
for i in range(n - 2, -1, -1):
|
||||
for j in range(i + 1, n):
|
||||
b[i] -= round(A[i][j] * b[j],round_num)
|
||||
b[i] /= round(A[i][i])
|
||||
b = [round(b[i], round_num) for i in range(n)]
|
||||
# 得到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 IterativeMethod(A, b, err, N, fake_round_num=15):
|
||||
b_c = [b[i] for i in range(len(b))]
|
||||
A_c = [[A[i][j] for j in range(len(A[0]))] for i in range(len(A))]
|
||||
P,L,U,x0 = SovleRowMain(A_c, b_c,fake_round_num)
|
||||
print(L)
|
||||
print(U)
|
||||
print(f"初始解为: {x0}")
|
||||
count = 0
|
||||
while count<N:
|
||||
r1 = [b[i] - sum([A[i][j] * x0[j] for j in range(len(A[0]))]) for i in range(len(A))]
|
||||
A_c = [[A[i][j] for j in range(len(A[0]))] for i in range(len(A))]
|
||||
d1 = SovleRowMain(A_c, r1,fake_round_num)[3]
|
||||
x0 = [x0[i] + d1[i] for i in range(len(x0))]
|
||||
print(f"第{count+1}次迭代, r{count+1} = {r1}, d{count+1} = {d1}, x{count+2} = {x0}")
|
||||
err_now = max(abs(r1[i]) for i in range(len(r1)))
|
||||
count += 1
|
||||
if err_now < err:
|
||||
break
|
||||
return x0,count
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
##########################################################################
|
||||
#把矩阵改成题干的矩阵，b改成题干结果，err精度要求修改##########################
|
||||
A = [
|
||||
[51,82],
|
||||
[151/3,81]
|
||||
]
|
||||
b = [235,232]
|
||||
|
||||
err = 1e-4 # 精度要求
|
||||
N = 1000 # 迭代次数上限
|
||||
fake_round_num = 4 # 模拟的四舍五入精度，根据题目情况或者凑过程修改
|
||||
|
||||
x = IterativeMethod(A, b, err, N, fake_round_num)[0]
|
||||
print(f"解为: {x}")
|
||||
|
||||
# 先求范数与逆矩阵（条件数）
|
||||
|
||||
69
按方法整理/矩阵-逆矩阵-条件数.py
Normal file
69
按方法整理/矩阵-逆矩阵-条件数.py
Normal file
@@ -0,0 +1,69 @@
|
||||
import math
|
||||
|
||||
#模 范数
|
||||
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 Det(A):
|
||||
if len(A) == 2:
|
||||
return A[0][0] * A[1][1] - A[0][1] * A[1][0]
|
||||
det = 0
|
||||
for c in range(len(A)):
|
||||
sub_matrix = [row[:c] + row[c+1:] for row in A[1:]]
|
||||
det += ((-1) ** c) * A[0][c] * Det(sub_matrix)
|
||||
return det
|
||||
# 计算矩阵的逆矩阵
|
||||
def Inverse(A):
|
||||
n = len(A)
|
||||
# 计算代数余子式矩阵
|
||||
B = [[0 for i in range(n)] for j in range(n)]
|
||||
for i in range(n):
|
||||
for j in range(n):
|
||||
minor = [row[:j] + row[j+1:] for row in (A[:i] + A[i+1:])]
|
||||
B[j][i] = ((-1) ** (i + j)) * sum(minor[k][l] * (-1) ** (k + l) for k in range(n - 1) for l in range(n - 1))
|
||||
det = Det(A)
|
||||
print("det(A):",det)
|
||||
if det == 0:
|
||||
print("矩阵不可逆")
|
||||
return None
|
||||
A_inv = [[B[i][j] / det for j in range(n)] for i in range(n)]
|
||||
return A_inv
|
||||
|
||||
# 计算矩阵的条件数
|
||||
def Cond(A,v):
|
||||
inv_A = Inverse(A)
|
||||
print(f"inv_A: {inv_A}, Norm(A, v): {Norm(A, v)}, Norm(inv_A, v): {Norm(inv_A, v)}")
|
||||
# print(inv_A,Norm(A, v), Norm(inv_A, v))
|
||||
return Norm(A, v) * Norm(inv_A, v)
|
||||
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
##########################################################################
|
||||
# 把矩阵换成题干的矩阵 #########################
|
||||
A = [
|
||||
[1,2],
|
||||
[1.001,2.001]
|
||||
]
|
||||
# 把范数的种类数换成题干的要求，inf是无穷范数 #########################
|
||||
print(f"矩阵A的条件数为: {Cond(A, float('inf')):.5f}") # 1 1范数；2 2范数；float('inf') 无穷范数
|
||||
# 把矩阵换成题干的矩阵 #########################
|
||||
A = [
|
||||
[1,2],
|
||||
[3,4]
|
||||
]
|
||||
# 把范数的种类数换成题干的要求，inf是无穷范数 ########################
|
||||
print(f"矩阵A的条件数为: {Cond(A, float('inf')):.5f}") # 1 1范数；2 2范数；float('inf') 无穷范数
|
||||
Reference in New Issue
Block a user