编程函数包含负虚轴的切割

Question

最终更新

我给论文的作者发了邮件，结果发现sigma的公式有个错误。我给pv的回答打了最高分，因为他们确实帮助我解决了问题。

第一次尝试

我正在尝试编写一个程序，用来数值表示下面这个函数：

$sigma$ ,

$Dm$

这里的'+'和'-'上标表示当z接近分支切割线时的极限，这条切割线位于负虚数半轴上。H和J分别是汉克尔函数和贝塞尔函数。其余的变量（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的图：

enter image description here

这个图应该是这样的（看右边的m = 11曲线）：

enter image description here

用户pv帮助我把汉克尔函数的切割线移到了负虚数半轴上，但我的sigma图还是不对。我在论文中注意到，sigma被描述为“纯虚数”（第五页顶部，第一列）。

这些公式和图来自这篇文章的第四页：http://arxiv.org/pdf/1302.0245v1.pdf

第二次尝试

文章的附录B说明了汉克尔函数在切割线两侧的差异：

enter image description here

根据这个关系，我们还可以找到汉克尔函数一阶导数在切割线两侧的差异：

enter image description here

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

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的结果与第一种方法是一样的。

数值计算数学建模极限贝塞尔函数函数绘图汉克尔函数复分析分支切割

编程函数包含负虚轴的切割

1 个回答

撰写回答