2
$\begingroup$

I have the following expression to be numerically integrated in a vector-based library (e.g. numpy, MATLAB, etc), $$ F(r_2) = \frac{1}{r_2^n} \int_0^{r_2} r_1^n f(r_1)\ \mathrm{d}r_1, $$ where $n$ is an integer that can have value up to $40$ and $r_1$ and $r_2$ can have values of many orders of magnitude (e.g. $r_2 \in [10^{-10},10^4]$).

The problem is, the numerical integration method that I use introduces a slight numerical instability which makes the integration result $\lim_{r_2\rightarrow 0}\int_0^{r_2}r_1^n f(r_1)\mathrm{d}r_1$ not exactly $0$, but something around $10^{-13}$. So when the integration result is divided by $r_2^n\sim10^{-10n}$ where $n \geq 3$, the result of $F(r_2\rightarrow 0)$ becomes really big.

I am looking for a way to allow the multiplication of $r_1^n/r_2^n$ takes place before the integration, so the numerical instability does not grow, or anything that prevents the numerical instability from the integration from growing really big.

The question is: how can I integrate the expression above numerically using numpy, MATLAB, or any vector-based library without growing the numerical instability?

P.S: I avoid for-loops if possible because I program in Python and for-loops can be really slow.

| cite | improve this question | |
$\endgroup$
  • $\begingroup$ Do you know $f(r_{1})$ analytically? If yes, What do you expect to have for $\lim_{r_{2} \rightarrow 0} \frac{1}{r_{2}^{n}} \int_{0}^{r_{2}} r_{1}^{n} f(r_{1}) dr_{1}$? By using L'Hôpital's rule you have: $$\lim_{r_{2} \rightarrow 0} \frac{1}{r_{2}^{n}} \int_{0}^{r_{2}} r_{1}^{n} f(r_{1}) dr_{1} = \lim_{r_{2} \rightarrow 0} \frac{r_{2}^{n} f(r_{2})}{n r_{2}^{n-1}} = 0$$ Unless $f(r_{2})$ does something funky at $r_{2} \rightarrow 0$, which I assume is not the case here. $\endgroup$ – Alone Programmer Feb 26 at 17:17
  • $\begingroup$ Also, something is not clear here. You said you have $r_{2} \in [10^{-10},10^{4}]$, but you are trying to evaluate $F(r_{2})$ for really small $r_{2}$ values close to zero ($r_{2} \rightarrow 0$). $\endgroup$ – Alone Programmer Feb 26 at 18:14
  • $\begingroup$ @AloneProgrammer using L'Hopital's rule is a good idea. Thanks! This is actually a numerical question, so in my case $r_2 = 10^{-10}$ is actually close to 0. $\endgroup$ – Firman Feb 26 at 20:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.