根据t_span使用solve_ivp求解微分方程问题

Question

我正在使用

solve_ivp

，但根据问题的设置，我得到了奇怪的结果，尽管我认为这更多是求解器的实现问题而不是编码问题，但我希望有人提供一个输入。 所以我的代码如下：

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

omega_ir = 1.0 * 2 * np.pi
gamma_ir = 0.5


def pulse(t):
    E0 = 1
    w = 0.35
    return -E0 * np.exp(-((t - 0) ** 2) / (w**2)) * np.sin((t - 0) / 0.33)


def equations(t, y):
    Q_ir, P_ir = y

    # Equations of motion
    dQ_ir = P_ir
    dP_ir = -(omega_ir**2) * Q_ir - gamma_ir * P_ir + pulse(t)

    return [dQ_ir, dP_ir]


initial_conditions = [0, 0]

# Time span for the simulation
t_span = (-1, 40)
t_eval = np.linspace(t_span[0], t_span[1], 1000)

solution = solve_ivp(
    equations,
    t_span,
    initial_conditions,
    t_eval=t_eval,
)

Q_ir = solution.y[0]
print(solution.message)


fig, ax = plt.subplots()
ax.plot(t_eval, Q_ir / max(Q_ir))
ax.plot(t_eval, pulse(t_eval) / max(pulse(t_eval)))
ax.set_xlabel("Time")
ax.set_ylabel("Normalised intensity Intesnity")

plt.show()

所以这是一个简单的钟摆/振荡器问题。如果我运行它，它可以正常工作，消息是：

The solver successfully reached the end of the integration interval.

太棒了。绘图如下（蓝色是振荡器的位置，橙色是初始激励，为了清晰起见，两者均已标准化）： 很不错。 但是如果我尝试改变

t_span=(-15, 40)

剧情如下： 请注意，它仍然是标准化的，蓝线实际上是 ~1e-50 或类似的东西。 我当然认为这是采样的问题，所以我尝试将其更改为更精细的采样（例如10000点）但没有成功。如果第一个时间点早于 -11，就会发生这种情况。

我认为这是数学或求解器实现的问题。我尝试更改为其他方法，但它似乎给出了类似的无意义（〜0）结果。 我对这些求解器的理论了解不多，所以如果有人能指出我解决这个问题的方向，或者让我知道这个问题是否无法解决，我会很高兴。谢谢

Answer 1

Solve_ivp 使用自适应时间步长。这是基于事情发生的速度。如果您从任何脉冲之前开始，那么似乎什么也没有发生，求解器会快速增加时间步长，然后……当脉冲发生时会错过它！

一个简单的解决方案是限制solve_ivp中的时间步长。添加可选参数 max_step=0.1 （或任何其他足够小的参数），你应该没问题。