python - python 中的振荡积分

Question

我编写了以下代码来绘制从光学组件射出的光的强度，它基本上是入射场上的球面傅立叶积分，因此它具有贝塞尔函数。其参数取决于积分变量 (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 如何处理这个问题？对于这个特定的应用，是否有更好的数值积分方法？

Answer 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)