编程功能包含切入负虚轴

Question

最终更新

我通过电子邮件发送了论文的作者,结果发现sigma的等式中有一个错误.我给了pv最好的答案,因为他们确实帮助回答了所说的问题.

第一次尝试 我试图编写下面函数的数值表示:

,

并且'+'/' - '上标表示当z接近分支切割时的限制,其沿着负假想的半轴.H和J是Hankel和Bessel函数.其余变量(n_r,m,R)取决于问题的几何形状.我希望沿着相对于k的负假想半轴绘制该函数.我当前的代码(增加了pv)如下)

import scipy as sp
import numpy as np
import matplotlib.pyplot as plt
from numpy import pi
from scipy.special import jv, iv, kv, jvp, ivp, kvp

m = 11  # Mode number
nr = 2  # Index of refraction
R = 1   # Radius of cylinder
eps = 10e-8


def yv_imcut(n, z):
    return -2/(pi*(1j)**n)*kv(n, -1j*z) + 2*(1j)**n/pi * (0.5j*pi) * iv(n,-1j*z)

def yvp_imcut(n, z):
    return (n/z)*yv_imcut(n,z) - yv_imcut(n+1,z)

def hankel1_imcut(n, z):
    return jv(n, z) + 1j*yv_imcut(n, z)

def h1vp_imcut(n, z):
    return jvp(n, z, 1) + 1j*yvp_imcut(n, z)

# Define the characteristic equation
def Dm(n, z):
    return nr*jvp(n, nr*z, 1) * hankel1_imcut(n, z) - jv(n, nr*z)*h1vp_imcut(n,z)

# Define the cut pole density function
def sigma(k,n):
    return  4*(nr**2 - 1)*jv(n,nr*k*R)/(pi**2 * k * ((Dm(n, k*R-eps).real)**2 + (Dm(n, k*R+eps).imag)**2))

k = np.linspace(-eps*1j, -15j,1000)
y = sigma(k,m)
x = np.linspace(0,15,1000)

plt.plot(x, y.imag)
plt.show()

这是我在负虚轴上的sigma.imag图:

这是情节应该是什么样的(看看右边的m = 11曲线):

用户光伏帮助我将Hankel函数的切割移动到负假想的半轴,但我的sigma图仍然是关闭的.我在论文中指出,它指出西格玛是"纯粹想象的"(第五页,第一栏)

这些方程式和图表来自本文第4页:http://arxiv.org/pdf/1302.0245v1.pdf

第二次尝试

文章的附录B说明了切割中Hankel函数的差异

从这个关系中,我们也可以找到切割中Hankel函数的一阶导数的差异:

我用这些公式编写了一个脚本:

def hankel1_minus(n,z):
    return hankel1(n,z) - 4*jv(n,z)

def h1vp_minus(n,z):
    return (n/z)*hankel1_minus(n,z) - hankel1_minus(n+1,z)

def Dm_plus(n, z):
    return nr *jvp(n, nr*z, 1) * hankel1(n, z) - jv(n, nr*z)*h1vp(n,z)

def Dm_minus(n, z):
    return nr *jvp(n, nr*z, 1) * hankel1_minus(n,z) - jv(n, nr*z)*h1vp_minus(n,z)

def sigma(k,n):
    return  4*(nr**2 - 1)*jv(n,nr*k*R)/(pi**2 * k * (Dm_plus(n, k*R) *     Dm_minus(n,k*R)).real)

绘制此sigma会得到与第一种方法相同的结果.

Answer 1

scipy 中 H1 的分支切割是 (-inf, 0) 而不是您引用的论文中预期的 (-1j*inf, 0)，这解释了为什么您得到不正确的结果。

正如我在上面的评论中指出的那样，可以通过对 Y_nu 的参数转换进行一些创造性的使用来解决这个问题。

让我们假设整数阶 n。我们有

hankel1(n, z) = jv(n, z) + 1j*yv(n, z)

jv 没有分支切割（整数顺序），但 yv 有。变换公式为

yv(n, 1j*z) = -2/(pi*(1j)**n)*kv(n, z) + 2*(1j)**n/pi * (log(1j*z) - log(z))*iv(n,z)

或者，换句话说，

yv(n, z) = -2/(pi*(1j)**n)*kv(n, -1j*z) + 2*(1j)**n/pi * (log(z) - log(-1j*z))*iv(n,-1j*z)

kv(n,z) 在 Scipy 中定义为在 (-inf, 0) 处有分支切割，而 iv(n,z) 没有分支切割（整数顺序）。因此，除了对数之外，RHS 上的分支切割位于 (-1j*inf,0 ) 中，恰好位于我们想要的位置。剩下要做的唯一的事情就是适当地选择对数项的分支切割。

在 (-1j*inf, 0) 中进行分支切入的正确分析延续是

def yv_imcut(n, z):
    return -2/(pi*(1j)**n)*kv(n, -1j*z) + 2*(1j)**n/pi * (0.5j*pi) * iv(n,-1j*z)

yv(n, z)这与4 个象限中的 3 个完全一致。它有一个切入 (-1j*inf,0) 的分支。而且，它显然是一个解析函数。因此，它与 相同yv，但具有不同的分支切割选择。

然后我们有

def hankel1_imcut(n, z):
    return jv(n, z) + 1j*yv_imcut(n, z)

显然这是在 (-1j*inf, 0) 处进行分支切入的 Hankel 函数。

在此基础上，还可以算出导数。