python - 如何在 Python 中对作为洛伦兹和余弦乘积的函数进行数值积分？

Question

我是 stackoverflow 的新手，也是 Python 的新手。所以，我希望以适当的方式提出我的问题。我正在运行一个类似于这个最小示例的 Python 代码，其中包含一个示例函数，该函数是洛伦兹函数的产物，我想对其进行数值积分:

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

#minimal example:
omega_loc = 15
gamma = 5

def Lorentzian(w):
    #print(w)
    return (w**3)/((w/omega_loc) + 1)**2*(gamma/2)/((w-omega_loc)**2+(gamma/2)**2)

def intRe(t):
    return quad(lambda w: w**(-2)*Lorentzian(w)*(1-np.cos(w*t)),0,np.inf,limit=10000)[0]


plt.figure(1)
plot_range = np.linspace(0,100,1000)
plt.plot(plot_range, [intRe(t) for t in plot_range])

独立于集成的上限，我从来没有让代码运行并给我一个结果。当我启用 #print(w) 行时，代码似乎只是在无限循环(？)中继续探测 w 的随机不同值的积分。控制台还让我检测到舍入错误。在 Python 中是否有一种不同的数值积分方法比四元函数更适合这种函数，还是我犯了更基本的错误？

Answer 1

观察

1. 接近于零 (1 - cos(w*t))/w**2 趋向于 0/0。我们可以采用泰勒展开 t**2(1/2 - (w*t)**2/24)。

2. 接近无穷大时，洛伦兹是常数，余弦项将导致输出无限振荡，积分可以通过将该项乘以缓慢下降的项来近似。

3. 您正在使用具有许多点的线性间隔比例尺。在对数刻度中使用 w 更容易可视化。

在对余弦项进行阻尼之前，该图看起来像这样

我引入了两个参数来调整振荡的衰减

def cosinus_term(w, t, damping=1e4*omega_loc):
    return np.where(abs(w*t) < 1e-6,  t**2*(0.5 - (w*t)**2/24.0), (1-np.exp(-abs(w/damping))*np.cos(w*t))/w**2)
def intRe(t, damping=1e4*omega_loc):
    return quad(lambda w: cosinus_term(w, t)*Lorentzian(w),0,np.inf,limit=10000)[0]

用下面的代码绘图

plt.figure(1)
plot_range = np.logspace(-3,3,100)
plt.plot(plot_range, [intRe(t, 1e2*omega_loc) for t in plot_range])
plt.plot(plot_range, [intRe(t, 1e3*omega_loc) for t in plot_range])
plt.xscale('log')

这里运行不到 3 分钟，两个结果非常接近，特别是对于大的 w，表明阻尼对结果影响不大。