【发布时间】:2014-08-31 02:14:12
【问题描述】:
最终更新
我给论文的作者发了电子邮件,结果发现 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 曲线):
用户 pv 帮助我将 Hankel 函数的切割移动到负虚半轴,但我的 sigma 绘图仍然关闭。我在论文中指出,sigma 是“纯虚构的”(第五页顶部,第一列)
这些方程式和图表来自本文第 4 页:http://arxiv.org/pdf/1302.0245v1.pdf
第二次尝试
文章的附录B将汉克尔函数的差异描述为:
从这个关系中,我们还可以找到 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 得到与第一种方法相同的结果。
【问题讨论】:
-
公式,如所写,产生你得到的曲线。问题似乎是该论文对 Hankel 函数的分支切割使用了非标准选择。 Scipy(和 Mathematica)都沿着负 real 轴而不是论文中假设的负虚轴放置切口。
-
您可以通过创造性地使用formula for Y_m(-i z)来解决分支切割问题
-
感谢关于分支切割的提醒。我已经用这些信息更新了这个问题。考虑到这一点,我会再试一次。
-
将 Neumann/Hankel 函数及其参数乘以 i 确实会移动切割,但该函数的值仍然与切割沿负实轴时的值相同。据我所知,SciPy 中的 Hankel 函数是从 ( -pi, pi ) 定义的,并且将切割移动到负虚轴意味着在 (-pi/2, 3pi/2) 上定义函数。我没有使用新切割在新的黎曼曲面上定义函数,而是将现有曲面旋转到另一个位置。
-
请参阅下面的答案,了解如何使用转换来移动分支切割。