90 lines
2.2 KiB
Python
90 lines
2.2 KiB
Python
#抛物线法解方程
|
|
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}") |