代码之家 › 专栏 › 技术社区 › mmonti

求解微分方程的solve_ivp问题取决于t_span

ode scipy python

mmonti · 技术社区 · 6 月前

我在用 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,这种情况似乎就会发生。