Integral expansion of Riemann sum

Question

Consider an equidistant partition of the interval $[0,1]$ with a step size of $\Delta$ and $N+1$ grids $\{x_i\}_{i=0:N}$ such that $N\Delta= 1$. For a Riemann integrable function $f$, we have the following standard result for the left Riemann sum: $$\sum_{i=0}^Nf(x_i)\cdot \Delta = \int_0^1 f(x)dx + O(\Delta).$$ I wonder whether I can express the $O(\Delta)$ in RHS of the above equation with a more precise result. For example, suppose I take the Taylor expansion of $f$ at any point $x_i\in[0,1]$: $$f(x)=f(x_i)+\frac{1}{2}f'(x_i)(x-x_i)+... $$ Integrating both sides gives: $$\int_0^1f(x)dx = f(x_i)+(\frac{1}{4}-\frac{x_i}{2})f'(x_i)+ ... $$ I can now take an average of the above equation over each $i$: $$\int_0^1f(x)dx =\sum_{i=0}^Nf(x_i)\Delta +\sum_{i=0}^N (\frac{1}{4}-\frac{x_i}{2})f'(x_i) \Delta+ ... $$ The second term then converges again as a Riemann sum to $\int_0^1(\frac{1}{4}-\frac{x}{2})f'(x)dx$. Thus it is tempting to write: $$\sum_{i=0}^Nf(x_i)\cdot \Delta \approx \int_0^1 f(x)dx + \int_0^1(\frac{x}{2}-\frac{1}{4})f'(x)dx.$$ Does this approximation work in the sense that it reduces the order of the residual term? Can we also generate this expansion to higher order by taking more terms from the Taylor expansion? I'd be grateful if anyone can provide me with some hints or directions.