内部工作scipy.integrate.od

Question

我在用scipy.integrate.ode我想知道，当我收到消息UserWarning: zvode: Excess work done on this call. (Perhaps wrong MF.) 'Unexpected istate=%s' % istate))时，内部会发生什么 当我调用ode.integrate(t1)太大t1时就会出现这种情况，所以我被迫使用for-循环并逐步积分我的方程，这会降低速度，因为解算器不能非常有效地使用自适应步长。我已经为积分器尝试了不同的方法和设置。最大步骤数nsteps=100000已经非常大了，但是使用这个设置，我仍然不能在一个调用中集成多达1000步，我想这样做。在

我使用的代码是：

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import ode

h_bar=0.658212      #reduced Planck's constant (meV*ps)
m0=0.00568563       #free electron mass (meV*ps**2/nm**2)
m_e=0.067*m0        #effective electron mass (meV*ps**2/nm**2)
m_h=0.45*m0         #effective hole mass (meV*ps**2/nm**2)
m_reduced=1/((1/m_e)+(1/m_h)) #reduced mass of electron and holes combined
kB=0.08617          #Boltzmann's constant (meV/K)
mu_e=-50            #initial chemical potential for electrons
mu_h=-100           #initial chemical potential for holes
k_array=np.arange(0,1.5,0.02) #a list of different k-values
n_k=len(k_array) #number of k-values

def derivative(t,y_list,Gamma,g,kappa,k_list,n_k):
    #initialize output vector
    y_out=np.zeros(3*n_k+1,dtype=complex)

    y_out[0:n_k]=-g*g*2*np.real(y_list[2*n_k:3*n_k])/h_bar

    y_out[n_k:2*n_k]=-g*g*2*np.real(y_list[2*n_k:3*n_k])/h_bar      

    y_out[2*n_k:3*n_k]=((-1.j*(k_list**2/(2*m_reduced))-(Gamma+kappa))*y_list[2*n_k:3*n_k]-y_list[-1]*(1-y_list[n_k:2*n_k]-y_list[0:n_k])+y_list[0:n_k]*y_list[n_k:2*n_k])/h_bar

    y_out[-1]=(2*np.real(g*g*sum(y_list[2*n_k:3*n_k]))-2*kappa*y_list[-1])/h_bar
    return y_out

def dynamics(t_list,N_ini=1e-3, T=300, Gamma=1.36,kappa=0.02,g=0.095):

    #initial values
    t0=0 #initial time
    y_initial=np.zeros(3*n_k+1,dtype=complex)
    y_initial[0:n_k]=1/(1+np.exp(((h_bar*k_array)**2/(2*m_e)-mu_e)/(kB*T))) #Fermi-Dirac distributions
    y_initial[n_k:2*n_k]=1/(1+np.exp(((h_bar*k_array)**2/(2*m_h)-mu_h)/(kB*T)))

    t_list=t_list[1:] #remove t=0 from list (not feasable for integrator)

    r=ode(derivative).set_integrator('zvode',method='adams', atol=10**-6, rtol=10**-6,nsteps=100000)   #define ode solver 
    r.set_initial_value(y_initial,t0)
    r.set_f_params(Gamma,g,kappa,k_array,n_k)

    #create array for output (the +1 accounts values at t0=0)
    y_output=np.zeros((len(t_list)+1,len(y_initial)),dtype=complex)

    #insert initial data in output array
    y_output[0]=y_initial

    #perform integration for time steps given by t_list (the +1 account for the initial values already in the array)
    for i in range(len(t_list)):
        print(r't = %s' % t_list[i])
        r.integrate(t_list[i])
        if not (r.successful()):
            print('Integration not successful!!')
            break
        y_output[i+1]=r.y

    return y_output

t_list=np.arange(0,100,5)
data=dynamics(t_list,N_ini=1e-3, T=300, Gamma=1.36,kappa=0.02,g=1.095)

Answer 1

该消息表示该方法已达到nsteps参数指定的步骤数。由于您询问了内部结构，我查看了Fortran source，它提供了以下解释：

-1 means an excessive amount of work (more than MXSTEP steps) was done on this call, before completing the requested task, but the integration was otherwise successful as far as T. (MXSTEP is an optional input and is normally 500.)

引发错误的条件语句是this "GO TO 500"。在

根据LutzL，对于你的ODE，解算器选择步长2e-4，这意味着5000000个步骤可以集成到1000个。您的选择是：

• 尝试这样大的值nsteps（在前面提到的Fortran例程中，它转换成MXSTEP）
• 降低误差容限
• 运行一个for循环，就像你已经做的那样。在