python中的振荡积分答案

【问题标题】：Oscillatory integral in pythonpython中的振荡积分
【发布时间】：2018-08-09 20:52:02
【问题描述】：

我编写了以下代码来绘制出光学组件的光强度，这基本上是入射场上的球面傅里叶积分，因此它具有贝塞尔函数。其参数取决于积分变量（x）和绘图变量（r）。

from sympy import *
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import quad
from scipy.special import jn

#constants
mm = 1
um = 1e-3 * mm
nm = 1e-6 * mm

wavelength = 532*nm
klaser = 2*np.pi / wavelength
waist = 3.2*mm
angle = 2 #degrees
focus = 125 * mm
ng = 1.5 # refractive index of axicon
upperintegration = 5

#integrals

def b(angle):
    radians = angle* np.pi/180
    return klaser * (ng-1) * np.tan(radians)

def delta(angle):
    return np.pi/(b(angle)*waist)    

def integrand(x, r): 
    return klaser/focus * waist**2 * np.exp(-x**2) * np.exp(-np.pi * 1j * x/delta(angle)) * jn(0, waist*klaser*r*x/focus) * x

def intensity1D(r):
    return np.sqrt(quad(lambda x: np.real(integrand(x, r)), 0, upperintegration)[0]**2 + quad(lambda x: np.imag(integrand(x, r)), 0, upperintegration)[0]**2)  

fig = plt.figure()
ax = fig.add_subplot(111)

t = np.linspace(-3.5, 3.5, 25)
plt.plot(t, np.vectorize(intensity1D)(t))

问题是，当我改变我在linspace 中使用的点数时，情节发生了巨大变化。

我怀疑这可能是因为积分的振荡性质，因此所采用的步长可以显着改变指数的值，从而改变积分的值。

quad 如何处理这个问题？有没有更好的方法来为这个特定的应用程序进行数值积分？

【问题讨论】：

您确定公式正确吗？通常你不会有像 mm 这样的单位出现在 exp 或 jn 函数中。（例如，平方毫米的指数是多少？？）
X 是无量纲的

标签： python python-2.7 matplotlib scipy

【解决方案1】：

在对quad 的调用中，将limit 参数设置为一个大数字。这增加了允许quad 用于估计积分的最大子区间数。当我使用

def intensity1D(r):
    re = quad(lambda x: np.real(integrand(x, r)), 0, upperintegration, limit=8000)[0]
    im = quad(lambda x: np.imag(integrand(x, r)), 0, upperintegration, limit=8000)[0]
    return np.sqrt(re**2 + im**2)

计算函数，数组t定义为

t = np.linspace(1.5, 3, 1000)

我得到以下情节：

(我还删除了from sympy import *这一行。sympy似乎没有用于你的脚本。）

您应该始终检查作为quad 的第二个返回值的错误估计。例如：

In [14]: r = 3.0

In [15]: val, err = quad(lambda x: np.real(integrand(x, r)), 0, upperintegration, limit=8000)

In [16]: val
Out[16]: 2.975500141416676e-11

In [17]: err
Out[17]: 1.4590630152807049e-08

如您所见，误差估计远大于近似积分。 quad 返回的估计值可能是保守的，但仍应谨慎对待具有如此大误差估计值的结果。我们来看看对应的虚部：

In [25]: val, err = quad(lambda x: np.imag(integrand(x, r)), 0, upperintegration, limit=8000)

In [26]: val
Out[26]: 0.0026492702707317257

In [27]: err
Out[27]: 1.4808416189183e-08

val 现在比估计的误差大几个数量级。因此，当在intensity1D() 中计算复数值的大小时，我们最终得到的估计相对误差约为 1e-5。这对您的计算可能已经足够了。

在 r=2.1825 附近的峰值处，误差估计的幅度仍然很小，并且远小于计算的积分：

In [32]: r = 2.1825

In [33]: quad(lambda x: np.real(integrand(x, r)), 0, upperintegration, limit=8000)
Out[33]: (6.435730031424414, 8.801375195176556e-08)

In [34]: quad(lambda x: np.imag(integrand(x, r)), 0, upperintegration, limit=8000)
Out[34]: (-6.583055286038913, 9.211333259956749e-08)

【讨论】：