Python数值计算：数值汉克尔变换

笨笨笨的打工仔

已于 2024-02-28 13:22:20 修改

阅读量1.3k

点赞数 13

分类专栏： Python数值计算Tips 文章标签： python 笔记经验分享线性代数傅立叶分析 scipy numpy

于 2024-02-28 13:06:51 首次发布

本文链接：https://blog.csdn.net/qq_41197833/article/details/136342613

版权

Python数值计算Tips 专栏收录该内容

3 篇文章 0 订阅

订阅专栏

零阶汉克尔变换（Hankel Transform）

函数 $f_1(r)$ 的零阶汉克尔变换定义为：
$f_2(r_2)=2\pi\int_0^{+\infty}f_1(r_1)J_0(2\pi r_1r_2)r_1d r_1$
其中 $J_0(\cdot)$ 为零阶贝塞尔(Bessel)函数。其逆变换形式一样：
$f_1(r_1)=2\pi\int_0^{+\infty}f_2(r_2)J_0(2\pi r_1r_2)r_2d r_2$
本质上0阶汉克尔变换就是柱对称体系上的二维傅里叶变换，只不过由于柱对称、二维变换只留下了对径向 $r\in[0, +\infty)$ 的积分，变换的积分核也由 $exp(ik_xx+ik_yy)$ 退化为贝塞尔函数。
该变换常见于柱对称体系下，傍轴光束的传播方程（衍射积分方程），对于该体系的数值计算至关重要；另外也见于一些线性响应体系啥的。

数值汉克尔变换

数值汉克尔变换我见到了两种主流的算法，它们对径向 $r$ 轴的离散方式不同：

$r$ 轴对数离散，即 $r\simeq \exp[(i-i_0)*\delta]$ ，其中 $i_0,\delta$ 为常数， $i$ 为整数、从0取到较大的值。
$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.
首先构建一个类，初始化的时候构建内部的基本数据，比如 $\rho$ 轴等：

"""
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_{max}$ 就越大。
而Hankel_qDHT.transform(self, f1)就对 $f_1(r_1)$ 函数作数值的汉克尔变换，返回结果 $f_2(r_2)$ 。需要特别注意的是， $r_1,r_2$ 必须是__init__中生成的给定值：Hankel_qDHT.r。
这个类使用的时候，首先实例化一个对象，然后用这个对象来对不同的 $f_1(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\{-\frac{\pi r^2}{1+0.2i}\}$ ，它的零阶汉克尔变换可以由理论给出： $f_2(r)=(1+0.2i)\exp\{-\pi (1+0.2i)r^2\}$ 。

    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)|图像，与理论曲线对比
可以看到，精度上qDHT > romb > simpson > trapezoid = sum。我们的Hankel_qDHT表现最好，数值误差在 $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