Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I am interested in implementing Hahn-Exton's q-Bessel Function of the first and second kinds, $J_\alpha^{(3)}(x;q)$ and $Y_\alpha^{(3)}(x;q)$ , in Mathematica.

There's no difficulty in defining the Hahn-Exton q-Bessel of the first kind, exactly as defined in the Wikipedia page, by using Mathematica's QHypergeometricPFQ:

QBesselJ[n_, q_, z_] := 
 QPochhammer[q^(n + 1), q]/
   QPochhammer[q]*(z)^n QHypergeometricPFQ[{0}, {q^(n + 1)}, 
   q, q*z^2]

Where $n$ is the order and $q$ is the basic parameter. However, defining a q-analogue of the Bessel function of the second kind is more difficult, because of the limit needed for integer index, $Y_n^{(3)}(x;q)=\lim_{\alpha \rightarrow n}Y_\alpha^{(3)}(x;q)$. I am currently using the following cumbersome expression:

QBesselY[n_, q_, z_] := 
 If[Element[n, Integers], 
  2*(q - 1)/(\[Pi]*Log[q])*Log[z/(1 - q)]*QBesselJ[n, q, z] + 
   S1[n, q, z] + S2[n, q, z] + S3[n, q, z], 
  QGamma[n, q] QGamma[1 - n, q]/(\[Pi])*
   q^(-n^2/2) (Cos[\[Pi]*n]*q^(n/2)*
      QBesselJ[n, q, z] - QBesselJ[-n, q, z*q^(-n/2)])]

Where S1[n, q, z], S2[n, q, z] and S3[n, q, z] are complicated sums, given by:

S1[n_, q_, z_] := 
 If[n == 0, 
  0, -(1 - q)/(\[Pi])*z^(-n)*
   Sum[QPochhammer[q, q, n - k - 1]*z^(2 k)/QPochhammer[q, q, k], {k, 
     0, n - 1}]]

S2[n_, q_, z_] := (1 - q)/(\[Pi]*Log[q])*z^n*
  Sum[(-1)^k*q^(k (k + 1)/2)*
    z^(2 k)/(QPochhammer[q, q, k]*
       QPochhammer[q, q, n + k])*(QPolyGamma[0, 1 + k + n, q] + 
      QPolyGamma[0, 1 + k, q]) , {k, 0, Infinity}]

S3[n_, q_, z_] := -(1 - q)/(\[Pi]*2)*z^n*
  Sum[(-1)^k*q^(k (k + 1)/2)*
    z^(2 k)*(2 k + 1)/(QPochhammer[q, q, k]*
       QPochhammer[q, q, n + k]), {k, 0, Infinity}]

As it stands, these sums can only be evaluated numerically and take a long calculating time. I would like to find an alternative implementation that avoids this issue. In the case of the ordinary Bessel of the second kind, I believe Mathematica avoids calculating the limit for integer index by using a representation in terms of Meijer G functions, however, implementing a basic analogue of these seems to me a hopeless task.

So the question is, what methods should I use in implementing these special functions, in a way that allows relatively quick evaluation?

share|improve this question
Welcome here. Would you consider to reformat your question? Please mark code-blocks and inline code appropriately as described here. –  halirutan Aug 14 '13 at 2:00
add comment

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Browse other questions tagged or ask your own question.