【问题标题】:Programming function containing cut in negative imaginary axis包含负虚轴切割的编程功能
【发布时间】: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) 上定义函数。我没有使用新切割在新的黎曼曲面上定义函数,而是将现有曲面旋转到另一个位置。
  • 请参阅下面的答案,了解如何使用转换来移动分支切割。

标签: python numpy scipy sympy


【解决方案1】:

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

正如我在上面的 cmets 中指出的那样,可以通过创造性地使用 argument transformation for 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)

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

然后我们有

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

显然是Hankel函数,分支切入(-1j*inf, 0)。

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

【讨论】:

  • 感谢您的回答,对您有很大帮助。我认为这可以解决整个问题,但我仍然无法绘制 sigma 函数。我很抱歉从帖子中删除“已回答”状态。
猜你喜欢
  • 1970-01-01
  • 1970-01-01
  • 2011-03-29
  • 2014-02-04
  • 2013-09-23
  • 1970-01-01
  • 2022-08-16
  • 2020-12-04
  • 1970-01-01
相关资源
最近更新 更多