Scaled Modified Bessel Function

The Theory

presets...

Off[InverseFunction::ifun, Solve::ifun]; << 
      Graphics`Graphics`
thick = 0.01;
{Red, Green, 
          Blue} = {{Thickness[thick], RGBColor[1, 0, 0]}, {Thickness[thick], 
            RGBColor[0, 1, 0]}, {Thickness[thick], 
            RGBColor[0, 0, 1]}};
SetOptions[Plot, 
        TextStyle -> {FontFamily -> "Times", 
            FontSize -> 12}, ImageSize -> 300, 
        PlotStyle -> {Blue, Red, Green}];
SetOptions[LogPlot, 
        TextStyle -> {FontFamily -> "Times", 
            FontSize -> 12}, ImageSize -> 300, 
        PlotStyle -> {Blue, Red, Green}];

The scaled modified Bessel function (e^(-x )) (Subscript[I, 0])(
          x) of the first kind of order zero is calculated with the help of an  asymptotic  expansion with error bounds that can be found at [Olv74, p. 269, Ex. 13.2]. Specifically, we obtain for x > 0 that

\(TraditionalForm\`e\^\(-x\)\ \(\(I\_0\)(
            x)\) = \(1\/\@\(2  π\ x\)\) \((∑\+\(k = 0\)\%\(n - 1\)\((\(\((2\ k - 1)\)!!\))\)\^2\/\(\(k!\)\ \((8\ x)\)\^k\) + R\_n)\)\)

where, with χ(
          n) = (π^(1/2)) Γ((1/2) n + 
                1)/Γ((1/2) n + (1/2)),

\(TraditionalForm\`\(\(|\)\(R\_n\)\(|\)\(\(≤\)\(\ \)\(2 \( χ(n)\) e\^\(π/\((8  x)\)\)\ \(\((2  n - 1)\)!!\)\^2\/\(\(n!\)\ \((8  x)\)\^n\)\)\)\) = 
        2\ e\^\(π/\((8  x)\)\)\ \(Γ(n + 1\/2)\)\^2\/\(\(Γ(\(1\/2\) \((n + 1)\))\)\^2\ \((4  x)\)\^n\)\).

derive the last equality...

FullSimplify[
        2  (π^(1/2)) Γ((1/2) n + 
                1)/Γ((1/2) n + 
                (1/2)) (E^(π/(8  x))) ((((2  n - 1)!!)^2)/(n! (((8  x))^n))), {x > 0, n > 0, n ∈ Integers}] // 
      TraditionalForm

Stirling's formula yields the asymptotics

\(TraditionalForm\`\(Γ(n + 1\/2)\)\^2\/\(\(Γ(\(1\/2\) \((n + 1)\))\)\^2\ \((4  x)\)\^n\) ∼ \((\(2  n\)\/e)\)\^n\) ,

check this asymptotics...

Limit[(((Γ(n + (1/2)))^2)/(((Γ((1/2) (n + 
                        1)))^2) ))/(((((2  n)/E)))^n), n -> ∞]

in fact by monotonicity

\(TraditionalForm\`\(Γ(n + 1\/2)\)\^2\/\(\(Γ(\(1\/2\) \((n + 1)\))\)\^2\ \((4  x)\)\^n\) ≤ \ \(\(e\ π\)\/8\) \((\(2  n\)\/e)\)\^n\),	n ≥ 1.

check this claim...

(((Γ(n + (1/2)))^2)/(((Γ((1/2) (n + 
                        1)))^2) ))/(((((2  n)/E)))^n) /. n -> 1

Plot[((E π)/8) (((((2  n)/E)))^n)/(((Γ(n + (1/2)))^2)/(((Γ((1/2) (n + 1)))^2) )),{n, 1, 100},AxesLabel->{"n","ratio"}]

Summarizing, we have

\(TraditionalForm\`\(\(|\)\(R\_n\)\(|\)\(\(≤\)\(\(\(e\ π\)\/4\) e\^\(π/\((8  x)\)\)\ \((n\/\(2  e\ x\))\)\^n\)\(≤\)\(\ \)\(e\ \((n\/\(2  e\ x\))\)\^n\)\)\)\),	x ≥ 2.

check the last inequality...

((E π)/4) (E^(π/(8  x))) ≤  E /. x -> 2

Thus, we can guarantee the error \(TraditionalForm\`\(\(|\)\(R\_n\)\(|\)\(\(≤\)\(\ \)\(e\ ·10\^\(-d\)\)\)\)\), which essentially guarantees d correct decimal digits, if for a given x the summation index n is chosen in accordance to the solution \(TraditionalForm\`\(\(η\_*\)\(\ \)\)\)of the equation

\(TraditionalForm\`\(\(\((η\/\(2  e\ x\))\)\^η\)\(==\)\(10\^\(-d\)\)\(\ \)\)\),	that is 	\(TraditionalForm\`η\ \ \(log(\(2  e\ x\)\/η)\) == d\ log\ 10\( = \&def\)\ 2  x\_0\) .

Generally, if any, there will be two solution that a grouped around a minimum error. This behavior is typical for asymptotic expansions that do not converge as n -> ∞.

have a look at a typical situation...

LogPlot[((((η/(2  E x))))^η) /. 
        x -> 60. , {η, 3. , 300. }, 
      AxesLabel -> {"η", "error"}]

It turns out that the minimum is located at (Subscript[η, min]) = 2  x.

check this claim...

Solve[
      D[((((η/(2  E x))))^η), η] == 
        0, η]

We obtain

\(TraditionalForm\`\(η\_*\) = \(2\ x\_0\)\/\(-\(W(\(-x\_0\)/\((e\ x)\))\)\)\),

calculate the solution \(TraditionalForm\`\(η\_*\)\)

Solve[η  log(((2  E x)/η)) == 2  (Subscript[x, 0]), η]

where W denotes Lambert's W-function. There are exactly two branches of -W(-z) that are real for z ∈ (0, (e^(-1))]. We choose the smaller solution \(TraditionalForm\`\(η\_*\)\), that is, the branch for which -W(-z) ≥  1,  z\(TraditionalForm\`\(\(∈\)\((0, e\^\(-1\)]\)\)\). (The other has -W(-z) ≤ 
        1, z∈(0, (e^(-1))] .) In particular, -W(-(e^(-1))) = 1 and -W(-z)  ->  +∞ as z ↘ 0.This branch is given in Mathematica by -ProductLog[-1,-z]. Since -W(-z) is complex for z > (e^(-1)), we must restrict ourselves to \(TraditionalForm\`x\_0/\((e\ x)\) ≤ e\^\(-1\)\), that is x ≥ (Subscript[x, 0]).

have a look at the branches...

Plot[{-ProductLog[-1, -z], -ProductLog[0, -z], 
          1},{z, 0, (E^(-1))},AxesLabel->{"z","-W(-z)"}]

With that choice of the branch we have (Subscript[η, * ]) ≤  
          2  (Subscript[x, 0]) ≤  2  x = (Subscript[η, min]). Thus, the error is descending here and we choose the summation index n as the ceiling of \(TraditionalForm\`\(η\_*\)\),

\(TraditionalForm\`n = ⌈\(2\ x\_0\)\/\(-\(W(\(-x\_0\)/\((e\ x)\))\)\)⌉\),		x ≥ (Subscript[x, 0]) = (d/2) log 10.

The Code

BesselITilde[x_?Positive] := 
      Module[{x0 = Precision[x] Log[10]/2, n}, 
        If[x < x0, Exp[-x] BesselI[0, x], 
          n = Ceiling[
              2  x0/-ProductLog[-1, -x0/(E x)]]; (1/(Sqrt[2  π x])) 
            Sum[(( ((((2  k - 1)!!))^2))/(k! (((8  x))^k))), {k, 0, n - 1}]]]

a test...

BesselITilde[63.`50]

(E^(-x)) (Subscript[I, 0])(x) /. 
      x ->  63.`52

Sitzung S. 180 (erzeugt Tabelle 6.1)

ReturnProbability[d_Integer] := 
      1 - 1/NIntegrate[BesselITilde[t/d]^d, {t, 0, ∞}, 
            WorkingPrecision -> 40, PrecisionGoal ->  17, 
            Method ->  DoubleExponential]
{#, NumberForm[ReturnProbability[#], {15, 15}]} & /@ 
        Join[Range[3, 20], {30, 40, 100}] // TableForm