import numpy as np
import sympy
from sympy import symbols
import random
sympy.init_printing(use_unicode=True)
kp = 1.96*10**9
kb = 2.7048*10**9
Mb = 3577
L = 1.3
E = 4.393079*10**8
# 定义符号变量
w, k, x = symbols('w k x')
L, E = symbols('L E')
def f(x):
w_val = x*2*np.pi/60
k_val = kp*(kb-Mb*w_val**2)/(kp+kb-Mb*w_val**2)
# 将矩阵元素转换为符号变量
T1 = sympy.Matrix([
[1+(2940*w**2-k)*L**3/(6*E),L,L**2/(2*E),L**3/(6*E)],
[(2940*w**2-k)*L**2/(2*E),1,L/E,L**2/(2*E)],
[(2940*w**2-k)*L,0,1,L],
[2940*w**2-k,0,0,1]
])
B2 = sympy.Matrix([
[1+(5880*w**2)*L**3/(6*E),L,L**2/(2*E),L**3/(6*E)],
[(5880*w**2)*L**2/(2*E),1,L/E,L**2/(2*E)],
[(5880*w**2)*L,0,1,L],
[5880*w**2,0,0,1]
])
B23 = B2*B2
B4 = sympy.Matrix([
[1+(5880*w**2-k)*L**3/(6*E),L,L**2/(2*E),L**3/(6*E)],
[(5880*w**2-k)*L**2/(2*E),1,L/E,L**2/(2*E)],
[(5880*w**2-k)*L,0,1,L],
[5880*w**2-k,0,0,1]
])
T24 = B4*B23
T27 = T24*T24
T20 = T27*T24
T12 = B23
T3 = sympy.Matrix([
[1,0,0,0],
[0,1,0,0],
[0,0,1,0],
[2940*w**2-k,0,0,1]
])
T13 = T3*T12
T10 = T20*T1
T = T13*T10
# 计算行列式并化简
Tm = sympy.Matrix([
[T[2,0],T[2,1]],
[T[3,0],T[3,1]]
])
g = Tm.det().simplify()
# 将符号表达式转换为函数并求解
g_func = sympy.lambdify((w, k, L, E, x), g)
return g_func(w_val, k_val, L, E, x)
# 使用牛顿迭代法求解方程
x0 = 1.0 # 初始解
eps = 1e-10 # 精度
max_iter = 100 # 最大迭代次数
for i in range(max_iter):
fx = f(x0)
if abs(fx) < eps:
break
dfx = sympy.diff(g, x).subs(x, x0)
if dfx == 0:
break
x0 = x0 - fx/dfx
print('方程的解为:', x0)