Your definition is
$(1)$ $f\in \mathscr R[a,b]$ with integral $\int_a^b f$ if and only if for each $\epsilon >0$ there exists $\delta >0$ such that for any partition $P$ with $\lVert P\rVert <\delta$ we have that $$\left|\int_a^b f-\sum(f,P)\right|<\epsilon$$
where $\sum(f,P)$ means the Riemann sum of $f$ with respect to $P$. The alternative definition is
$(2)$ $f\in \mathscr R[a,b]$ with integral $\int_a^b f$ if and only if for each $\epsilon >0$ there exists a partition $P_\epsilon$ such that for any partition $P$ a refinement of $P_\epsilon$ we have that $$\left|\int_a^b f-\sum(f,P)\right|<\epsilon$$
We would like to show one implies the other. First
$(1) \implies (2)$ Suppose $(1)$ holds. Since refinements can only decrease the mesh, and $\delta$ depends on $\epsilon$, the claim follows: for each $\epsilon>0$ take a partition $P_\epsilon$ with mesh $\delta'<\delta$ given by the above. Then any refiniment will have mesh at most $\delta'$ which will be less than $\delta$, and hence $(2)$ will hold.
$(2)\implies (1)$ This one is the tricky one. Here you can find a proof
An easy characterization of Riemann integrability is the following, which doesn't requiere that we know what the value of integral is:
$f:[a,b]\to\Bbb R$ is Riemann integrable over $[a,b]$ if and only if for each $\epsilon >0$ there exist a partition $P_\epsilon$ such that for any refinement $P$ of $P_\epsilon$ $$U(f,P)-L(f,P)<\epsilon$$
Moreover, this is equivalent to $$\overline{\int_a^b} f=\underline{\int_a^b}f$$
so if you're able to prove the above and evaluate any upper or lower integral, you're done.
In fact, this applies to the Riemann Stieljes integral whenever the integrator $\alpha$ is montone.
Let $\alpha:[a,b]\to\Bbb R$ be monotone. Let $f:[a,b]\to\Bbb R$. Then the following are equivalent:
$(1)$ The function $f$ is Riemann-Stieltjes integrable with respect to $\alpha$ over $[a,b]$
$(2)$ For each $\epsilon >0$ there exists a partition $P_\epsilon$ such that for any refinement $P$ of $P_\epsilon$ $$U(f,\alpha,P)-L(f,\alpha,P)<\epsilon$$
$(3)$ $$\overline{\int_a^b} f=\underline{\int_a^b}f$$
The idea is pretty clear, but we just have to write a lot, and use a lot of notation.
Definitions. A partition of $[a,b]$ is a finite sequence $\pi = (a_j)_{j=0}^J$ with
$a=a_0 < a_1 < \cdots < a_J = b$.
The norm of $\pi$ is $\|\pi\| = \max_{1 \le j \le J} |a_j-a_{j-1}|$.
A tagged partition of $[a,b]$ is a pair $(\pi,\gamma)$ where $\pi = (a_j)_{j=0}^J$
is a partition of $[a,b]$ and the finite squence $\gamma = (c_j)_{j=1}^J$ satisfies
$a_{j-1} \le c_j \le a_j$ for all $j$.
Let $f : [a,b] \to \mathbb R$. Let $\pi = (a_j)_{j=0}^J$ be a partition of $[a,b]$.
The upper sum is
$$
U(f,\pi) = \sum_{j=1}^J M_j(f)\;\cdot(a_j-a_{j-1})\qquad\text{where}\qquad
M_j(f) = \sup \{f(x) : a_{j-1} \le x \le a_j\} .
$$
The lower sum is
$$
L(f,\pi) = \sum_{j=1}^J m_j(f)\;\cdot(a_j-a_{j-1})\qquad\text{where}\qquad
m_j(f) = \inf \{f(x) : a_{j-1} \le x \le a_j\} .
$$
Let $(\pi,\gamma)$ be a tagged partition. The Riemann sum is
$$
R(f,\pi,\gamma) = \sum_{j=1}^J f(c_j)\cdot (a_j-a_{j-1})
$$
We say that $f$ is Riemann integrable, and its integral is $V$ iff
for every $\epsilon > 0$ there exists $\delta > 0$ such that for every
tagged partition $(\pi,\gamma)$ of $[a,b]$, if $\|\pi\| < \delta$, then
$$
\big |R(f,\pi,\gamma) - V \big| < \epsilon .
$$
Alternatively (from a Cauchy criterion), $f$ is Riemann integrable if and only if:
for every $\epsilon > 0$ there exists $\delta > 0$ such that for every
partition $\pi$ of $[a,b]$, if $\|\pi\| < \delta$ then
$$
U(f,\pi) - L(f,\pi) < \epsilon .
$$
If so, then the integral $V$ is the unique number with
$L(f,\pi) \le V \le U(f,\pi)$ for all partitions $\pi$ of $[a,b]$.
Theorem A
Let $f_1,f_1,\dots, f_K$ be functions on $[a,b]$ such that $f_1, f_2,\dots, f_K$
and the product $F = f_1 f_2\cdots f_K$
are all Riemann integrable.
Claim: for any $\epsilon > 0$ there exists $\delta > 0$ such that: if $\pi$ is a partition of $[a,b]$ and $\|\pi\| < \delta$ and $\gamma^{(k)}=(c_j^{(k)})_{j=1}^J$ , $1 \le k \le K$, are $K$ choices of tags so that for all $k$, the pair $(\pi,\gamma^{(k)})$ is a tagged partition of $[a,b]$, then
$$
\left|U(F,\pi) - \sum_{j=1}^J f_1(c_j^{(1)}) f_2(c_j^{(2)})\cdots
f_K(c_j^{(K)})\cdot (a_j-a_{j-1})\right| < \epsilon
\tag1$$
and
$$
\left|\sum_{j=1}^J f_1(c_j^{(1)}) f_2(c_j^{(2)})\cdots
f_K(c_j^{(K)})\cdot (a_j-a_{j-1})-L(F,\pi)\right| < \epsilon
\tag2$$
and
$$
U(F,\pi) - L(F,\pi) < \epsilon.
\tag3$$
Case 1: All $f_k \ge 0$.
So, $\epsilon > 0$ is given. Since all the functions $f_k$ are Riemann integrable, they are bounded. There is a single bound $A > 0$ with
$$
|f_k(x)|\le A \qquad\text{for}\qquad k=1,\dots,K,\qquad a\le x \le b .
$$
Let $$\epsilon' = \frac{\epsilon}{K A^{K-1}}.$$
Since $f_1,\dots,f_K$ are all Riemann integrable, there is $\delta > 0$ so that
for any partition $\pi = (a_j)_{j=0}^J$ of $[a,b]$, if $\|\pi\| < \delta$, then
$$
U(f_k,\pi) - L(f_k,\pi) < \epsilon'\qquad\text{for } k=1,\dots,K .
$$
Let $(\pi, \gamma^{(k)})$ be tagged partitions as above. Assume $\|\pi\| < \delta$. Write $\widetilde{M}_j = M_j(f_1)M_j(f_2)\cdots M_j(f_K)$ and
$\widetilde{m}_j = m_j(f_1)m_j(f_2)\cdots m_j(f_K)$.
Write
$$
\widetilde{U} = \sum_{j=1}^J \widetilde{M}_j\cdot(a_j-a_{j-1}),\qquad
\widetilde{L} = \sum_{j=1}^J \widetilde{m}_j\cdot(a_j-a_{j-1})
$$
Claim 1: $\widetilde{U}-\widetilde{L} < \epsilon$
Claim 2: All three of
$$
U(F,\pi), \quad L(F,\pi), \quad
\sum_{j=1}^J f_1(c_j^{(1)}) f_2(c_j^{(2)})\cdots
f_J(c_j^{(K)})\cdot (a_j-a_{j-1})
$$
are between $\widetilde{U}$ and $\widetilde{L}$.
Consequently, we obtain $(1), (2), (3)$, as required.
Proof of Claim 1.
The difference $\widetilde{U}-\widetilde{L}$ is written as a sum of $L$ terms. First, fix $j$ between $1$ and $J$. Then
\begin{align}
\widetilde{M}_j &= M_j(f_1)M_j(f_2)M_j(f_3)\cdots M_j(f_K)
\\ & \ge m_j(f_1)M_j(f_2)M_j(f_3)\cdots M_j(f_K)
\\ & \ge m_j(f_1)m_j(f_2)M_j(f_3)\cdots M_j(f_K)
\\ & \ge \cdots
\\ & \ge m_j(f_1)m_j(f_2)\cdots m_n(f_K) = \widetilde{m}_j
\end{align}
From each row to the next, one $M_j$ changes to an $m_j$. Then
\begin{align}
\widetilde{M}_j - \widetilde{m}_j &=
(M_j(f_1)-m_j(f_1)) M_j(f_2)M_j(f_3)\cdots M_j(f_K)
\\ &+ m_j(f_1)(M_n(f_2)-m_j(f_2)))M_j(f_3)\cdots M_j(f_K)
\\ &+ \cdots
\\&+ m_j(f_1) m_j(f_2)\cdots m_j(f_{K-1})(M_j(f_K)-m_j(f_K))
\\ &\le A^{K-1}\sum_{k=1}^K (M_j(f_k) - m_j(f_k))
\end{align}
Now multiply by $a_j-a_{j-1}$ and sum on $j$. The left side is
$$
\sum_{j=1}^J (\widetilde{M}_j - \widetilde{m}_j)\cdot(a_j-a_{j-1})
=\widetilde{U} - \widetilde{L}.
$$
The right side is
\begin{align}
A^{K-1}\sum_{k=1}^K &\sum_{j=1}^J (M_j(f_k) - m_j(f_k))\cdot(a_j-a_{j-1})
= A^{K-1}\sum_{k=1}^K \big(U(f_k,\pi) - L(f_k,\pi)\big)
\\ & < A^{K-1} K \epsilon' =\epsilon
\end{align}
This completes the proof of Claim 1.
Proof of Claim 2.
To prove $\widetilde{U} \ge U(F,\pi) \ge \widetilde{L}$ note that
$\widetilde{M}_j \ge M_j(F) \ge \widetilde{m}_j$.
To prove $\widetilde{U} \ge L(F,\pi) \ge \widetilde{L}$ note that
$\widetilde{M}_j \ge m_j(F) \ge \widetilde{m}_j$.
To prove
$$
\widetilde{U} \ge \sum_{j=1}^J f_1(c_j^{(1)}) f_2(c_j^{(2)})\cdots
f_J(c_j^{(K)})\cdot (a_j-a_{j-1}) \ge \widetilde{L}
$$
note
\begin{align}
\widetilde{M}_j &= M_j(f_1)M_j(f_2)\cdots M_j(f_K)
\\ &\ge f_1(c_j^{(1)}) f_2(c_j^{(2)})\cdots f_K(c_j^{(K)})
\\ &\ge m_j(f_1)m_j(f_2)\cdots m_j(f_K) = \widetilde{m}_j
\end{align}
multiply by $a_j-a_{j-1}$ and sum on $j$.
Case 2. arbitrary signs.
Write $f_k = f_k^+ - f_k^-$ in positive and negative parts. Then
$f_1 f_2\cdots f_k$ is a linear combination of $2^K$ such sums with nonnegative functions. Both
$$
\sum_{j=1}^J f_1(c_j^{(1)}) f_2(c_j^{(2)})\cdots
f_K(c_j^{(K)})\cdot (a_j-a_{j-1}) \qquad\text{and}\qquad
\int_a^b f_1(x)f_2(x)\cdots f_K(x)\;dx
$$
are linear combinations of $2^K$ terms of the same form, but involving only nonnegative functions. Theorem A for $f_1f_2\cdots f_k$ follows from those $2^k$ cases of Theorem A for nonnegative functions.
Best Answer
First example
Suppose we wish to calculate $\int_0^1 \sqrt{x} dx$. This is trivial if we already have the fundamental theorem of calculus in our disposal, but let's see what happens if we try to calculate this integral directly from the definition.
If we take a partition with equal intervals, that is, for every $n\in\mathbb{N}$ we take $x_k = \frac{k}{n}$ for $k\in \{0,1,\ldots,n\}$, we get the Riemann sum $$\sum_{k=1}^n \Delta x_k f(x_k) = \sum_{k=1}^n \frac{1}{n}\sqrt{\frac{k}{n}}=\frac{1}{n\sqrt{n}}\sum_{k=1}^n \sqrt{k}$$ and since there's no formula for the sum of square roots it appears we are stuck.
Since we're bothered by the the square roots, it would be a nice idea to try a partition with unequal intervals given by $x_k = \frac{k^2}{n^2}$, because now the Riemann sum becomes
$$\sum_{k=1}^n \Delta x_k f(x_k) = \sum_{k=1}^n (x_k - {x_{k-1}})\sqrt{\frac{k^2}{n^2}}=\sum_{k=1}^n \left(\frac{k^2}{n^2}-\frac{(k-1)^2}{n^2}\right)\frac{k}{n}=\frac{1}{n^3}\sum_{k=1}^n 2k^2 - k$$ so we got rid of the square roots. By the well known formulas for $\sum k$ and $\sum k^2$ (see Faulhaber's formula) we get that the Riemann sum equals $$\frac{1}{n^3}\left(2\cdot\frac{n(n+1)(2n+1)}{6}-\frac{n(n+1)}{2}\right)$$ so that taking the limit as $n\to\infty$ we get $\int_0^1 \sqrt{x} dx=\frac{2}{3}$.
Second example
Suppose we wish to find
$$\lim_{n\to\infty} \frac{1}{n^2}\sum_{k=1}^n (2k-1) \cos\left(\frac{\pi k^2}{2n^2}\right)$$
Noticing that $2k-1 = k^2 - (k-1)^2$, the limit equals
$$\lim_{n\to\infty} \sum_{k=1}^n \frac{k^2 - (k-1)^2}{n^2} \cos\left(\frac{\pi k^2}{2n^2}\right) = \frac{2}{\pi} \lim_{n\to\infty} \sum_{k=1}^n \left(\frac{\pi k^2}{2n^2} - \frac{\pi (k-1)^2}{2n^2}\right)\cos\left(\frac{\pi k^2}{2n^2}\right)$$ Denoting $x_k = \frac{\pi k^2}{2n^2}$ we see that this is $\frac{2}{\pi} \lim_{n\to\infty} \sum_{k=1}^n \left(x_k - x_{k-1}\right)\cos\left(x_k \right)$. Aside for $\frac{2}{\pi}$ this is the limit of a Riemann sum for $f(x)=\cos x$ on the interval $[0,\frac{\pi}{2}]$, so it equals $\int_0^{\frac{\pi}{2}} \cos x dx$ since we know $f(x)=\cos x$ is integrable. By the fundamental theorem of calculus $\int_0^{\frac{\pi}{2}} \cos(x)dx = 1$, so that the limit of the whole expression equals $\frac{2}{\pi}$.
Without knowing that Riemann sums of integrable functions converge to the same value (the integral) regardless of how we partition the domain of integration, this limit would probably be very hard to calculate.