[Python] Integration limits when integrating to infinity using quad and vectorize

I have a rather complicated physics-model (Thomas Hopfield model from https://www.researchgate.net/public..._at_Randomly_Distributed_Donors_and_Acceptors).

For the sake of later fitting it to experimental data, I first simply want to graphically plot it and am already running into difficulties. The model reacts very sensible to changes in the integration limit and I am not sure, whether the vectorization-function creates problems with the integration.

The model consists of two analytically not-solvable integrals from 0 to infinity, and I want to plot them using the scipy.integrate.quad function.

(example code below)

As expected, integrating to infinity does not work, so I integrated up to a set value (r_max). However, the resulting function is very sensible to changes in the r_max. When increasing it (e.g. from 250e-9 to 250e-8), the resulting function gets stretched. That poses a problem, as it seems unjustified to then take the highest computable r_max and assume it similar to infinity. (What is even weirder, is that there is a theoretical value for t=0, which is most closely produced using r_max = 250e-9 and not using higher values "closer to infinity".)

#### Script for numerically integrating and plotting the Thomas-Hopfield-Model


#-------------- IMPORTING LIBRARIES --------------
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad


#-------------- DEFINING MODEL FUNCTION ----------------
# The physical model has three parameters (W_0, R_0, n),
# in order to fit it onto experimental data later, two constants (A, B) have been added.
# it contains two integrals over the variable r from 0 to infinity, but due to overflow issues have here been replaced by integrals from 0 to a set value r_max.
# The overall model is a function of time t.

def W(r,W_0,R_0): # exponential decay function for the recombination rate of electron-hole pair
return W_0 * np.exp( (-2 * r) / R_0 )

def integral1(r,W_0,R_0,t): # First integral (over r) for the Thomas-Hopfield model
return W(r,W_0,R_0) * np.exp( -W(r,W_0,R_0) * t ) * r**2

def integral2(r,W_0,R_0,t): # second integral (over r) for the Thomas-Hopfield model (in the exponent)
return ( np.exp( -W(r,W_0,R_0) *t ) - 1 ) *r**2

def thomas_hopfield(t,W_0,R_0,n,A,B): # comined function for the Thomas-Hopfield model
r_max = 250e-7 #np.inf #UPPER INTEGRATION BOUND
int1 = quad( integral1, 0, r_max, args=(W_0,R_0,t) )[0]
int2 = quad( integral2, 0, r_max, args=(W_0,R_0,t) )[0]
lum_int = A*( (4 * np.pi * n * int1) * np.exp(4 * np.pi * n * int2) ) + B
return lum_int

thomas_hopfield_v = np.vectorize(thomas_hopfield) # this has been added, as trying to plot it without vectorization resulted in the error: "TypeError: only size-1 arrays can be converted to Python scalars"


#--------------- EXAMPLE VALUES FOR PARAMETERS ---------------:
W_0 = 1e6
R_0 = 1.2e-9
n = 3e-19
A = 1
B = 0


#---------------------------------------------------------------
#--------- SIDE QUEST: theoretical value for t=0 --------------

theoretical_value = n*np.pi*W_0*R_0**3 # literature theoretical value for t=0
thomas_hopfield_t0 = thomas_hopfield(0,W_0,R_0,n,A,B) # calculate the theoretical value for t=0
print('Theoretical value for t=0:', theoretical_value)
print('Calculated value for t=0:', thomas_hopfield_t0) # print the calculated value for t=0

#-----------------------------------------------------------------
#----------------------------------------------------------------

#--------------- GRAPHICAL REPRESENTATION-------------

t = np.linspace(0,100000000, 100000) # create time value array
y_model = thomas_hopfield_v(t,W_0,R_0,n,A,B) # apply the model on the time t
plt.plot(t, y_model, color='r',linestyle='-',linewidth=0.7, alpha=1, label='Thomas Hopfield Model',zorder=1)
plt.yscale('log')
plt.xscale('log')
plt.show()


As someone using scipy.integrate.quad for the first time, I am asking myself how quad is computing the result (in which order it takes in the array, computes the integrals and returns the values) and whether or not the vectorization np.vectorize() of the function is the (or a) problem. I had to include np.vectorize(), as I would otherwise receive "TypeError: only size-1 arrays can be converted to Python scalars", but it makes it almost impossible for me to follow and comprehend the individual steps of the programme.

Does someone have experience with that? It would already help me to know, whether the vectorisation creates the problems and / or how I could plot the function without using it. Thank you very much in advance!

Continue reading...