scipy.integrate.ode 放弃整合

Question

我在使用 scipy.integrate.ode 时遇到了一个奇怪的问题。这是一个最小的工作示例：

import sys
import time
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
from scipy.integrate import ode, complex_ode

def fun(t):
    return np.exp( -t**2 / 2. )

def ode_fun(t, y):
    a, b = y
    f = fun(t)
    c = np.conjugate(b)
    dt_a = -2j*f*c + 2j*f*b
    dt_b = 1j*f*a
    return [dt_a, dt_b]

t_range = np.linspace(-10., 10., 10000)

init_cond = [-1, 0]
trajectory = np.empty((len(t_range), len(init_cond)), dtype=np.complex128)

### setup ###
r = ode(ode_fun).set_integrator('zvode', method='adams', with_jacobian=False)
r.set_initial_value(init_cond, t_range[0])
dt = t_range[1] - t_range[0]

### integration ###
for i, t_i in enumerate(t_range):
    trajectory[i,:] = r.integrate(r.t+dt)

a_traj    = trajectory[:,0]
b_traj    = trajectory[:,1]
fun_traj  = fun(t_range)

### plot ###
plt.figure(figsize=(10,5))
plt.subplot(121, title='ODE solution')
plt.plot(t_range, np.real(a_traj))
plt.subplot(122, title='Input')
plt.plot(t_range, fun_traj)
plt.show()

此代码运行正常，输出图为（ODE 明确依赖于输入变量，右面板显示输入，左面板显示第一个变量的 ode 解）。

所以原则上我的代码是有效的。奇怪的是，如果我简单地替换积分范围

t_range = np.linspace(-10., 10., 10000)

来自

t_range = np.linspace(-20., 20., 10000)

我得到输出

不知何故，集成商只是放弃了集成并将我的解决方案保留为常量。 为什么会这样？我该如何解决？

我测试过的一些东西：这显然不是分辨率问题，集成步骤已经很小了。相反，似乎集成器在几步之后甚至不再费心调用 ode 函数。我已经通过在 ode_fun().

中包含一个打印语句来测试它

我目前的怀疑是集成商认为我的函数在前几个集成步骤中没有发生显着变化后是恒定的。我可能必须在某处设置一些容差水平吗？

感谢任何帮助！

Answer 1

“我目前的怀疑是积分器在前几个积分步骤中没有显着变化后决定我的函数是恒定的。”你的怀疑是正确的。 ODE 求解器通常有一个内部步长，它根据求解器计算的误差估计进行自适应调整。这些步长与请求输出的时间无关；请求时间的输出是使用在内部步骤计算的点处对解进行插值来计算的。

当您在 t = -20 处启动求解器时，显然输入变化如此缓慢，以至于求解器的内部步长变得足够大，以至于当求解器接近 t = 0 时，求解器会向右跳过在输入脉冲上。

您可以使用 set_integrator 方法的选项 max_step 来限制内部步长。如果我将 max_step 设置为 2.0（例如），

r = ode(ode_fun).set_integrator('zvode', method='adams', with_jacobian=False,
                                max_step=2.0)

我得到了你期望的输出。

scipy.integrate.ode 放弃整合

scipy.integrate.ode gives up on integration

python

scipy

ode

integrator