It's not a rounding error, you just don't pick the best choice for integration points. Change the initialisation to
x[j]=(j+0.5)*dx;
so that you take the midpoint of each integration strip to calculate the value of the integrand. If you always take the left or right endpoint, you will get a systematically too large error for monotonic functions.
If you approximate the integral of a sufficiently smooth function f
by a Riemann sum,
b n
∫ f(x) dx ≈ ∑ f(y_k)*(b-a)/n
a k=1
the choice of y_k
in the interval [x_(k-1), x_k] = [a+(k-1)*(b-a)/n, a+k*(b-a)/n]
influences the error and the speed of convergence. Writing
f(x) = f(y_k) + f'(y_k)*(x-y_k) + 1/2*f''(y_k)*(x-y_k)² + O((x-y_k)³)
in that interval, you find that
x_k x_k x_k
∫ f(x) dx = f(y_k)*(b-a)/n + f'(y_k)* ∫ (x-y_k) dx + 1/2*f''(y_k) * ∫ (x-y_k)² dx + O(1/n^4)
x_(k-1) x_(k-1) x_(k-1)
= f(y_k)*(b-a)/n + 1/2*f'(y_k)*(b-a)/n*((x_k-y_k)-(y_k-x_(k-1))) + O(1/n³)
and the first and largest error term with respect to the approximation f(y_k)*(b-a)/n
vanishes for
y_k = (x_k + x_(k-1))/2
giving you an overall O(1/n³) error for that strip, and a total O(1/n²) error for the entire Riemann sum.
If you choose y_k = x_(k-1)
(or y_k = x_k
), the first error term becomes
±1/2*f'(y_k)*[(b-a)/n]²
leading to an O(1/n) total error.