scipy.integrate.ode 放弃整合 [英] scipy.integrate.ode gives up on integration

问题描述

我在使用

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

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

由

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

我得到输出

所以不知何故，积分器只是放弃了积分，并将我的解决方案保留为常数.为什么会这样?我该如何解决?

我测试过的一些事情:这显然不是解决问题，集成步骤已经很小了.相反，似乎集成器在几步之后甚至不再打扰调用 ode 函数.我已经通过在 ode_fun() 中包含一个打印语句来测试这一点.

我目前的怀疑是积分器认为我的函数在最初的几个积分步骤中没有显着变化后是恒定的.我可能需要在某处设置一些容忍度吗?

感谢任何帮助！

解决方案

我目前的怀疑是，积分器认为我的函数在最初的几个积分步骤中没有显着变化后是恒定的." 你的怀疑是正确的.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)

我得到了你期望的输出.

I am encountering a strange problem with scipy.integrate.ode. Here is a minimal working example:

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()

This code works correctly and the output figure is (the ODE is explicitly dependent on the input variable, right panel shows the input, left panel the ode solution for the first variable).

So in principle my code is working. What is strange is that if I simply replace the integration range

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

by

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

I get the output

So somehow the integrator just gave up on the integration and left my solution as a constant. Why does this happen? How can I fix it?

Some things I've tested: It is clearly not a resolution problem, the integration steps are really small already. Instead, it seems that the integrator does not even bother calling the ode function anymore after a few steps. I've tested that by including a print statement in the ode_fun().

My current suspicion is that the integrator decided that my function is constant after it did not change significantly during the first few integration steps. Do I maybe have to set some tolerance levels somewhere?

Any help appreciated!

解决方案

"My current suspicion is that the integrator decided that my function is constant after it did not change significantly during the first few integration steps." Your suspicion is correct. ODE solvers typically have an internal step size that is adaptively adjusted based on error estimates computed by the solver. These step sizes are independent of the times for which the output is requested; the output at the requested times is computed using interpolation of the solution at the points computed at the internal steps.

When you start your solver at t = -20, apparently the input changes so slowly that the solver's internal step size becomes large enough that by the time the solver gets near t = 0, the solver skips right over the input pulse.

You can limit the internal step size with the option max_step of the set_integrator method. If I set max_step to 2.0 (for example),

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

I get the output that you expect.

这篇关于scipy.integrate.ode 放弃整合的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！