Ограничения интеграции при интеграции в бесконечность с использованием Quad и Vectorize ⇐ Python

[Anonymous]

У меня довольно сложная физика-модель (модель Thomas Hopfield от https://www.researchgate.net/publicatio ... _accecters) Сначала я просто хочу графически построить его и уже сталкиваюсь с трудностями. Модель реагирует очень разумно на изменения в пределе интеграции, и я не уверен, создает ли векторизационная функция проблемы с интеграцией. < /P>
Модель состоит из двух аналитических неразрешаемых интегралов от 0 до бесконечности, и я хочу построить их с использованием код SCIPY.Integrate.quad. Ожидается, что интеграция в бесконечность не работает, поэтому я интегрировал до установленного значения (

Код: Выделить всё

r_max

). Однако полученная функция очень разумна для изменений в R_MAX . При его увеличении (например, с 250e-9 до 250e-8), полученная функция растягивается. Это создает проблему, так как кажется неоправданным, а затем взять самый высокий вычисляемый R_MAX и предположить, что это похоже на бесконечность. (Что даже более страннее, так это то, что существует теоретическое значение для t = 0 , которое наиболее тесно создается с использованием R_MAX = 250E-9 и не использует более высокие значения «ближе к бесконечности».)

Код: Выделить всё

#### 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()

Как кто -то, использующий scipy.integrate.quad Впервые я спрашиваю себя, как Quad вычисляет результат (в каком порядке он принимает в массиве, вычисляет интегралы и возвращает значения) и независимо от того, является ли векторизация np.Vectorize () функции (или (или проблема. Я должен был включить np.vectorize () , как я бы в противном случае получил «TypeError: только массивы размера 1 могут быть преобразованы в скаларс Python», но для меня практически невозможно следовать и понимать отдельные шаги программы. Это уже поможет мне узнать, создает ли векторизация проблемы и / или как я мог бы построить функцию, не используя ее. Большое спасибо заранее!

Подробнее здесь: https://stackoverflow.com/questions/797 ... -vectorize