【发布时间】：2014-11-21 16:24:40
【问题描述】：
我试图模拟强激光的振荡电场如何推动接近 +1 离子库仑势的电子。激光场是
E = Eo sin(wt)，在 y 方向。
库仑势是
F = ke q1*q2/r^2，在 r 方向。
强电场使电子向tunnel ionize，所以电子的初始条件是在y方向上从原子上位移。然后，电子被激光场来回推动，并有机会与库仑势相互作用。我想模拟库仑势如何影响电子的飞行。模拟需要在三个维度上进行，因为我最终想要包括更复杂的激光场，这些激光场在 x 和 y 方向上推动电子，并且电子可以从 z 方向上的动量开始。
起初，我认为这很容易。下面是我用来以非常小的步骤（1e-18 秒）逐步完成时间的代码。当电子不在离子附近时，这可以正常工作。然而，对于通过靠近离子的电子，结果很大程度上取决于模拟中使用的时间步长。如果我使时间步长变小，则计算需要很长时间。
所以，我认为在这种情况下，我应该使用自适应时间步长。此外，根据我的阅读，龙格-库塔方法应该优于我使用的简单方法。但是，我不认为 scipy.odeint 适用于三维。关于如何提高这些模拟的准确性和速度的任何想法？
下图显示了时间步长如何对结果产生巨大影响（一件坏事）：
这是我的代码：
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()
【问题讨论】：
-
有趣的问题；但几乎不是编程问题...
-
问题基本上是我如何在python中实现自适应时间步长runge-kutta。可能已经有一个 python 包可以为我做这件事，所以这是一个编程问题。我只是提供了太多背景信息。
-
“但是，我认为 scipy.odeint 不适用于三维。”这不是真的。
scipy.odeint设计用于求解 n 维系统。例如，请参阅 wiki.scipy.org/Cookbook/CoupledSpringMassSystem 或 wiki.scipy.org/Cookbook/KdV -
要使用
scipy.odeint，您的首要任务是将微分方程编写为一阶系统（就像我之前评论中链接的 SciPy Cookbook 中的一对弹簧质量示例一样）。在这种情况下，您的系统有六个维度；变量是 x、y、z、vx、vy 和 vz。将系统的“右侧”作为函数实现后，您可以使用scipy.odeint。
标签： python scipy physics ode differential-equations