4

I'm using scipy.integrate.ode and would like to know, what happens internally when I get the message UserWarning: zvode: Excess work done on this call. (Perhaps wrong MF.) 'Unexpected istate=%s' % istate)) This appears when I call ode.integrate(t1) for too big t1, so I'm forced to use a for-loop and incrementally integrate my equation, what lowers the speed since the solver can not use adaptive step size very effectively. I already tried different methods and setting for the integrator. The maximal number of steps nsteps=100000 is very big already but with this setting I still can't integrate up to 1000 in one call, which I would like to do.

The code I use is:

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)
Fred_UB
  • 49
  • 2
  • The automatically generated step size seems to be about `2e-4` or `5000` steps per time difference `1`. Thus there seems no penalty for integration in smaller increments. – Lutz Lehmann Nov 24 '16 at 19:35
  • Is `y_out[0:n_k]=-g*g*2*np.real(y_list[2*n_k:3*n_k])/h_bar` correct or should it be `y_out[0:n_k]=-g*g*2*np.real(y_list[1*n_k:2*n_k])/h_bar`, i.e., the derivative of the first third is the second third and the derivative of the second third is the third third. – Lutz Lehmann Nov 24 '16 at 19:40
  • I checked that again...The code is correct. The issue seems to be with the internal maximal step size (which was too small for my problem). – Fred_UB Nov 28 '16 at 10:44
  • Do you have a link to the equations that you implement to check against the equations I read off the code? – Lutz Lehmann Nov 28 '16 at 11:31

1 Answers1

3

The message means that the method reached the number of steps specified by nsteps parameter. Since you asked about internals, I looked into the Fortran source, which offers this explanation:

-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.)

The conditional statement that brings up the error is this "GO TO 500".

According to LutzL, for your ODE the solver chooses step size 2e-4, which means 5000000 steps to integrate up to 1000. Your options are:

  • try such a large value of nsteps (which translates to MXSTEP in aforementioned Fortran routine)
  • reduce error tolerance
  • run a for loop, as you already do.
  • I used a maximal step size of `nsteps=1e8`...that indeed solves the problem. The calculation takes very long though. I think that is due to the small but finite oscillations due to the imaginary part in one of the equations. This inhibits adaptive (large) step sizes...I think. Thanks for the help. – Fred_UB Nov 28 '16 at 10:43
  • 1
    Look into operator splitting, exponential methods of numerical integration, ... essentially, you solve the fast oscillations exactly via matrix exponentials and the slower components numerically, which should take the fast oscillations out of the step size computation. – Lutz Lehmann Nov 28 '16 at 11:05