python - scipy.integrate.ode 的内部工作

Question

我正在使用 scipy.integrate.ode 并且想知道，当我收到消息 UserWarning: zvode: Excess work done on this call 时内部发生了什么。 (可能是错误的 MF。)'Unexpected istate=%s' % istate))当我为太大的 t1 调用 ode.integrate(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 循环，就像您已经做的那样。