【发布时间】：2018-11-28 15:12:30
【问题描述】：
我正在尝试使用带有 zvode 积分器的 scipy.ode 在 python 中解决耦合复杂 ODE 的系统。但是一旦我运行代码，就会出现这个错误消息。
ZVODE-- At T(=R1) and step size H(=R2), the error test failed repeatedly or with abs(H) = HMIN. In above, R1 = 0.1018805870139D-15 R2 = 0.2392554739952D-22
我确实查看了 FORTRAN 源代码，但无法弄清楚它的含义。
感谢您对此提供任何帮助。
编辑：代码已包含在内。 我还尝试打印出一些值，并为使用简单 Euler 方法的集成编写了单独的代码。从这些我有一种感觉，错误可能是由于值超出范围，即大于 10^308。 （可能是由于某些参数的错误）。谁能确认一下？
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import ode
# Constants
h_cross = 6.5821e-16
charge = 1
a = 5.65e-10
gamma = 1/50e-15
E_g = 1.43
temp = 300
k_b = 8.6173e-5
k = np.linspace(-1, 1, 32) * np.pi / a
k = k[1:]
grad_k = k[1] - k[0]
delta_c_1 = 6.9
delta_v = 1
e_c_1 = delta_c_1/2.0 * (1 - np.cos(np.abs(k * a))) + E_g/2.0
e_v = -delta_v/2.0 * (1 - np.cos(np.abs(k * a))) - E_g/2.0
t_0 = 1e-14 # Initial time
dt = 5e-15 # Time interval
t_f = t_0 + 50e-14 # Final time
t_mid = (t_0 + t_f) / 2.0
steps = int((t_f - t_0) / dt)
t = np.linspace(t_0,t_f,steps)
d = np.ones(k.size) * 3.336e-30
w_0 = 0.1 / h_cross
f_0 = w_0 / (2*np.pi)
y_0 = np.zeros([k.size,3],dtype = complex)
y_input = y_0.flatten()
solution = y_input # Inserting initial condition as the first entry in the solution
def E(times):
pulse = np.cos(w_0 * times )
fwhm = 10e-14
sigma = fwhm / 2.35
envelope = ( 1 / (2 * np.pi *sigma**2)**0.5 ) * np.exp( -((times-t_mid)/sigma)**2 / 2.0)
waveform = pulse * envelope
return waveform
# NORMALIZE VALUE OF E(t) USING VALUE AT PEAK VALUE OF E(t)
E_peak_req = 1e8
E_peak = E(t).max()
normalisation = np.abs(E_peak_req / E_peak) * (1/1.6e-19)
def dynamics(t,y):
dydt = np.zeros([k.size,3],dtype = complex)
if(solution.size == (k.size*3)):
prev_y = solution
prev_y = np.reshape(prev_y,(k.size,3))
prev_prev_y = prev_y
prev_prev_y = np.reshape(prev_prev_y,(k.size,3))
else:
last_step = solution.shape[0] - 1
prev_y = solution[last_step,:]
prev_y = np.reshape(prev_y,(k.size,3)) # Extracting the latest values of the density matrix elements obtained in the last time step
prev_prev_y = solution[last_step - 1, :]
prev_prev_y = np.reshape(prev_prev_y,(k.size,3))
for index in range(k.size):
grad_p = prev_y[index][0] - prev_prev_y[index][0]
grad_f_c = prev_y[index][1] - prev_prev_y[index][1]
grad_f_v = prev_y[index][2] - prev_prev_y[index][2]
dipole_contr = d[index] * (E(t) * normalisation)
grad_contr_1 = 1j * charge * (E(t) * normalisation) * grad_p / grad_k
grad_contr_2 = charge * (E(t) * normalisation) * grad_f_c / grad_k
grad_contr_3 = charge * (E(t) * normalisation) * grad_f_v / grad_k
dpdt = (-1j / h_cross) * ( (e_c_1[index] + e_v[index] - 1j*h_cross*gamma) * prev_y[index][0] - (1 - prev_y[index][1] - prev_y[index][2]) * dipole_contr + grad_contr_1 )
dfcdt = (1 / h_cross) * ( -2 * np.imag( dipole_contr * np.conjugate(prev_y[index][0]) ) + grad_contr_2 )
dfvdt = (1 / h_cross) * ( -2 * np.imag( dipole_contr * np.conjugate(prev_y[index][0]) ) + grad_contr_3 )
dydt[index] = np.array([dpdt, dfcdt, dfvdt])
dydt = dydt.flatten()
return dydt
solver = ode(dynamics, jac = None).set_integrator('zvode', method ='bdf')
solver.set_initial_value(y_input, t_0) #.set_f_params()
while (solver.successful() and solver.t + dt <= t_f):
solver.integrate(solver.t + dt)
solution = np.vstack((solution,solver.y))
sol = np.reshape(solution,(solution.shape[0],k.size,3))
【问题讨论】：
-
如果您提供minimal, complete and verifiable example，别人会更容易帮助您。没有它，你将得到的只是猜测。它们可能是非常好的猜测，但如果你展示代码，你会得到更好的答案。
-
您的方程式中是否存在不连续性？您正在使用
ZVODE；你的方程是analytic吗？请参阅odedocstring 中“zvode”部分下的注释。 -
@WarrenWeckesser 感谢您的建议。我已经包含了代码和我在此期间发现的一些额外细节。