零阶汉克尔变换(Hankel Transform)
函数
f
1
(
r
)
f_1(r)
f1(r)的零阶汉克尔变换定义为:
f
2
(
r
2
)
=
2
π
∫
0
+
∞
f
1
(
r
1
)
J
0
(
2
π
r
1
r
2
)
r
1
d
r
1
f_2(r_2)=2\pi\int_0^{+\infty}f_1(r_1)J_0(2\pi r_1r_2)r_1d r_1
f2(r2)=2π∫0+∞f1(r1)J0(2πr1r2)r1dr1
其中
J
0
(
⋅
)
J_0(\cdot)
J0(⋅)为零阶贝塞尔(Bessel
)函数。其逆变换形式一样:
f
1
(
r
1
)
=
2
π
∫
0
+
∞
f
2
(
r
2
)
J
0
(
2
π
r
1
r
2
)
r
2
d
r
2
f_1(r_1)=2\pi\int_0^{+\infty}f_2(r_2)J_0(2\pi r_1r_2)r_2d r_2
f1(r1)=2π∫0+∞f2(r2)J0(2πr1r2)r2dr2
本质上0阶汉克尔变换就是柱对称体系上的二维傅里叶变换,只不过由于柱对称、二维变换只留下了对径向
r
∈
[
0
,
+
∞
)
r\in[0, +\infty)
r∈[0,+∞)的积分,变换的积分核也由
exp
(
i
k
x
x
+
i
k
y
y
)
\exp(ik_xx+ik_yy)
exp(ikxx+ikyy)退化为贝塞尔函数。
该变换常见于柱对称体系下,傍轴光束的传播方程(衍射积分方程),对于该体系的数值计算至关重要;另外也见于一些线性响应体系啥的。
数值汉克尔变换
数值汉克尔变换我见到了两种主流的算法,它们对径向 r r r轴的离散方式不同:
- r r r轴对数离散,即 r ≃ exp [ ( i − i 0 ) ∗ δ ] r\simeq \exp[(i-i_0)*\delta] r≃exp[(i−i0)∗δ],其中 i 0 , δ i_0,\delta i0,δ为常数, i i i为整数、从0取到较大的值。
- r r r轴近似线性离散。
python
的数值计算库scipy.fft.fht
提供了第一种方法的程序,可以在这里找到示例和引文。
第二种方式更适合于衍射方程的计算、均匀的线性的离散方式也更加符合数值计算的惯例。
接下来,我们参考相关文章,实现零阶汉克尔变换,并探讨数值精度。
引文是:
Li Yu, et.al, Quais-discrete Hankel transform, Optics Letters, 23, 409, 1998.
更高阶的汉克尔变换可以参考其后续工作,以及这篇文章:
Manuel Guizar-Sicairos and Julio C. Gutie´rrez-Vega, Computation of quasi-discrete Hankel transforms of integer order for propagating optical wave fields, J. Opt. Soc. Am. A, Vol 21, 53, 2004.
另外也有第三方库PyHank: Quasi-Discrete Hankel Transforms for Python作为一个参考。
零阶数值汉克尔变换:代码实现
根据引文Li Yu, et.al, Quais-discrete Hankel transform, Optics Letters, 23, 409, 1998.
首先构建一个类,初始化的时候构建内部的基本数据,比如
r
,
ρ
r, \rho
r,ρ轴等:
"""
an implementation of quasi-discrete Hankel transform.
Reference: Li Yu, et.al, Quais-discrete Hankel transform, Optics Letters, 23, 409, 1998.
@author: XXX
@date: 2021/12/19
"""
import numpy as np
from scipy.special import jn_zeros, j0, j1, jv
class Hankel_qDHT:
"""
an implementation of quasi-discrete Hankel transform.
R1 = R2 == R is assumed.
Reference: Li Yu, et.al, Quais-discrete Hankel transform, Optics Letters, 23, 409, 1998.
"""
def __init__(self, N=1000):
"""
N: total nodes used, of Bessel 0-order function.
generate:
.r: 1D-array, the corresponding radius axis data.
._S: _S = j_(N), the (N+1)-th zero point of J_0(x) function.
.js: the N+1 positive zeros of J_0(x) function.
.j1_inv: the function value of 1/(J_1(x)) for x at each .js.
.R: = (._S/2/pi)**0.5, then
.r = [j_1, j_2, j_3, ..., j_N ] /(2pi*R)
.C: 2D matrix, symmetric, C_ij = 2/S*J_0(j_i*j_j/S) / |J_1(j_i)| / |J_1(J_j)|
.F_multiplier" 1D array, F_multiplier_i = R/|J_1(j_i)|
"""
self.N = N
self.js = jn_zeros(0, self.N + 1)
self._S = self.js[-1] # j_(N+1)
self.R = (self._S/2.0/np.pi)**0.5 # r's cutoff position.
self.r = self.js[0:-1] / (2.0*np.pi*self.R) # the r axis for the field.
self.j1_inv = 1.0 / np.abs(j1(self.js))
self.F_multiplier = self.j1_inv[0:-1]*self.R
self.F_multiplier_inv = 1.0 / self.F_multiplier
# 1 / |J_1(n)*J_1(m)|:
J1_inv_mesh_x, J1_inv_mesh_y = np.meshgrid(self.j1_inv[0:-1]*(2.0/self._S), self.j1_inv[0:-1])
#self.Cmn = J1_inv_mesh_x * J1_inv_mesh_y * j0(np.outer(self.js[0:-1], self.js[0:-1]/self._S) ) #
self.Cmn = np.outer(self.j1_inv[0:-1]*(2.0/self._S), self.j1_inv[0:-1]) * j0(np.outer(self.js[0:-1], self.js[0:-1]/self._S) ) #
def transform(self, f1):
"""
perform 0-order Hankel transform. f2(r2) = 2pi*\int f1(r1)*J_0(2pi*r1*r2)*r1*dr1
return f2
Note: f1 must defined in self.r axis.
len(f1)==len(self.r)==self.N, and each r-point corresponding to J_0(x)'s zeros jn by:
rn = jn/(2*pi*R), wher 2*pi*R^2==j_(N+1)~(N+1)*pi
==> R ~ \sqrt((N+1)/2)
f1: 1D-array, float or complex valued. length == self.N, defined on self.r axis.
"""
if len(f1) != self.N:
print("invalid f1")
return None
f2 = self.F_multiplier_inv * np.matmul(self.Cmn, self.F_multiplier*f1)
# transform,
return f2
这里的Hankel_qDHT.__init__(self, N=1000)
函数内部构建了变换所需的数值环境,具体的定义可以参考文献;这里的N越大,所涵盖的径向范围
r
m
a
x
r_{max}
rmax就越大。
而Hankel_qDHT.transform(self, f1)
就对
f
1
(
r
1
)
f_1(r_1)
f1(r1)函数作数值的汉克尔变换,返回结果
f
2
(
r
2
)
f_2(r_2)
f2(r2)。需要特别注意的是,
r
1
,
r
2
r_1,r_2
r1,r2必须是__init__
中生成的给定值:Hankel_qDHT.r
。
这个类使用的时候,首先实例化一个对象,然后用这个对象来对不同的
f
1
(
r
)
f_1(r)
f1(r)函数作变换即可。
精度和效率
为了体现这种方法的精度,我们对比几种传统的数值积分方式。
首先,定义函数的定时器、数值积分的函数,作为对比:
# test and compare:
if __name__ == "__main__":
import time
import matplotlib.pyplot as plt
from scipy.integrate import trapezoid, simpson, romb
# ----------------------------------------------------------------------- #
# 一些辅助函数:计时器、作为对比的数值积分函数等等
def timming(func, total=100, name=""):
start = time.time()
for _ in range(total):
func()
end = time.time()
print(name, "total time [sec]: ", end-start)
#
def part_dot(cMN, cN, dtype=np.complex_):
"""
return an array shaped (M, N), with elements,
cij = cMN(i,j)*cN(j)
cMN: 2d array, shape==(M, N)
cN: 1d array, shape==(N,)
"""
# assert len(cN) >= cMN.shape[1]
res = np.zeros(shape=cMN.shape, dtype=dtype)
for i in range(cMN.shape[0]):
res[i, :] = cMN[i, :] * cN[:]
return res
def dot_integrate(cMN, cN, integrator):
"""
like matmul(cMN, cN), but using an integrator rather than sum.
cMN: 2d array, shape==(M, N)
cN: 1d array, shape==(N,)
integrator: callable, support complex-valued integration
integrator(y, x=None, dx=1.0, axis=-1)
supports: scipy.integrate.trapezoid, simpson and romb
return: 1d array, shape==(M, ),
values close to matmul(cMN, cN)
"""
part_cMNcN = part_dot(cMN, cN)
# res = np.zeros(shape=(cMN.shape[0], ))
# for i in range(len(res)):
# res[i] = integrator(part_cMNcN[i, :], dx=1.0)
# return res
return integrator(part_cMNcN, dx=1.0, axis=1)
# ----------------------------------------------------------------------- #
函数timming
对传入的函数调用N次并计时,用来评价效率;函数part_dot
和dot_integrate
一起用来计算数值积分,可以结合后面的调用来理解。
同样在if __name__ == "__main__":
区块下,构建函数
f
1
(
r
)
=
exp
{
−
π
r
2
1
+
0.2
i
}
f_1(r)=\exp\{-\frac{\pi r^2}{1+0.2i}\}
f1(r)=exp{−1+0.2iπr2},它的零阶汉克尔变换可以由理论给出:
f
2
(
r
)
=
(
1
+
0.2
i
)
exp
{
−
π
(
1
+
0.2
i
)
r
2
}
f_2(r)=(1+0.2i)\exp\{-\pi (1+0.2i)r^2\}
f2(r)=(1+0.2i)exp{−π(1+0.2i)r2}。
N = 1024
# qDHT:
ht = Hankel_qDHT(N) # <--- HERE
r = ht.r # <--- HERE
# sum or numerical-integration as integral:
# length of rho = 4096 + 1, ~(2^k+1),
# so that trapezoid, simpson, romb are supported at the same time.
rho = np.linspace(0., ht.R, N+1) # start from 0.0 or ht.r[0], doesn't matter much.
drho = (ht.R - 0.0) / (N+1 - 1.0)
c1 = 2.0*np.pi*drho * rho # 2pi*r*dr
c2 = j0(np.outer(2.*np.pi*rho, rho) ) # J_0(r*r')
# field:
a, c = 1.0, 0.2
field = lambda r: np.exp(-np.pi*r**2/(a+1.0J*c))
ft_theory = lambda r: (a+1.0J*c)*np.exp(-np.pi*(a+1.0J*c)*r**2)
# ----------------------------------------------------------------------- #
接下来,用4种积分的方式对比Hankel_qDHT
的结果:
- 直接求和,
sum
代替积分 - 使用
scipy.integrate.trapezoid
- 使用
scipy.integrate.simpson
- 使用
scipy.integrate.romb
并将它们和理论结果theory
画到一起:
# numerical integration results:
# 1.0 by numerical Hankel transform:
f_qDHT = ht.transform(field(ht.r))
# 2.0 by direct sum:
f_sum = np.matmul(c2, c1*field(rho))
# 3.0 by trapezoid integration:
f_trapezoid = dot_integrate(c2, c1*field(rho), trapezoid)
# 4.0 by simpson integration:
f_simpson = dot_integrate(c2, c1*field(rho), simpson)
# 5.0 by romb integration:
f_romb = dot_integrate(c2, c1*field(rho), romb)
# display:
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(ht.r, np.abs(f_qDHT), "-", c='b', label="qDHT", lw=1.5)
ax.plot(rho, np.abs(f_sum), "-", c='g',label="sum", lw=1.5)
ax.plot(rho, np.abs(f_trapezoid), "-", c="gray", label="trapezoid", lw=1.5)
ax.plot(rho, np.abs(f_simpson), "-", c="pink", label="simpson", lw=1.5)
ax.plot(rho, np.abs(f_romb), "-", c="yellow", label="romb", lw=1.5)
ax.plot(rho, np.abs(ft_theory(rho)), "r--", label="theory", lw=1)
ax.set_yscale("log")
# ax.set_xscale("log")
plt.title("Hankel_0{$e^{\\frac{-\pi r^2}{%.2f+%.2fJ}}$ }"%(a, c))
plt.legend()
plt.show()
# r points are amost linearly space:
# fig = plt.figure()
# ax = fig.add_subplot(111)
# ax.plot(np.linspace(0, ht.r.size, ht.r.size), ht.r, "k*-")
# ax.set_xlabel("index")
# ax.set_ylabel("r points of qDHT")
# plt.show()
# ----------------------------------------------------------------------- #
输出的
∣
f
2
(
r
)
∣
|f_2(r)|
∣f2(r)∣图像是:
可以看到,精度上qDHT > romb > simpson > trapezoid = sum
。我们的Hankel_qDHT
表现最好,数值误差在
1
0
−
17
10^{-17}
10−17量级,是float64
的极限。
另一方面,200次计算耗费时间对比:
print(np.amax(np.matmul(ht.Cmn, ht.Cmn)-np.eye(ht.N)))
# Timming:
t1 = lambda : ht.transform(field(ht.r))
t2 = lambda : np.matmul(c2, c1*field(rho))
Nt = 200
timming(t1, Nt, name="qDHT")
timming(t2, Nt, name="sum")
timming(lambda x=0:dot_integrate(c2, c1*field(rho), trapezoid),
Nt, name="trapezoid")
timming(lambda x=0:dot_integrate(c2, c1*field(rho), simpson),
Nt, name="simpson")
timming(lambda x=0:dot_integrate(c2, c1*field(rho), romb),
Nt, name="romb")
# ----------------------------------------------------------------------- #
可能的输出如下:
3.2132629890213593e-13
qDHT total time [sec]: 1.061779499053955
sum total time [sec]: 0.9278731346130371
trapezoid total time [sec]: 3.752824068069458
simpson total time [sec]: 3.1270384788513184
romb total time [sec]: 1.8081090450286865
耗时上Hankel_qDHT
与sum
接近,远快于trapezoid, simpson, romb
等。
当然,Hankel_qDHT
需要初始化参数、该过程会需要一点时间(<1sec);好在全局只需要初始化一个实例即可。
总结
实现了零阶汉克尔变换的数值版本。既快速、又精确。