在python中用自适应步长模拟电子运动微分方程

Question

我正试图模拟一个强激光的振荡电场如何推动一个接近a+1离子库仑势的电子。激光场是

E=Eo sin（wt），在y方向。在

库仑势是

F=ke q1*q2/r^2，在r方向。在

强电场使电子tunnel ionize，因此电子的初始条件是在y方向上从原子中移位。然后，电子被激光场来回推动，有机会与库仑势相互作用。我想模拟库仑势如何影响电子的飞行。模拟必须是三维的，因为我最终想要包括更复杂的激光场，它能在x和y方向推动电子，电子可以从z方向的动量开始。在

一开始，我以为这很容易。下面是我用来以非常小的步骤（1e-18秒）逐步穿越时间的代码。当电子不在离子附近时，这是正常的。然而，对于通过离子附近的电子，结果强烈依赖于模拟中使用的时间步长。如果我把时间步长变小，计算要花很长时间。在

所以，我认为在这个例子中我应该使用一个自适应的时间步长。另外，从我所读到的内容来看，Runge-Kutta方法应该优于我使用的简单方法。不过，我不这么认为西皮·奥德林适用于三维。关于如何提高这些模拟的准确性和速度有什么想法吗？在

下图显示了时间步长如何对结果产生巨大影响（一件坏事）：

我的密码是：

import numpy as np
import matplotlib.pyplot as plt

q = 1.602e-19       #Coulombs   Charge of electron
h_bar = 1.054e-34   #J*s        Plank's Constant div by 2Pi
c = 3.0e8           #m/s        Speed of light
eo = 8.8541e-12     #C^2/(Nm^2) Permittivity of vacuum
me = 9.109e-31      #kg         Mass of electron
ke = 8.985551e9     #N m^2 C-2  Coulomb's constant

def fly_trajectory(wavelength,intensity,tb=0,pulseFWHM=40.e-15,
                   final_time=100e-15,time_step=.001e-15,Ip=15.13,v0=(2e4,0,0)):
    #Intensity is in w/cm2. Ip is in eV. Otherwise it's SI units throughout.
    #The electric field of the later is in the y-direction

    Ip = 15.13 * q  #Calculate the ionization potential of the atom in Joules
    Eo = np.sqrt(2*intensity*10**4/(c*8.85e-12)) # Electric field in V/m
    w = c/wavelength * 2. * np.pi #Angular frequency of the laser

    times = np.arange(tb,final_time,time_step)
    Ey = Eo*np.sin(w*times) * np.exp(-times**2/(2*(pulseFWHM / 2.35482)**2)) 
    Eb = Ey[0] #E-field at time of birth (time of tunneling)

    if Eb == 0: return 0,0 #No field --> no electrons
    tunnel_position = -Ip / (Eb*q)    
    x,y,z = 0,tunnel_position,0
    vx,vy,vz = v0

    y_list = np.zeros(len(times)) #store the y-values for later

    for index in range(0,len(times)):
        r = np.sqrt(x**2+y**2+z**2)
        rx = x/r; ry = y/r; rz=z/r

        Fcy = -q**2 * ke/(r**2) * ry
        ay = Ey[index]*(-q)/me + Fcy/me #only y includes the laser
        vy = vy + ay*time_step 
        y = y + vy * time_step

        Fcx = -q**2 * ke/(r**2) * rx
        ax = (-q)/me + Fcx/me
        vx = vx + ax*time_step  
        x = x + vx * time_step

        Fcz = -q**2 * ke/(r**2) * rz
        az = (-q)/me + Fcz/me
        vz = vz + az*time_step 
        z = z + vz * time_step

        y_list[index] = y

    return times,y_list

for tb in np.linspace(0.25*2.66e-15,0.5*2.66e-15,5):
    print tb
    times,ys = fly_trajectory(800e-9,2e14,tb=tb,time_step=.01e-15)
    plt.plot(times,ys,color='r')

    times,ys = fly_trajectory(800e-9,2e14,tb=tb,time_step=.001e-15)
    plt.plot(times,ys,color='b')


#plot legend and labels:
plt.plot(0,0,color='r',label='10e-18 sec step')
plt.plot(0,0,color='b',label='1e-18 sec step')
plt.xlim(0,10e-15); plt.ylim(-1e-8,1e-8)
leg = plt.legend(); leg.draw_frame(False)
plt.xlabel('Time (sec)')
plt.ylabel('Y-distance (meters)')
plt.show()

Answer 1

正如Warren Weckesser建议的，我可以简单地遵循耦合质量-弹簧系统的Scipy食谱。首先，我需要将“右侧”方程式写成：

x'  = vx
y'  = vy
z'  = vz
vx' = Ac*x/r
vy' = Ac*y/r + q*E/m
vz' = Ac*z/r 

式中Ac=keq^2/（mr^2）是库仑势引起的加速度值，E是激光的随时间变化的电场。然后，我可以使用scipy.integrate.odeint来找到解决方案。这比我以前使用的方法更快、更可靠。在

这是用odeint得到的电子轨迹。现在没有一个疯狂地飞走了：

代码如下：

^{pr2}$