# Function to Represent Recursive Integral

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.61350081547709315.954589770191005*^27949 and 3.61350081547709315.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?

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.

• 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! – Joel Kulesza Mar 1 '15 at 19:23

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.