Python

Question

我正在求解系统与脉冲相互作用时的动力学，这基本上是在求解 time-dependent 微分方程。一般来说它工作正常，但每当我将脉冲的带宽取小，即大约统一时，求解器取决于脉冲从 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,0,a0,0],\
                [0,-1,0,-a0],\
                [a0,0,-1,0],\
                [0,-a0,0,-1]])*kappa
    Dr=np.array([[1,1j,0,0],\
                [1j,1,0,0],\
                [0,0,1,1j],\
                [0,0,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,1j,0,0],\
                [1j,1,0,0],\
                [0,0,1,1j],\
                [0,0,1j,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 开始，我不完全明白为什么。不应该，不是那样的。

是求解器的问题？如果是这样，我该怎么做才能解决它？我现在担心这个问题会出现在我运行但我没有意识到的其他模拟中。

谢谢，琼.

Answer 1

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

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

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

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

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

Python - time-dependent 系数的微分方程求解器针对不同的时间偏移给出不同的动态

Python - Differential equation solver for time-dependent coefficients gives different dynamics for different time offsets

numpy

scipy

solver

differential-equations