慢收敛积分的数值积分

Question

基本上，这种类型的数值积分有一个问题：

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad

# Parameter
T = 94
Eps_inf = 3.05
Eps_s = 19.184544603857724
tau_0 = 1.27*10**-16
tau_CD = 0.34580390675331274
gamma = 0.49
C_pek = 1/Eps_inf - 1/Eps_s
hbar = 6.582119569*10**-16

# Functions
def func(w):
    return -4 * 1/(np.pi * C_pek) * Eps(w).imag/abs(Eps(w))**2

def Eps(w):
    return Eps_inf + (Eps_s - Eps_inf)/(1 + (1j * w * tau_CD))**gamma

w = np.logspace(-9,80,100000)
y = func(w)

IntegrandS = quad(lambda w: func(w),0,np.inf,limit=100000)[0]
print(f'quadResult: {IntegrandS}')

这给了我一个警告:积分可能是发散的，或者是缓慢收敛的。

这一功能确实在缓慢地收敛：

​

如果我在积分的上限上加上一个像1e15这样的大数字，它会给我一个结果，但是这个结果永远不会收敛到越来越高的集成极限。

是否有处理这个函数的方法，以便quad函数(或任何其他的集成方法，我也尝试了trapz，给我同样的问题)可以处理这个函数？

谢谢!

Answer 1

积分是发散的。这可以通过研究func(w)作为w→∞的渐近行为来解释。

func(w)中依赖于w的部分是商Eps(w).imag/abs(Eps(w))**2。

当w→∞时，Eps(w)的真实部分接近Eps_inf，虚部变为0，w→∞，abs(Eps(w))**2→Eps_inf**2也是如此。也就是说，分母接近一个正常数。所以这个商的重要部分是分子，Eps(w).imag。

当w→∞时，Eps(w).imag收敛到0，但这并不意味着func(w)的积分将收敛。为了收敛，Eps(w).imag必须收敛到0“足够快”。作为w→∞，Eps(w).imag的行为类似于D * w**-gamma，其中D是一个独立于w的常数。当p<-1时，x**p形式从x0到∞的积分(对于某些x0 > 0)是收敛的。用你的函数来说，这个功率是-gamma，是-0.49。所以你的函数衰减得太慢，以至于积分不能收敛。

如果将gamma更改为大于1的值，则可以检查quad是否收敛。例如。

In [65]: gamma = 1.49

In [66]: quad(func, 0, np.inf)
Out[66]: (28.792185960802843, 3.501437717545741e-07)

In [67]: gamma = 1.2

In [68]: quad(func, 0, np.inf)
Out[68]: (87.82367385721193, 4.4632464835103747e-07)

In [69]: gamma = 1.01

In [70]: quad(func, 0, np.inf)
Out[70]: (2274.435541035491, 2.4941527954069898e-06)

当伽玛为1时，积分是发散的。这里有两次使用quad的尝试。

In [71]: gamma = 1.0

In [72]: quad(func, 0, np.inf)
<ipython-input-72-b9978f1fcdcd>:1: IntegrationWarning: The maximum number of subdivisions (50) has been achieved.
  If increasing the limit yields no improvement it is advised to analyze 
  the integrand in order to determine the difficulties.  If the position of a 
  local difficulty can be determined (singularity, discontinuity) one will 
  probably gain from splitting up the interval and calling the integrator 
  on the subranges.  Perhaps a special-purpose integrator should be used.
  quad(func, 0, np.inf)
Out[72]: (882.2704810872806, 188.89503399200746)

In [73]: quad(func, 0, np.inf, limit=100000)
Out[73]: (inf, inf)

当gamma <1时，quad报告(正确地)积分可能是发散的。

In [74]: gamma = 0.98

In [75]: quad(func, 0, np.inf, limit=100000)
<ipython-input-75-005c1a83e644>:1: IntegrationWarning: The integral is probably divergent, or slowly convergent.
  quad(func, 0, np.inf, limit=100000)
Out[75]: (-1202.853690120471, 1.0195566801485256e-05)