scipy.integrate.ode 的内部工作

Posted 2023-03-12

【中文标题】scipy.integrate.ode 的内部工作

【英文标题】：Internal working of scipy.integrate.ode

【发布时间】：2017-04-08 21:30:18

【问题描述】：

我正在使用 scipy.integrate.ode 并想知道当我收到消息 UserWarning: zvode: Excess work done on this call. (Perhaps wrong 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)

【问题讨论】：

自动生成的步长似乎约为2e-4 或5000 每时差步数1。因此，以较小的增量集成似乎没有任何惩罚。

是y_out[0:n_k]=-g*g*2*np.real(y_list[2*n_k:3*n_k])/h_bar 正确还是应该是y_out[0:n_k]=-g*g*2*np.real(y_list[1*n_k:2*n_k])/h_bar，即第一个三分之一的导数是第二个三分之一，第二个三分之一的导数是第三个三分之一。

我再次检查了...代码是正确的。问题似乎在于内部最大步长（对于我的问题来说太小了）。

您是否有链接到您实施的方程式以检查我从代码中读取的方程式？

【参考方案1】：

该消息表示该方法已达到nsteps 参数指定的步骤数。既然你问的是内部结构，我查看了Fortran source，它提供了这样的解释：

-1 表示在完成请求的任务之前，在此调用上完成了过多的工作（超过 MXSTEP 步骤），但就 T 而言，集成是成功的。（MXSTEP 是可选输入，通常为 500 .)

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

根据 LutzL，求解器为您的 ODE 选择步长 2e-4，这意味着 5000000 步可以积分到 1000。您的选择是：

尝试使用如此大的 nsteps 值（在上述 Fortran 例程中转换为 MXSTEP） 降低容错率 运行for 循环，就像你已经做的那样。

【讨论】：

我使用了最大步长nsteps=1e8...确实解决了这个问题。不过计算时间很长。我认为这是由于方程之一中的虚部引起的小而有限的振荡。这抑制了自适应（大）步长......我认为。感谢您的帮助。

研究算子分裂，数值积分的指数方法，...本质上，您可以通过矩阵指数和数值上的较慢分量精确解决快速振荡，这应该将快速振荡排除在步长之外计算。