Numerical integration with integrable singularity

Question

Suppose I have a numerical estimation of discrete samples of a smooth function $C(t)$ at $t = a, \dots, T = Na$ and I want to (numerically) compute the integral of $f(t) = \frac{C(t)}{\sqrt{t}}$. In particular, I am interested in computing $$ I = \int_{t=0}^T dt \, f(t) \approx a \sum_{n=1}^{N} f(na) \equiv I_N \,. $$ If I apply standard error analysis in the computation of integrals, essentially I know that the error on the integral is bounded by $$ \vert I-I_N \vert \leq \frac{M T}{12} a^2 $$ where $$ M \equiv \max{ \vert f''(t) \vert }_{t \in [0,T]} . $$ In this case, the function $f(t)$ has an integrable singularity due to $1/\sqrt{t}$, which should not affect the numerical computation of the integral, while $M$ seems to diverge.

How can I compute a bound on $\vert I-I_N \vert$ for this case?

Answer 1

(i) Your purported error bound is of course incorrect: consider e.g. $f(t)=t$.

(ii) To get rid of the singularity, make the substitution $u=\sqrt t$, so that $$\int_0^T dt\,f(t)=\int_0^T \frac{dt}{\sqrt t}\,C(t)=2\int_0^{\sqrt T} du\,C(u^2)=2\int_0^{\sqrt T} du\,g(u),$$ where $g(u):=C(u^2)$, and then approximate $\int_0^{\sqrt T} du\,g(u)$ by an integral sum.

---

Let us show that the error $|I-I_N|$ of your approximation will be on the order of $1/\sqrt N$ for large $N$ unless $C(0)=0$. Indeed, by rescaling, without loss of generality $T=1$ and then $a=1/N$. By the linearity of both $I$ and $I_N$ in $f$, it suffices to consider two cases:

• Case 1, when $C(0)=0$

• Case 2, when $C(t)=1$ for all $t$

In Case 1, the singularity disappears, and we have $|I-I_N|=O(1/N)$, the error bound for the rectangle method.

In Case 2, $$I-I_N=\int_0^1\frac{dt}{\sqrt t}-\frac1N\,\sum_{n=1}^N\frac1{\sqrt{n/N}}\sim-\frac{\zeta(1/2)}{\sqrt N},$$ and $\zeta(1/2)=-1.4603\dots$.

Thus, in Case 2 and whenever $C(0)\ne0$, we have $|I-I_N|\asymp1/\sqrt N$. $\quad\Box$

---

If the value of $C(0)$ is available as well, then one can use the "sample" values $C(n/N)$ much more effectively. Indeed, letting $f_0(t):=\frac{C(t)-C(0)}t$, one has a smooth function $f_0$ on $(0,1]$ with a bounded second derivative, and $$I=\int_0^1\frac{dt}{\sqrt t}\,C(t) =C(0)\int_0^1\frac{dt}{\sqrt t}+\int_0^1 dt\,f_0(t) =2C(0)+J \approx 2C(0)+J_N,$$ where $$J:=\int_0^1 dt\,f_0(t),\quad J_N:=\frac1N\,\sum_{n=1}^N f_0(n/N).$$ So, $$|I-(2C(0)+J_N)|=|J-J_N|=O(1/N),$$ again the error bound for the rectangle method without singularities. So, we get an $O(1/N)$ error, instead of $\asymp1/\sqrt N$.

If the value of $f_0(0+)=C'(0)$ is available as well, then, using high-order quadratures, we can get an $O(1/N^p)$ error for any natural $p$.