6
$\begingroup$

I am looking for the approximation of the following function:

$$\rho(a,b)=1-e^{-(a+b)}\sum_{m=1}^{\infty}\left(\sqrt{\frac{b}{a}}\right)^m I_m(2\sqrt{ab})$$ where $I_m(x)$ is the modified Bessel function. Since $I_m(x)$ is again an infinite series, up to how many terms I need to do summation in both the cases above in order to get some better approximation? Is there any rule of thumb?

  • 0
    Actually, you can use the trapezoidal rule to compute the modified Bessel function of the first kind. That seems to me simpler than using infinite series. But don't do forward recursion over the Bessel function, as that is unstable.2011-09-09
  • 2
    Bessel function $I_m(2 z) \sim \frac{z^m}{\Gamma(m+1)}$ for large $m \gg \vert z \vert$. This should give a pessimistic bound on how many terms you need to sum.2011-09-09
  • 0
    Here's an idea: for large $m$, $I_m(z)$ is well approximated by $\frac1{\sqrt{2\pi m}}\left(\frac{e z}{2m}\right)^m$. Find the value of $m$ that makes this approximation tiny (perhaps via Newton-Raphson), compute the modified Bessel function there, and recurse back. On the other hand, you might also want to take into account the other stuff multiplying the modified Bessel function when you're estimating $m$.2011-09-09
  • 0
    Thanks for your comments. But up to what value of $m$ shall I perform summation?2011-09-09
  • 0
    @shaikh: How explicitly accurate of any estimate do you want?2011-09-09
  • 0
    any estimate up to 5-6 decimal places would be sufficient2011-09-09
  • 1
    @shaikh, what typical $(a,b)$-regime are you interested in?2011-09-10

1 Answers 1

5

If one is interested in expressing $\rho(a,b)$ with a single summation/integration operation, one can use, for every positive $a$ and $b$ such that $b

Case $b< a:$ This uses only the integral formula for modified Bessel functions with integer parameter, which says that for every integer $m$ and every $x$ with nonnegative real part, $$ I_m(x)=\int\limits_0^\pi e^{x\cos\theta}\cos(m\theta)\frac{\mathrm{d}\theta}\pi. $$ Writing the cosine as the real part of a complex exponential and summing this over $m\ge1$ yields, for every real number $s$ such that $|s|<1$, $$ \sum\limits_{m\ge1}s^mI_m(x)=\text{Re}\int\limits_0^\pi e^{x\cos\theta}\frac{se^{i\theta}}{1-se^{i\theta}}\frac{\mathrm{d}\theta}\pi. $$ Multiplying the numerator and the denumerator by the conjugate quantity $1-se^{-i\theta}$ and using the relation $|1-se^{i\theta}|^2=1+s^2-2s\cos\theta$, one sees that $$ \sum\limits_{m\ge1}s^mI_m(x)=\int\limits_0^\pi e^{x\cos\theta}\frac{s\cos\theta-s^2}{1+s^2-2s\cos\theta}\frac{\mathrm{d}\theta}\pi. $$ Plugging in $x=2\sqrt{ab}$ and $s=\sqrt{b/a}$ (hence the restriction $b< a$) yields $\rho(a,b)$ as above.

Case $a=b$: One sees that $$ \rho(a,a)=1-e^{-2a}\sum\limits_{m\ge1}I_m(2a). $$ Using the fact that $I_m(x)=I_{-m}(x)$ for every integer $m$ and the formula (written somewhere in the WP page already mentioned) for the generating function $$ \sum\limits_{m=-\infty}^{+\infty}I_m(x)\cos(m\theta)=e^{x\cos\theta}, $$ at $\theta=0$, one sees that $$ \rho(a,a)=\frac12+\frac12e^{-2a}I_0(2a)=\frac12+\frac12e^{-2a}\int\limits_0^\pi e^{2a\cos\theta}\frac{\mathrm{d}\theta}\pi. $$ Case $a Comparing $\rho(a,b)+\rho(b,a)$ with the generating function of the family $(I_m(x))$ that we recalled above, one gets $$ \rho(a,b)+\rho(b,a)=1+e^{-(a+b)}I_0(2\sqrt{ab}), $$ for every $a$ and $b$. Hence, when $a< b$, $$ \rho(a,b)=e^{-(a+b)}\int\limits_0^\pi e^{2\sqrt{ab}\cos\theta}\frac{b-\sqrt{ab}\cos\theta}{a+b-2\sqrt{ab}\cos\theta}\frac{\mathrm{d}\theta}\pi. $$ Remark: The result is discontinuous around $a=b$ since, for every fixed $a$, $\rho(a,b)\to\rho(a,a)+\frac12$ when $b\to a^-$ and $\rho(a,b)\to\rho(a,a)-\frac12$ when $b\to a^+$.

  • 0
    I'd like to see how this formula was derived.2011-09-10
  • 0
    This seems to work nicely if $a > b$.2011-09-10
  • 2
    @Matt, here you are. (But a "please" would have been nice.)2011-09-10
  • 0
    @J.M. Why the restriction?2011-09-10
  • 0
    When I tried it out, the results of the integral and the sum came out different if $a \leq b$. Here's a *Mathematica* example: `With[{a = 1/4, b = 1/4}, (Sqrt[b]/Pi)*NIntegrate[Exp[2*Sqrt[a*b]*Cos[t]]*((Sqrt[a]*Cos[t] - Sqrt[b])/(a - 2*Sqrt[a*b]*Cos[t] + b)), {t, 0, Pi}, WorkingPrecision -> 30]]` versus `With[{a = 1/4, b = 1/4}, NSum[Sqrt[b/a]^m*BesselI[m, 2*Sqrt[a*b]], {m, 1, Infinity}, WorkingPrecision -> 30]]` . I'm still trying to figure out why this is so.2011-09-10
  • 0
    @J.M. See my Edit for the case $a=b$ (which your test with *Mathematica* is based on, if I am not mistaken). As they say, *the plot thickens*... :-))2011-09-10
  • 0
    Didier, your derivation looks to be fine to me. :) That's why I'm scratching my head about the different behavior when $a \leq b$. I'm currently running more numerical tests to look into this and get back to you (unless someone turns up an explanation for this strangeness first).2011-09-10
  • 1
    @J. M. ...Except that, to sum the series, I assumed that $|s|<1$, that is, that $b2011-09-10
  • 0
    Okay, I thought I was going nuts. I was working on an entirely different algorithm, and yours seemed a bit simpler. Anyway, he now has a simple method for the case $a > b$. :)2011-09-10
  • 0
    For this version: great analysis! :D2011-09-10
  • 2
    @Didier Piau: Thank you so much! I tried not to put pressure on you. I'm sorry I didn't say please. Please forgive me! :-)2011-09-10
  • 0
    @Matt, no problem.2011-09-10
  • 0
    Thanks a lot to all of you for your time and efforts to derive the above formulas. Can you suggest which approximating rule i.e. The Trapezoidal Rule or Simpson Rule etc. would be best to approximate the above integral? or is it possible to get the polynomial of the above integral which can be used in computer programs to approximate the integrals in the answer? My question may seem stupid but I am not good at Integration. Thanks a lot.2011-09-11
  • 0
    @shaikh, what typical (a,b)-regime are you interested in? (Bis repetita...)2011-09-11
  • 0
    @Didier: a and b can be any positive real numbers. i.e. $0 \leq a,b \leq \infty $2011-09-11
  • 0
    @shaikh: if $a$ and $b$ are reasonably sized, I would suppose plain Romberg quadrature should be fine (which is easily done as a modification of a program for implementing the trapezoidal rule).2011-09-12
  • 0
    Is it possible to write the case a=b, i.e., $\rho(a,a)=\frac12+\frac12e^{-2a}I_0(2a)$ without Bessel function or in terms of $Cos\theta$ as in the case of ab?2011-09-13
  • 1
    Yes, it is. See alternative formula for $\rho(a,a)$ in the revised version.2011-09-13
  • 0
    @Didier: Thanks a lot.2011-09-13
  • 0
    @Didier: In case a < b, is it $\rho(a,b)=e^{-(a+b)}\int\limits_0^\pi e^{2\sqrt{ab}\cos\theta}\frac{a-\sqrt{ab}\cos\theta}{a+b-2\sqrt{ab}\cos\theta}\frac{\mathrm{d}\theta}\pi$ or $\rho(a,b)=e^{-(a+b)}\int\limits_0^\pi e^{2\sqrt{ab}\cos\theta}\frac{b-\sqrt{ab}\cos\theta}{a+b-2\sqrt{ab}\cos\theta}\frac{\mathrm{d}\theta}\pi$. The first expression is not giving me right result but the second on seems fine???2011-09-13
  • 0
    Right. $ $ $ $ $ $2011-09-13