I am trying to rewrite one code from Mathematica to Python. In Mathematica I have:
Solve[(K - w*M) . a == 0]
where K and M are square matrices, a = {x1, x2, x3, x4, x5, x6, x7}; and w will give the values in this case of 7 frequencies of the normal modes of vibrations of a system. Then, for each value of w, I get the relative positions of each x_i.
EDIT: I am inserting my full code:
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))
I do get the same numerical values for w that I get in Mathematica, although for some numerical values in the matrices K and M, the calculation does not finish or after two minutes I do not get any response. For other numerical values, I do get the solution in less than 10 seconds.
Now, when I try to get the position of x_i for each numerical value of w with the following lines:
.
.
.
wn=omegan**2
D1=K-wn*M
solve(D1*xn,xn)
I get for all cases of wn {x1: 0.0, x2: 0.0, x3: 0.0, x4: 0.0, x5: 0.0, x6: 0.0, x7: 0.0}
But with Mathematica I do get for example x1 -> 0.541196 x2, x3 -> 1.30656 x2, x4 -> 1.41421 x2, x5 -> 1.30656 x2, x6 -> x2, x7 -> 0.541196 x2.
I guess I am using a method not advisable in Python.
Now that the question is updated with complete code let's first see what those matrices look like:
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⎦
It looks like you are trying to solve a generalised eigenvalue problem i.e. to find scalars w_i^2 and nonzero vectors x_i such that
K x_i = w_i * M * x_i
SymPy does not have a dedicated solver for generalized eigenvalue problems like this but SciPy's eig function can do it for matrices of floating point numbers:
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]]))
Alternatively since the matrix M is invertible we can trivially turn this into an ordinary eigenvalue problem:
(M**-1 * K) x_i = w_i * x_i
Then we can do it with 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]])])
That is just showing the first (evs[0]) eigenvector and eigenvalue and 1 there is the multiplicty of the eigenvalue.
The advantage of using SymPy rather than SciPy in this calculation would be if you wanted to compute the results exactly. To do that first turn your floats into rationals:
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 ⎦⎦⎠
I won't show the rest because it is long but for completeness these are the approximate and exact eigenvalues as computed by 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 ⎪
↪ ⎭
These are the eigenvalues as computed by 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])