我正在计算电子-正电子对产生的狄拉克方程的自由基展开的解。为此，我需要解一个方程组，看起来像这样:
配对生产方程，出自Mocken at al.
编辑:通过将y0作为复杂类型传递到求解器中已经解决了这个问题。正如在这个问题中所述:https://github.com/scipy/scipy/issues/8453我肯定会认为这是一个bug，但它似乎已经在岩石下至少4年了
我使用SciPy的solve_ivp积分器，方法如下:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
from scipy.integrate import solve_ivp
import scipy.constants as constants
#Impulse
px, py = 0 , 0
#physics constants
e = constants.e
m = constants.m_e # electronmass
c = constants.c
hbar = constants.hbar
#relativistic energy
E = np.sqrt(m**2 *c**4 + (px**2+py**2) * c**2) # E_p
#adiabatic parameter
xi = 1
#Parameter of the system
w = 0.840 #frequency in 1/m_e
N = 8 # amount of amplitudes in window
T = 2* np.pi/w
#unit system
c = 1
hbar = 1
m = 1
#strength of electric field
E_0 = xi*m*c*w/e
print(E_0)
#vectorpotential
A = lambda t,F: -E_0/w *np.sin(t)*F
def linearFenster2(t):
conditions = [t <=0, (t/w>=0) and (t/w <= T/2), (t/w >= T/2) and (t/w<=T*(N+1/2)), (t/w>=T*(N+1/2)) and (t/w<=T*(N+1)), t/w>=T*(N+1)]
funcs = [lambda t: 0, lambda t: 1/np.pi *t, lambda t: 1, lambda t: 1-w/np.pi * (t/w-T*(N+1/2)), lambda t: 0]
return np.piecewise(t,conditions,funcs)
#Coefficient functions
nu = lambda t: -1j/hbar *e*A(w*t,linearFenster2(w*t)) *np.exp(2*1j/hbar * E*t) *(px*py*c**2 /(E*(E+m*c**2)) + 1j*(1- c**2 *py**2/(E*(E+m*c**2))))
kappa = lambda t: 1j*e*A(t,linearFenster2(w*t))* c*py/(E * hbar)
#System to solve
def System(t, y, nu, kappa):
df = kappa(t) *y[0] + nu(t) * y[1]
dg = -np.conjugate(nu(t)) * y[0] + np.conjugate(kappa(t))*y[1]
return np.array([df,dg], dtype=np.cdouble)
def solver(tmin, tmax,teval=None,f0=0,g0=1):
'''solves the system.
@tmin: starttime
@tmax: endtime
@f0: starting percentage of already present electrons of positive energy usually 0
@g0: starting percentage of already present electrons of negative energy, usually 1, therefore full vaccuum
'''
y0=[f0,g0]
tspan = np.array([tmin, tmax])
koeff = np.array([nu,kappa])
sol = solve_ivp(System,tspan,y0,t_eval= teval,args=koeff)
return sol
#Plotting of windowfunction
amount = 10**2
t = np.arange(0, T*(N+1), 1/amount)
vlinearFenster2 = np.array([linearFenster2(w*a) for a in t ], dtype = float)
fig3, ax3 = plt.subplots(1,1,figsize=[24,8])
ax3.plot(t,E_0/w * vlinearFenster2)
ax3.plot(t,A(w*t,vlinearFenster2))
ax3.plot(t,-E_0 /w * vlinearFenster2)
ax3.xaxis.set_minor_locator(ticker.AutoMinorLocator())
ax3.set_xlabel("t in s")
ax3.grid(which = 'both')
plt.show()
sol = solver(0, 70,teval = t)
ts= sol.t
f=sol.y[0]
fsquared = 2* np.absolute(f)**2
plt.plot(ts,fsquared)
plt.show()
窗口函数的图看起来像这样(并且是正确的)窗口函数然而，解的图是这样的:配对生产概率图根据论文的图表(以及使用mathematica进行的进一步测试)，这是不正确的。
当运行'sol = solver(..)'这行时，它说:
numpycore_asarray.py:102: ComplexWarning: Casting complex values to real discards the imaginary part
return array(a, dtype, copy=False, order=order)
我简直不知道为什么solve_ivp抛弃了虚部。这是绝对必要的。
有谁能告诉我谁知道的更多或看到的错误?
根据文献，传递给solve_ivp的y0必须是complex类型，以便在复域上进行积分。确保这一点的一个健壮的方法是在代码中添加以下内容:
def solver(tmin, tmax,teval=None,f0=0,g0=1):
'''solves the system.
@tmin: starttime
@tmax: endtime
@f0: starting percentage of already present electrons of positive energy usually 0
@g0: starting percentage of already present electrons of negative energy, usually 1, therefore full vaccuum
'''
f0 = complex(f0) # <-- added
g0 = complex(g0) # <-- added
y0=[f0,g0]
tspan = np.array([tmin, tmax])
koeff = np.array([nu,kappa])
sol = solve_ivp(System,tspan,y0,t_eval= teval,args=koeff)
return sol