0

$$\int_{r_{min}}^{\infty} \frac{dr}{r^2 \sqrt{1-\frac{V(r)}{E_c}-\frac{p^2}{r^2}}} $$ where $r_{min}$ is the root of the denominator.

$$V(r)=\frac{Z_1 Z_2 q^2}{4 \pi \epsilon_0 r} \Phi(r)$$

enter image description here

where $a_U$, $q$, $\epsilon_0$, $a_0$, $Z_1$, $Z_2$, $p$ are constant.

I tried some naive solutions, but the problem is that the thing to integrate (let's call it $Y$) approaches infinity when $r$ approaches $r_{min}$. So mathematically, this integral is supposed to converge (or is it?), but numerically it's ill-posed.

In the naive approach, I find $r_{min}$ numerically, but this gives a very high $Y$ at the beginning of the integration (that should be $\pm \infty$ if $r_{min}$ was exact), and then I try to integrate with very small trapezoids, but it seems this approach is fundamentally flawed.

So... I should probably transform this integral, changing variables and stuff, but I don't know where to begin...

Is there a known method to compute this ?

Yuriy S
  • 31,474

1 Answers1

0

I've googled "exp-sinh quadrature" thanks to user14717 comment.

I tried boost::math::quadrature::exp_sinh and this works flawlessly, thanks a lot!

  • 1
    I would use boost::math::quadrature::exp_sinh, rather than tanh-sinh. Tanh-sinh will work, but it's not optimal in this case. Also, if you hit a bug, let me know. Those routines are very new. – user14717 Feb 28 '18 at 05:37
  • http://www.boost.org/doc/libs/1_66_0/libs/math/doc/html/math_toolkit/double_exponential/de_exp_sinh.html – user14717 Feb 28 '18 at 05:39
  • @user14717 I changed tanh_sinh to exp_sinh, and this still works well, maybe a little faster. I should benchmark this when I find the time. – ThreeStarProgrammer57 Feb 28 '18 at 17:16