1
$\begingroup$

I'd like to represent the following recursive integral equation to evaluate/graph (for $n\leq3$):

$$K_{i,n}\left(x\right)=\int_{x}^{\infty}K_{i,n-1}\left(y\right)dy$$

where

$$K_{i,0}\left(x\right)=K_{0}\left(x\right)=\int_{0}^{\infty}\frac{\exp\left(-\sqrt{\tau^{2}+x^{2}}\right)}{\sqrt{\tau^{2}+x^{2}}}d\tau$$

The evaluation can be analytic or numeric to serve my purposes. I've tried this:

Kin[x_, 0] := NIntegrate[Exp[-Sqrt[τ^2 + x^2]]/Sqrt[τ^2 + x^2], {τ, 0, ∞}]
Kin[x_, n_] := NIntegrate[Kin[y, n-1], {y, x, ∞}]

which gives:

In[48]:= Kin[1.0, 0]
Out[48]= 0.421024

but also

In[49]:= Kin[1.0, 1]

During evaluation of In[49]:= NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small. >>
During evaluation of In[49]:= NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections in y near {y} = {8.16907*10^224}. NIntegrate obtained 3.613500815477093`15.954589770191005*^27949 and 3.613500815477093`15.954589770191005*^27949 for the integral and error estimates. >>

Out[49]= 3.613500815477093*10^27949

These are Bickley-Nayler functions and have definitely been evaluated numerically. The plots from those evaluations look sane and don't indicate that divergence should occur. What have I done wrong in my function definition?

$\endgroup$
2
$\begingroup$

I noticed one mistake: your second definition had $K_{n-1}(x)\,dy$ in the integrand, and not $K_{n-1}(y)\,dy$. Here is the corrected definition:

Kin[x_?NumericQ, 0] := 
 NIntegrate[Exp[-Sqrt[τ^2 + x^2]]/Sqrt[τ^2 + x^2], {τ, 0, ∞}]
Kin[x_?NumericQ, n_] := NIntegrate[Kin[y, n - 1], {y, x, ∞}]

Kin[1, 0] (* 0.421024 *)
Kin[1, 1] (* 0.328286 *)

Note the use of NumericQ to prevent NIntegrate from trying to operate on a symbolic quantity.

$\endgroup$
  • $\begingroup$ Outstanding, thanks for the quick reply! Indeed, I saw the typo and corrected it just before you posted your answer. The NumericQ is likely what I was missing all along as I was trying different things: it seemed like it was complaining about the bounds but I didn't know how to tell it to treat them appropriately. I'll definitely keep that ability in mind for the future! $\endgroup$ – Joel Kulesza Mar 1 '15 at 19:23
2
$\begingroup$

Here is another recursive implementation of the Bickley(-Nayler) function $\mathop{\mathrm{Ki}_n}(z)$, using Leonid's method from here:

SetAttributes[BickleyKi, Listable];
BickleyKi[n_Integer?NonPositive, z_] := (-1)^n Derivative[0, -n][BesselK][0, z]; 
BickleyKi[1, z_] := π (1 - z BesselK[{0, 1}, z].StruveL[{-1, 0}, z])/2;
BickleyKi[n_Integer, z_] := Module[{zl}, 
          Set @@ Hold[BickleyKi[n, zl_],
                      (zl (BickleyKi[n - 3, zl] - BickleyKi[n - 1, zl]) +
                       (n - 2) BickleyKi[n - 2, zl])/(n - 1)];
          BickleyKi[n, z]]

Examples:

N[BickleyKi[{0, 1}, 1]]
   {0.42102443824070834, 0.328286478171118}

With[{n = 7, z = 8/3},
     {N[BickleyKi[n, z], 20],
      NIntegrate[BesselK[0, t] (t - z)^(n - 1), {t, z, ∞},
                 WorkingPrecision -> 20]/Gamma[n]}]
   {0.028396677600849973457, 0.028396677600849973363}

The Bickley function $\mathop{\mathrm{Ki}_\alpha}(z)$ of arbitrary order $\alpha$ can in fact be expressed in terms of the Meijer $G$ function:

BickleyKi[α_, z_] := MeijerG[{{}, {(α + 1)/2}}, {{0, 1/2, α/2}, {}}, z/2, 1/2]/2

With this definition, one can now consider things like

BickleyKi[1/2, z] // FunctionExpand
   (π^(3/2) Sqrt[z] BesselI[-1/4, z/2]^2)/(2 Sqrt[2]) -
   (π^(3/2) Sqrt[z] BesselI[-1/4, z/2] BesselI[1/4, z/2])/Sqrt[2] +
   (π^(3/2) Sqrt[z] BesselI[1/4, z/2]^2)/(2 Sqrt[2])

or

Plot3D[Im[BickleyKi[1 - I, x + I y]], {x, -2, 2}, {y, -2, 2}]

plot of imaginary part of an arbitrary order Bickley function

Nevertheless, for integer order, the recursive implementation given above is still the most efficient evaluation method.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.