Function to Represent Recursive Integral

Question

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?

Answer 1

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.

Answer 2

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}]

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