aggiornamento finalefunzione Programmazione contenenti tagliare in asse immaginario negativo
Ho contattato l'autore della carta e come risulta c'era un errore nell'equazione per sigma. Ho dato la migliore risposta al PV perché hanno aiutato a rispondere al problema come affermato.
primo tentativo sto cercando di programmare una rappresentazione numerica della funzione di seguito:
,
e il '+'/'-' apici indicano i limiti mentre z si avvicina al taglio del ramo, che si trova lungo il mezzo asse immaginario negativo. Le H e J sono funzioni di Hankel e Bessel. Il resto delle variabili (n_r, m, R) dipende dalla geometria del problema. Vorrei tracciare questa funzione lungo il semiasse immaginario negativo rispetto a k. Il mio codice corrente (con l'aggiunta del pv) è la seguente)
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()
Ecco il mio appezzamento di sigma.imag lungo l'asse immaginario negativo:
Ecco cosa la trama dovrebbe simile (guardare il m = 11 della curva a destra):
L'utente pv mi ha aiutato a spostare il taglio della funzione di Hankel sul semiasse immaginario negativo, ma la mia trama di sigma è ancora disattivata. Ho osservato nella carta che afferma che sigma è "puramente immaginario" (in alto cinque, prima colonna)
Queste equazioni e il grafico provengono da pagina 4 in questo articolo: http://arxiv.org/pdf/1302.0245v1.pdf
secondo tentativo
Appendice B del articolo afferma la differenza della funzione di Hankel attraverso il taglio come
Fro m questo rapporto si può anche trovare la differenza della prima derivata della funzione Hankel attraverso il taglio:
ho scritto uno script con queste formule:
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)
Tracciando questo sigma dà la stesso risultato del primo metodo.
La formula, come scritto, produce la curva che si ottiene. Il problema sembra essere che la carta usa una scelta non standard per il taglio di ramo della funzione di Hankel. Scipy (e Mathematica) posizionano entrambi il taglio lungo l'asse * reale * negativo piuttosto che l'asse immaginario negativo assunto nel foglio. –
Probabilmente è possibile aggirare il problema del ramo con un uso creativo della [formula per Y_m (-iz)] (http://functions.wolfram.com/Bessel-TypeFunctions/BesselY/16/01/01/) –
Grazie per il promemoria sul taglio del ramo. Ho aggiornato la domanda con queste informazioni. Darò al problema un altro tentativo con questo in mente. – PeteyCoco