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

Question

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

3

我是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中，是否有更适合这种函数的数值积分的不同方法，比quad函数更好，或者我犯了更基本的错误？

- Frederik

2

如果你将积分上限设为100而不是np.inf，积分结果会更快，并且与np.inf的结果相差不大。不确定原因，但可能是无穷大转换的检测对你的函数效果不佳。我会把这个问题留给专家们来解决。 - hesham_EE

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Bob · Accepted Answer

观察：

接近零时，(1-cos(w*t))/w**2 趋向于 0/0。我们可以采用泰勒展开式 t**2(1/2-(w*t)**2/24)。
当趋近于无穷大时，Lorentzian是一个常数，余弦项会导致输出无限振荡，积分可以通过将该项乘以缓慢下降的项来近似。
您正在使用具有许多点的线性间隔刻度。使用对数刻度的 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，这表明阻尼对结果的影响不太大。