Решение проблемы с дифференциальным уравнением с помощьюsolve_ivp в зависимости от t_span ⇐ Python

[Anonymous]

Я использую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 или что-то в этом роде. .
Я, конечно, подумал, что это проблема с выборкой, поэтому безуспешно попытался изменить ее на более мелкую выборку (например, 10 000 точек). Кажется, это происходит, если первая точка времени раньше -11.
Я думаю, что это проблема с математикой или реализацией решателя. Я пробовал перейти на другие методы, но, похоже, они дают такие же бессмысленные (~0) результаты.
Я мало что знаю о теории этих решателей, поэтому был бы рад, если бы кто-нибудь указал мне на направление решения этой проблемы или дайте мне знать, если это неразрешимо. Спасибо

Подробнее здесь: https://stackoverflow.com/questions/793 ... -on-t-span