There are two kind of ways to extend the Riemann integral to unbounded sets. First is the one that you would probably guess first $$ \int_0^\infty f(x) dx = \lim_{b \to \infty}\int_0^b f(x)dx $$ as long as the integrals on the right exist for all $b$ and the limit on the right hand side exists. The other way is to say a sequence of unions of disjoint compact intervals $E_i = \cup_{k=1}^{n_i}[a_k,b_k]$ exhausts a set $E$ if $E_{i+1} \supseteq E_i$ and $\cup_{i=1}^\infty E_i = E$. (I will stick to unions of disjoint unions of intervals in order to avoid discussion about Peano-Jordan measurability, which you may look up if you wish.) Then you can define $$ \int_0^\infty f(x)dx = \lim_{n\to \infty} \int_{E_i} f(x)dx $$ whenever the limit on the right exists and is independent of the exhausting sequence $E_i$.

At first these definitions may look similar, but consider $$\int_0^\infty \frac{\sin(x)}{x} dx$$ If you take the limit as $b \to \infty$ as in the first definition, you get $\frac{\pi}{2}$. However, if you take the second definition, you will find that the integral does not exist because the integral of the positive portion and the negative portion are both $\infty$.

One thing that is easy to see is that if $f \geq 0$ is continuous, then the two kinds of integrals with both exist and equal eachother. Thus, of $f^+,f^-\geq 0$ are the positive and negative parts of $f$ and $$\int_0^\infty f^+(x)dx,\int_0^\infty f^-(x)dx < \infty$$ or equivalently $\int_0^\infty |f(x)| dx < \infty$, then the two kinds of integrals with both exist, be finite, and equal eachother.

Now to answer your question. Let's suppose that $f$ is continuous and $\int_0^\infty |f(x)|dx < \infty$. Let's investigate when $$\lim_{n\to\infty}\frac{1}{n}\sum_{k \geq 0} f\left(\frac{k}{n}\right) = \int_0^\infty f(x)dx.$$

Write $$ \frac{1}{n}\sum_{k \geq 0} f\left(\frac{k}{n}\right) = \sum_{m=0}^\infty \frac{1}{n}\sum_{k = mn}^{(m+1)n-1}f\left(\frac{k}{n}\right) = \sum_{m=0}^\infty S(f;1/n;[m,m+1]) $$ Where $S(f;\Delta x,[a,b])$ denotes the left Riemann sum of $f$ of width $\Delta x$ over the interval $[a,b]$. Since $f$ is continuous on $[m,m+1]$ for every $m$, we have $S(f;1/n;[m,m+1]) \to \int_{m}^{m+1} f(x)dx$ as $n \to \infty$. Then what we would like to say is that $$ \lim_{n\to \infty}\sum_{m=0}^\infty S(f;1/n;[m,m+1]) = \sum_{m=0}^\infty\lim_{n\to \infty} S(f;1/n;[m,m+1]) $$ $$ =\sum_{m=0}^\infty \int_{m}^{m+1} f(x)dx = \int_{0}^\infty f(x)dx $$ so what let's us justify these steps? The second equality is just definition, and since we assumed $\int_0^\infty |f(x)|dx < \infty$ we don't have to worry about the last equality since we know that $\{\cup_{m=0}^k[m,m+1]\}_k$ is an exhausting sequence of $[0,\infty]$ (single point overlap is okay, don't worry), so all that we have to justify is the first one, exchanging the limit and sum.

Let $S^*$ and $S_*$ represent the upper and lower Riemann sums respectively. Then of course $$ S_*(f;1/n;[m,m+1]) \leq S(f;1/n;[m,m+1]) \leq S^*(f;1/n;[m,m+1]) $$ and of more importance for this problem $$ |S_*(f;1/n;[m,m+1])|, |S^*(f;1/n;[m,m+1]) | \leq S^*(|f|;1/n;[m,m+1]) $$ Then one sufficient case to be able to exchange the sum and the limit is that for some $n$ $$\sum_{m=0}^\infty S^*(|f|;1/n;[m,m+1]) < \infty.$$ The justification for this is that you are bounding the partial sums by an integrable function and thus may apply Lebesgue's dominated convergence theorem. As a corollary, we get the friendly

Corollary. If $|f(x)| \leq g(x)$ where $g(x)$ is decreasing and $\int_0^\infty g(x)dx < \infty$, then $$\lim_{n\to\infty}\frac{1}{n}\sum_{k \geq 0} f\left(\frac{k}{n}\right) = \int_0^\infty f(x)dx. $$

Proof. $$\sum_{m=0}^N S^*(|f|;1;[m,m+1]) \leq S^*(g;1;[0,N+1]) \leq g(0) + \int_{0}^{N+1} g(x)dx \leq g(0) + \int_0^\infty g(x)dx$$ for all $N$. Thus $$\sum_{m=0}^\infty S^*(|f|;1;[m,m+1]) \leq g(0) + \int_0^\infty g(x)dx < \infty. $$

Also note, that you only really need $|f(x)| \leq g(x)$ for all $x>x_0$ for some fixed $x_0$.


For some $x\in[a,b]$, we have $$ \int_a^b f(x)\,\mathrm{d}x-f(a)(b-a)=f'(x)(b-a)^2 $$ Therefore, $$ \lim_{n\to\infty}n\left|\,\sum_{k=0}^\infty f\left(\frac kn\right)\frac1n-\int_0^\infty f(x)\,\mathrm{d}x\,\right| \lesssim\int_0^\infty\left|\,f'(x)\,\right|\,\mathrm{d}x $$ This says that if the variation of a function is finite, then the Riemann sums converge to the integral. This makes sense, because any function of bounded variation can be written as the difference of two monotonic decreasing functions.