2014-07-09 7 views
8

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:

sigma,

Dm

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:

enter image description here

Ecco cosa la trama dovrebbe simile (guardare il m = 11 della curva a destra):

enter image description here

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

enter image description here

Fro m questo rapporto si può anche trovare la differenza della prima derivata della funzione Hankel attraverso il taglio:

enter image description here

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.

+2

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. –

+1

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/) –

+0

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

risposta

1

Il taglio di ramo di H1 in scipy è (-inf, 0) e non in (-1j * inf, 0) come previsto nella carta citata, il che spiega perché si ottengono risultati errati.

Il problema può essere risolto con un uso creativo di uno argument transformation for Y_nu, come indicato sopra nei commenti.

Supponiamo di assumere l'ordine intero n. Abbiamo

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

jv non ha diramazioni (in ordine intero), ma yv ha. La formula trasformazione legge

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

o, in altre parole,

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) è definito in SciPy avere ramo tagliato a (-inf, 0), e iv (n , z) non ha tagli di diramazione (in ordine intero). A parte i logaritmi, i diramazioni sul RHS sono quindi in (-1j * inf, 0), esattamente dove vogliamo che siano. L'unica cosa che rimane da fare è scegliere il taglio del ramo del termine logaritmo in modo appropriato.

Il prolungamento analitico corretto con taglio ramo in (-1j * inf, 0) è quindi

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) 

Questo coincide esattamente con yv(n, z) in 3 dei 4 quadranti. Ha un ramo tagliato in (-1j * inf, 0). Inoltre, è ovviamente una funzione analitica. Pertanto, è lo stesso di yv, ma con una scelta di taglio di ramo diversa.

Abbiamo poi

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

È ovviamente la funzione Hankel con taglio ramo in (-1j * inf, 0).

Sulla base di questo, è anche possibile elaborare le derivate.

+0

Grazie per la tua risposta, aiuta molto. Ho pensato che questo avrebbe risolto il tutto, ma sto ancora avendo problemi a tracciare la funzione sigma. Le mie scuse per aver rimosso lo stato "risposto" dal post. – PeteyCoco