python - Scipy.ode "Error test failed repeatedly"

Question

我正在尝试使用 scipy.ode 和 zvode 积分器来求解 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 源代码，但无法弄清楚它的含义。

任何有关此问题的帮助都将不胜感激。

编辑:代码已包含在内。 我还尝试打印出一些值，并为采用简单欧拉方法的集成编写了单独的代码。从这些我有一种感觉，错误可能是由于值超出范围，即大于 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))

Answer 1

这意味着您的系统是僵硬的，步长 Controller 的启发式计算表明它需要非常小的步长来保证所需的误差范围，但是步长已经变得如此之小，因此所需的步数如此之大浮点噪声的积累变得更加占主导地位，这意味着 Controller 失去了对误差积累的控制。似乎为了避免这种情况， Controller 将步长限制在大约 2e-7，比 sqrt(mu) 稍大，是 t 值的一部分.