我正在尝试将一段代码从 Mathematica 重写为 Python。在 Mathematica 中我有：
求解[(K - w*M) 。 a == 0]
其中 K 和 M 是方阵，a = {x1, x2, x3, x4, x5, x6, x7}；在这种情况下，w 将给出系统正常振动模式的 7 个频率的值。然后，对于 w 的每个值，我得到每个 x_i 的相对位置。
编辑：我正在插入完整的代码：
from sympy.solvers import solve
from sympy import *
import numpy as np
m=3.0
k=1.50
w,x1,x2,x3,x4,x5,x6,x7 = symbols("w,x1,x2,x3,x4,x5,x6,x7")
M = Matrix([[m,0,0,0,0,0,0], [0,m,0,0,0,0,0],[0,0,m,0,0,0,0],[0,0,0,m,0,0,0],[0,0,0,0,m,0,0],[0,0,0,0,0,m,0],[0,0,0,0,0,0,m]])
K = Matrix([[2*k,-k,0,0,0,0,0], [-k,2*k,-k,0,0,0,0],[0,-k,2*k,-k,0,0,0],[0,0,-k,2*k,-k,0,0],[0,0,0,-k,2*k,-k,0],[0,0,0,0,-k,2*k,-k],[0,0,0,0,0,-k,2*k]])
xn = Matrix([x1,x2,x3,x4,x5,x6,x7])
D1=K-w*M
print("omega squared")
A=solve(D1.det(), w)
print(A)
res =(np.array(A))**0.5
for i in range(7):
omegan=float(res[i])
wn=round(omegan**2,5)
D1=K-wn*M
print(D1*xn)
print("\n** Positions x_i ", i+1, "for omega= ", np.sqrt(wn), " son; ",solve(D1*xn,xn))
我确实得到了与在 Mathematica 中得到的 w 相同的数值，尽管对于矩阵 K 和 M 中的某些数值，计算未完成或两分钟后我没有得到任何响应。对于其他数值，我确实在不到 10 秒的时间内得到了答案。
现在，当我尝试使用以下几行获取 w 的每个数值的 x_i 的位置时：
.
.
.
wn=omegan**2
D1=K-wn*M
solve(D1*xn,xn)
我得到了所有 wn 的情况 {x1：0.0，x2：0.0，x3：0.0，x4：0.0，x5：0.0，x6：0.0，x7：0.0}
但是使用 Mathematica，我确实得到了 x1 -> 0.541196 x2、x3 -> 1.30656 x2、x4 -> 1.41421 x2、x5 -> 1.30656 x2、x6 -> x2、x7 -> 0.541196 x2。
我想我使用的方法在Python中是不可取的。
现在问题已更新为完整的代码，让我们首先看看这些矩阵是什么样的：
In [2]: K
Out[2]:
⎡3.0 -1.5 0 0 0 0 0 ⎤
⎢ ⎥
⎢-1.5 3.0 -1.5 0 0 0 0 ⎥
⎢ ⎥
⎢ 0 -1.5 3.0 -1.5 0 0 0 ⎥
⎢ ⎥
⎢ 0 0 -1.5 3.0 -1.5 0 0 ⎥
⎢ ⎥
⎢ 0 0 0 -1.5 3.0 -1.5 0 ⎥
⎢ ⎥
⎢ 0 0 0 0 -1.5 3.0 -1.5⎥
⎢ ⎥
⎣ 0 0 0 0 0 -1.5 3.0 ⎦
In [3]: M
Out[3]:
⎡3.0 0 0 0 0 0 0 ⎤
⎢ ⎥
⎢ 0 3.0 0 0 0 0 0 ⎥
⎢ ⎥
⎢ 0 0 3.0 0 0 0 0 ⎥
⎢ ⎥
⎢ 0 0 0 3.0 0 0 0 ⎥
⎢ ⎥
⎢ 0 0 0 0 3.0 0 0 ⎥
⎢ ⎥
⎢ 0 0 0 0 0 3.0 0 ⎥
⎢ ⎥
⎣ 0 0 0 0 0 0 3.0⎦
看起来您正在尝试解决广义特征值问题，即找到标量
w_i^2x_iK x_i = w_i * M * x_i
SymPy 没有专门的求解器来解决此类广义特征值问题，但 SciPy 的 eig 函数可以用于浮点数矩阵：
https://docs.scipy.org/doc/scipy/reference/ generated/scipy.linalg.eig.html
In [5]: from scipy.linalg import eig
In [11]: import numpy as np
In [12]: Kn = np.array(K, dtype=float)
In [13]: Mn = np.array(M, dtype=float)
In [14]: type(Kn), type(Mn)
Out[14]: (numpy.ndarray, numpy.ndarray)
In [15]: Kn
Out[15]:
array([[ 3. , -1.5, 0. , 0. , 0. , 0. , 0. ],
[-1.5, 3. , -1.5, 0. , 0. , 0. , 0. ],
[ 0. , -1.5, 3. , -1.5, 0. , 0. , 0. ],
[ 0. , 0. , -1.5, 3. , -1.5, 0. , 0. ],
[ 0. , 0. , 0. , -1.5, 3. , -1.5, 0. ],
[ 0. , 0. , 0. , 0. , -1.5, 3. , -1.5],
[ 0. , 0. , 0. , 0. , 0. , -1.5, 3. ]])
In [16]: Mn
Out[16]:
array([[3., 0., 0., 0., 0., 0., 0.],
[0., 3., 0., 0., 0., 0., 0.],
[0., 0., 3., 0., 0., 0., 0.],
[0., 0., 0., 3., 0., 0., 0.],
[0., 0., 0., 0., 3., 0., 0.],
[0., 0., 0., 0., 0., 3., 0.],
[0., 0., 0., 0., 0., 0., 3.]])
In [17]: eig(Kn, Mn)
Out[17]:
(array([0.07612047+0.j, 0.29289322+0.j, 0.61731657+0.j, 1. +0.j,
1.92387953+0.j, 1.70710678+0.j, 1.38268343+0.j]),
array([[ 1.91341716e-01, 3.53553391e-01, -4.61939766e-01,
5.00000000e-01, -1.91341716e-01, 3.53553391e-01,
4.61939766e-01],
[ 3.53553391e-01, 5.00000000e-01, -3.53553391e-01,
-1.93816282e-16, 3.53553391e-01, -5.00000000e-01,
-3.53553391e-01],
[ 4.61939766e-01, 3.53553391e-01, 1.91341716e-01,
-5.00000000e-01, -4.61939766e-01, 3.53553391e-01,
-1.91341716e-01],
[ 5.00000000e-01, -7.68527201e-16, 5.00000000e-01,
-1.91221591e-16, 5.00000000e-01, -4.23414829e-16,
5.00000000e-01],
[ 4.61939766e-01, -3.53553391e-01, 1.91341716e-01,
5.00000000e-01, -4.61939766e-01, -3.53553391e-01,
-1.91341716e-01],
[ 3.53553391e-01, -5.00000000e-01, -3.53553391e-01,
-3.60198282e-16, 3.53553391e-01, 5.00000000e-01,
-3.53553391e-01],
[ 1.91341716e-01, -3.53553391e-01, -4.61939766e-01,
-5.00000000e-01, -1.91341716e-01, -3.53553391e-01,
4.61939766e-01]]))
或者，由于矩阵 M 是可逆的，我们可以简单地将其转化为普通的特征值问题：
(M**-1 * K) x_i = w_i * x_i
然后我们可以用 SymPy 来做：
In [22]: evs = (M.inv() * K).eigenvects()
In [23]: print(evs[0])
(0.0761204674887132, 1, [Matrix([
[-0.191341716182545],
[-0.353553390593274],
[-0.461939766255643],
[ -0.5],
[-0.461939766255643],
[-0.353553390593274],
[-0.191341716182545]])])
这只是显示第一个（
evs[0]1在此计算中使用 SymPy 而不是 SciPy 的优点是如果您想精确计算结果。要做到这一点，首先将你的浮点数变成有理数：
In [24]: Mr = M.applyfunc(nsimplify)
In [25]: Kr = K.applyfunc(nsimplify)
In [26]: Mr
Out[26]:
⎡3 0 0 0 0 0 0⎤
⎢ ⎥
⎢0 3 0 0 0 0 0⎥
⎢ ⎥
⎢0 0 3 0 0 0 0⎥
⎢ ⎥
⎢0 0 0 3 0 0 0⎥
⎢ ⎥
⎢0 0 0 0 3 0 0⎥
⎢ ⎥
⎢0 0 0 0 0 3 0⎥
⎢ ⎥
⎣0 0 0 0 0 0 3⎦
In [27]: Kr
Out[27]:
⎡ 3 -3/2 0 0 0 0 0 ⎤
⎢ ⎥
⎢-3/2 3 -3/2 0 0 0 0 ⎥
⎢ ⎥
⎢ 0 -3/2 3 -3/2 0 0 0 ⎥
⎢ ⎥
⎢ 0 0 -3/2 3 -3/2 0 0 ⎥
⎢ ⎥
⎢ 0 0 0 -3/2 3 -3/2 0 ⎥
⎢ ⎥
⎢ 0 0 0 0 -3/2 3 -3/2⎥
⎢ ⎥
⎣ 0 0 0 0 0 -3/2 3 ⎦
In [28]: evs = (Mr.inv() * Kr).eigenvects()
In [29]: evs[0]
Out[29]:
⎛ ⎡⎡-1⎤⎤⎞
⎜ ⎢⎢ ⎥⎥⎟
⎜ ⎢⎢0 ⎥⎥⎟
⎜ ⎢⎢ ⎥⎥⎟
⎜ ⎢⎢1 ⎥⎥⎟
⎜ ⎢⎢ ⎥⎥⎟
⎜1, 1, ⎢⎢0 ⎥⎥⎟
⎜ ⎢⎢ ⎥⎥⎟
⎜ ⎢⎢-1⎥⎥⎟
⎜ ⎢⎢ ⎥⎥⎟
⎜ ⎢⎢0 ⎥⎥⎟
⎜ ⎢⎢ ⎥⎥⎟
⎝ ⎣⎣1 ⎦⎦⎠
In [30]: evs[1]
Out[30]:
⎛ ⎡⎡-1 ⎤⎤⎞
⎜ ⎢⎢ ⎥⎥⎟
⎜ ⎢⎢-√2⎥⎥⎟
⎜ ⎢⎢ ⎥⎥⎟
⎜ ⎢⎢-1 ⎥⎥⎟
⎜ √2 ⎢⎢ ⎥⎥⎟
⎜1 - ──, 1, ⎢⎢ 0 ⎥⎥⎟
⎜ 2 ⎢⎢ ⎥⎥⎟
⎜ ⎢⎢ 1 ⎥⎥⎟
⎜ ⎢⎢ ⎥⎥⎟
⎜ ⎢⎢√2 ⎥⎥⎟
⎜ ⎢⎢ ⎥⎥⎟
⎝ ⎣⎣ 1 ⎦⎦⎠
In [31]: evs[2]
Out[31]:
⎛ ⎡⎡ 1 ⎤⎤⎞
⎜ ⎢⎢ ⎥⎥⎟
⎜ ⎢⎢ ________ ⎥⎥⎟
⎜ ⎢⎢ ╲╱ 2 - √2 ⎥⎥⎟
⎜ ⎢⎢ ⎥⎥⎟
⎜ ⎢⎢ 2 ⎥⎥⎟
⎜ ⎢⎢ ⎛ ________⎞ ⎥⎥⎟
⎜ ⎢⎢ ⎜ ╲╱ 2 - √2 ⎟ ________ ⎥⎥⎟
⎜ ⎢⎢ -5 + 4⋅⎜1 - ──────────⎟ + 4⋅╲╱ 2 - √2 ⎥⎥⎟
⎜ ⎢⎢ ⎝ 2 ⎠ ⎥⎥⎟
⎜ ⎢⎢ ⎥⎥⎟
⎜ ________ ⎢⎢ 3 2⎥⎥⎟
⎜ ╲╱ 2 - √2 ⎢⎢ ⎛ ________⎞ ⎛ ________⎞ ⎥⎥⎟
⎜1 - ──────────, 1, ⎢⎢ ⎜ ╲╱ 2 - √2 ⎟ ________ ⎜ ╲╱ 2 - √2 ⎟ ⎥⎥⎟
⎜ 2 ⎢⎢-16 - 8⋅⎜1 - ──────────⎟ + 10⋅╲╱ 2 - √2 + 24⋅⎜1 - ──────────⎟ ⎥⎥⎟
⎜ ⎢⎢ ⎝ 2 ⎠ ⎝ 2 ⎠ ⎥⎥⎟
⎜ ⎢⎢ ⎥⎥⎟
⎜ ⎢⎢ 2 ⎥⎥⎟
⎜ ⎢⎢ ⎛ ________⎞ ⎥⎥⎟
⎜ ⎢⎢ ⎜ ╲╱ 2 - √2 ⎟ ________ ⎥⎥⎟
⎜ ⎢⎢ -5 + 4⋅⎜1 - ──────────⎟ + 4⋅╲╱ 2 - √2 ⎥⎥⎟
⎜ ⎢⎢ ⎝ 2 ⎠ ⎥⎥⎟
⎜ ⎢⎢ ⎥⎥⎟
⎜ ⎢⎢ ________ ⎥⎥⎟
⎜ ⎢⎢ ╲╱ 2 - √2 ⎥⎥⎟
⎜ ⎢⎢ ⎥⎥⎟
⎝ ⎣⎣ 1 ⎦⎦⎠
我不会显示其余部分，因为它很长，但为了完整性，这些是由 SymPy 计算的近似和精确特征值：
In [34]: (M.inv() * K).eigenvals()
Out[34]:
{0.0761204674887132: 1, 0.292893218813452: 1, 0.61731656763491: 1, 1.0: 1, 1.38268343236509: 1, 1.7 ↪
↪ 0710678118655: 1, 1.92387953251129: 1}
In [35]: (Mr.inv() * Kr).eigenvals()
Out[35]:
⎧ ________ ________ ________ ___ ↪
⎪ √2 ╱ 1 √2 ╱ √2 1 √2 ╱ 1 √2 ╱ √2 ↪
⎨1: 1, 1 - ──: 1, 1 - ╱ ─ - ── : 1, 1 - ╱ ── + ─ : 1, ── + 1: 1, ╱ ─ - ── + 1: 1, ╱ ── ↪
⎪ 2 ╲╱ 2 4 ╲╱ 4 2 2 ╲╱ 2 4 ╲╱ 4 ↪
⎩ ↪
↪ _____ ⎫
↪ 1 ⎪
↪ + ─ + 1: 1⎬
↪ 2 ⎪
↪ ⎭
这些是由 SciPy 计算的特征值：
In [36]: from scipy.linalg import eigvals
In [37]: eigvals(Kn, Mn)
Out[37]:
array([0.07612047+0.j, 0.29289322+0.j, 0.61731657+0.j, 1. +0.j,
1.92387953+0.j, 1.70710678+0.j, 1.38268343+0.j])