python - scipy odeint 结果取决于输入时间数组

Question

我正在尝试使用 scipy odeint 函数求解以下一维 ODE(二阶):

z'' = 2*F*cos(vt-kz)*(z*cos(vt-k*z) - k*(1+z^2)*sin(vt-kz))/(1+z^2)^2

我的 python 脚本是

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

def func(var,t,F,k,v):
    z,Vz = var
    tmp = v*t-k*z
    cos_tmp = np.cos(tmp)
    denom = 1+z**2

    return [Vz,
            2*F*cos_tmp*(z*cos_tmp - k*denom*np.sin(tmp))/(denom*denom)]

我知道 odeint 使用自适应时间步长，它返回在输入 t 数组给定的时间点评估的解决方案。但是，我观察到当时间数组不够精细时解决方案会出现分歧，如果 odeint 自动选择内部时间步长，这对我来说没有意义。

k= 12184.
F = -0.000125597154176
v = 0.3
y0 = [-1.0,0.]

t1 = np.linspace(0,10,10000)
t2 = np.linspace(0,10,100)
t3 = np.linspace(0,10,10)

result1 = odeint(func,y0,t1,args=(F,k,v))
result2 = odeint(func,y0,t2,args=(F,k,v))
result3 = odeint(func,y0,t3,args=(F,k,v))

fig,ax = plt.subplots(1,2,figsize=(14,5))

ax[0].plot(t1,result1[:,0],label='t1')
ax[0].plot(t2,result2[:,0],'r.',label='t2')
ax[1].plot(t3,result3[:,0],'r.',label='t3')

ax[0].legend(frameon=False,loc='upper left',numpoints=1)
ax[1].legend(frameon=False,loc='upper left',numpoints=1)

plt.show()

解决方案的图片在这里:

Answer 1

当您运行代码时，您是否注意到打印的警告？它说

 lsoda--  at current t (=r1), mxstep (=i1) steps   
       taken on this call before reaching tout     
      in above message,  i1 =       500
      in above message,  r1 =  0.3622415764602D+00
[...]/scipy/integrate/odepack.py:218: ODEintWarning: Excess work done on this call (perhaps wrong Dfun type). Run with full_output = 1 to get quantitative information.
  warnings.warn(warning_msg, ODEintWarning)

有点含糊，但它的意思是odeint在达到请求的时间之前达到 500 个内部时间步长的限制。

您可以通过设置 mxstep 参数来增加该最大值。我刚刚使用

尝试了您的代码

result3 = odeint(func,y0,t3,args=(F,k,v), mxstep=2000)

它奏效了。