# Evaluating integral $F(r_2) = \frac{1}{r_2^n} \int_0^{r_2} r_1^n f(r_1)\ \mathrm{d}r_1$ without growing instability

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.

• 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. – Alone Programmer Feb 26 at 17:17
• 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$). – Alone Programmer Feb 26 at 18:14
• @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. – Firman Feb 26 at 20:40