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.

的二维方阵

问题是否与矩阵中的 large/very 小值有关？还是与复数值有关？

Answer 1

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

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

其中 h_ = 1.054e-34 是约化普朗克常数。重新缩放等式并删除 h_ 可解决问题。

solve_ivp 永远被这个问题卡住了

solve_ivp gets stuck, forever, with this problem

python

numpy

scipy

numerical-integration