Блог

Мне нужно вычислить собственные значения матрицы 8×8 и построить каждое из собственных значений для символьной переменной, встречающейся в матрице. Для матрицы, которую я использую, я получаю 8 разных собственных значений, где каждое представляет функцию в «W», которая является моей символической переменной.

Используя Python, я попытался вычислить собственные значения с помощью Scipy и Sympy, которые вроде как сработали, но результаты хранятся странным образом (по крайней мере, для меня, как новичка, пока мало разбирающегося в программировании), и я не нашел способа извлечь только одно собственное значение, чтобы построить его.

 import numpy as np
import sympy as sp

W = sp.Symbol('W')

w0=1/780
wl=1/1064

# This is my 8x8-matrix
A= sp.Matrix([[w0 3*wl, 2*W, 0, 0, 0, np.sqrt(3)*W, 0, 0],
    [2*W, 4*wl, 0, 0, 0, 0, 0, 0],
    [0, 0, 2*wl w0, np.sqrt(3)*W, 0, 0, 0, np.sqrt(2)*W],
    [0, 0, np.sqrt(3)*W, 3*wl, 0, 0, 0, 0],
    [0, 0, 0, 0, wl w0, np.sqrt(2)*W, 0, 0],
    [np.sqrt(3)*W, 0, 0, 0, np.sqrt(2)*W, 2*wl, 0, 0],
    [0, 0, 0, 0, 0, 0, w0, W],
    [0, 0, np.sqrt(2)*W, 0, 0, 0, W, wl]])

# Calculating eigenvalues
eva = A.eigenvals()
evaRR = np.array(list(eva.keys()))
eva1p = evaRR[0]   # <- this is my try to refer to the first eigenvalue
  

В конце я надеюсь получить график над «W», где интересный диапазон равен [-0.002 0.002]. Для тех, кто интересуется, это касается атомной физики, а W относится к частоте раби, и я смотрю на так называемые одетые состояния.

Вы ничего не делаете неправильно — я думаю, вы просто пойманы, поскольку ваши собственные значения выглядят такими запутанными и сложными.

 import numpy as np
import sympy as sp
import matplotlib.pyplot as plt

W = sp.Symbol('W')

w0=1/780
wl=1/1064

# This is my 8x8-matrix
A= sp.Matrix([[w0 3*wl, 2*W, 0, 0, 0, np.sqrt(3)*W, 0, 0],
    [2*W, 4*wl, 0, 0, 0, 0, 0, 0],
    [0, 0, 2*wl w0, np.sqrt(3)*W, 0, 0, 0, np.sqrt(2)*W],
    [0, 0, np.sqrt(3)*W, 3*wl, 0, 0, 0, 0],
    [0, 0, 0, 0, wl w0, np.sqrt(2)*W, 0, 0],
    [np.sqrt(3)*W, 0, 0, 0, np.sqrt(2)*W, 2*wl, 0, 0],
    [0, 0, 0, 0, 0, 0, w0, W],
    [0, 0, np.sqrt(2)*W, 0, 0, 0, W, wl]])

# Calculating eigenvalues
eva = A.eigenvals()
evaRR = np.array(list(eva.keys()))
# The above is copied from your question
# We have to answer what exactly the eigenvalue is in this case
print(type(evaRR[0])) # >>> Piecewise
# Okay, so it's a piecewise function (link to documentation below).
# In the documentation we see that we can use the .subs method to evaluate
# the piecewise function by substituting a symbol for a value. For instance,

print(evaRR[0].subs(W, 0)) # Will substitute 0 for W
# This prints out something really nasty with tons of fractions.. 
# We can evaluate this mess with sympy's numerical evaluation method, N
print(sp.N(evaRR[0].subs(W, 0))) 
# >>> 0.00222190090611143 - 6.49672880062804e-34*I

# That's looking more like it! Notice the e-34 exponent on the imaginary part...
# I think it's safe to assume we can just trim that off.
# This is done by setting the chop keyword to True when using N:
print(sp.N(evaRR[0].subs(W, 0), chop=True)) # >>> 0.00222190090611143


# Now let's try to plot each of the eigenvalues over your specified range
fig, ax = plt.subplots(3, 3) # 3x3 grid of plots (for our 8 e.vals)
ax = ax.flatten() # This is so we can index the axes easier

plot_range = np.linspace(-0.002, 0.002, 10) # Range from -0.002 to 0.002 with 10 steps
for n in range(8):
    current_eigenval = evaRR[n]
    # There may be a way to vectorize this computation, but I'm not familiar enough with sympy.
    evaluated_array = np.zeros(np.size(plot_range))
    # This will be our Y-axis (or W-value). It is set to be the same shape as
    # plot_range and is initally filled with all zeros.

    for i in range(np.size(plot_range)):
        evaluated_array[i] = sp.N(current_eigenval.subs(W, plot_range[i]), 
                                  chop=True)
        # The above line is evaluating your eigenvalue at a specific point,
        # approximating it numerically, and then chopping off the imaginary.
    ax[n].plot(plot_range, evaluated_array, "c-")
    ax[n].set_title("Eigenvalue #{}".format(n))
    ax[n].grid()

plt.tight_layout()
plt.show()
  

И, как и было обещано, кусочная документация.