Scipy.ode "Ошибка проверки неоднократно"

Question

Я пытаюсь решить систему связанных сложных ODE в Python, используя scipy.ode с интегратором zvode. Но это сообщение об ошибке появляется, когда я запускаю код.

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

Это означает, что ваша система жесткая, эвристика контроллера размера шага вычисляет, что ему нужны очень маленькие размеры шага, чтобы гарантировать требуемые границы ошибок, но размер шага стал настолько маленьким и, таким образом, количество требуемых шагов настолько большим, что накопление шума с плавающей точкой становится более доминирующим, что означает, что контроллер теряет контроль над накоплением ошибок. Похоже, что во избежание этого контроллер ограничивает размер шага примерно 2e-7чуть больше чем sqrt(mu), части стоимости t,