solve_ivp 永远卡在这个问题上

Question

我正在尝试用scipy.integrate.solve_ivp 求解一个微分方程组。系统依赖于一个真正的自变量t和因变量，c_n(t) -s 通常是复杂的。问题是，求解器总是卡住，无论系统的维度如何（由n_max 确定）。这是设置：

# Constants
from scipy.constants import hbar as h_
n_max = 2
t_max = 1

# The derivative function
def dcdt(t, c):
    return (-1.0j/h_)*((V_mat*np.exp(1.0j*w_mat*t)) @ c)

# Initial conditions
c_0 = np.zeros(n_max, dtype = complex)
c_0[0] = 1.0

#  Solving the deal
t = np.linspace(0, t_max, 10)
c = solve_ivp(dcdt, [0, t_max], c_0, t_eval = t)

就这样，永远不会停止运行。
以下是V_mat 和w_mat 的示例矩阵：

>>> V_mat
array([[1.0000000e-09, 1.8008153e-56],
       [1.8008153e-56, 1.0000000e-09]])

>>> w_mat
array([[      0.        , -156123.07053024],
       [ 156123.07053024,       0.        ]])

您会注意到，V_mat 和 w_mat 是维度为 n_max 的二维方阵。

问题是否与矩阵中的大/非常小的值有关？还是与复杂值有关？

Answer 1

正如我已经猜到的那样，问题与我试图解决的微分方程中的大值有关，特别是由于-1.0j/h_ in

def dcdt(t, c):
    return (-1.0j/h_)*((V_mat*np.exp(1.0j*w_mat*t)) @ c)

其中h_ = 1.054e-34 是约简普朗克常数。重新调整方程并删除 h_ 可以解决问题。