我如何对 Python 中洛伦兹和余弦乘积的函数进行数值积分?
How do I numerically integrate a function thats a product of a lorentzian and a cosinus in Python?
我是 Whosebug 的新手,也是 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 函数更适合这种函数的数值积分方法,还是我犯了更基本的错误?
观察
接近于零(1 - cos(w*t)) / w**2
趋于0/0
。我们可以采用泰勒展开 t**2(1/2 - (w*t)**2/24).
走向无穷大洛伦兹是常数,余弦项会导致输出无限振荡,积分可以通过用该项乘以一个缓慢递减的项来近似。
您使用的是具有许多点的线性间隔刻度。在对数刻度中使用 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
,说明阻尼对结果影响不大
我是 Whosebug 的新手,也是 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 函数更适合这种函数的数值积分方法,还是我犯了更基本的错误?
观察
接近于零
(1 - cos(w*t)) / w**2
趋于0/0
。我们可以采用泰勒展开 t**2(1/2 - (w*t)**2/24).走向无穷大洛伦兹是常数,余弦项会导致输出无限振荡,积分可以通过用该项乘以一个缓慢递减的项来近似。
您使用的是具有许多点的线性间隔刻度。在对数刻度中使用
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
,说明阻尼对结果影响不大