Python - 时间相关系数的微分方程求解器为不同的时间偏移提供不同的动态

如何解决Python - 时间相关系数的微分方程求解器为不同的时间偏移提供不同的动态

我正在求解系统与脉冲相互作用时的动力学问题，这基本上是在求解与时间相关的微分方程。一般来说，它工作正常，但每当我将脉冲的带宽取小，即大约为 1 时，求解器取决于脉冲从 t0 开始的位置。让我给你代码和一些图片

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import matplotlib as mpl
from scipy.integrate import solve_ivp
import scipy as scipy

## pulses 

def Box(t):
    return (abs(t-t0) < tau)
    #return  np.exp(-(t - t0)**2.0/(tau ** 2.0))


## Differential eq.
def Leq(t,u,pulse): 
    v=u[:16].reshape(4,4)
    a0=u[16] 
    da0=-a0+pulse(t)*E_0
    
    M=np.array([[-1,a0,0],\
                [0,-1,-a0],\
                [a0,-a0,-1]])*kappa
    Dr=np.array([[1,1j,\
                [1j,1,1j],1]])*kappa/2.0

    
    dv=M.dot(v)+v.dot(M)+Dr
    
    return np.concatenate([dv.flatten(),[da0]])
## Covariance matrix

cov0=np.array([[1,1]])/4 ##initial vector
cov0 = cov0.reshape(-1);   ## vectorize initial vector

a0_in=np.zeros((1,1),dtype=np.complex64) ##initial vector
a0_in = a0_in.reshape(-1);   ## vectorize initial vector
u_0=np.concatenate([cov0,a0_in])

### Parameters 
kappa=1.0 ##adimenstional kappa: kappa0/kappa 
tau=1.0  ##bandwidth pump

#Eth=kappa0*kappa/g  ##threshold intensity 
E_0=0.8 # pump intensity normalized to threshold 

t0=2.0   ##off set 

Tmax=10 ##max value for time
Nmax=100000 ##number of steps
dt=Tmax/Nmax  ##increment of time

t=np.linspace(0.0,Tmax,Nmax+1)

Box_sol=solve_ivp(Leq,[min(t),max(t)],u_0,t_eval= t,args=(Box,))
print(Box_sol.y[0:16,int(Nmax/2)].reshape(4,4))

然后只是一些期望值和绘图。

从图片中可以看出，参数是一样的，唯一不同的是时移t0。似乎只要脉冲的带宽很小，我就必须从非常接近 t=0 开始，我不完全理解为什么。不应该这样。

是求解器的问题吗？如果是这样，我能做些什么来解决它？我现在担心这个问题会出现在我正在运行的其他模拟中，而我没有意识到。

谢谢，琼。

解决方法

这是一个众所周知的行为，关于这个话题有几个问题。

简而言之，就是步长控制器。其背后的假设是 ODE 函数在高阶上是平滑的，并且局部行为在中等范围内通知全局行为。因此，如果您开始平坦，随着更高的导数消失，步长会迅速增加。如果情况不幸，这将在积分器步骤的所有阶段跳过不平滑的凸起。

有两种策略可以避免这种情况

将 max_step 设置为 2*tau 以便肯定会遇到凹凸，步长选择将确保跳跃周围的采样是密集的。
对跳转的片段使用单独的积分，这样每个片段内的控制输入都是一个常数。

第一个变体在某种程度上是对步长控制器的滥用，因为它在跳跃周围的一些启发式紧急模式下在规范之外工作。第二种变体需要更多的编码工作，但内部步骤较少，因为每个段的 ODE 函数再次完全平滑。