慢速收敛积分的数值积分

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函数（或者其他任何积分方法，我也试过scipys 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 的常数。从 x0 到 ∞（对于某些 x0 > 0）的 x**p 形式的积分仅在 p < -1 时收敛。使用您的函数，该幂为 -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)

当gamma为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)

慢速收敛积分的数值积分

Numerical integration for a slow converging integral

python

numpy

scipy

numerical-integration