Code:

import sympy

x = sympy.symbols('x')
f = 3*x + 1
a, b = 0, 2

# Exact area using definite integral
exact_area = sympy.integrate(f, (x, a, b))

# Riemann sum approximation (Right Riemann Sum with n subintervals)
n = sympy.symbols('n', integer=True, positive=True)
delta_x = (b - a) / n
x_i = a + sympy.symbols('i') * delta_x
f_xi = f.subs(x, x_i)
riemann_sum_expr = sympy.Sum(f_xi * delta_x, (sympy.symbols('i'), 1, n))
riemann_sum_simplified = riemann_sum_expr.doit().simplify()

# Limit as n approaches infinity
limit_area = sympy.limit(riemann_sum_simplified, n, sympy.oo)

print(f"{exact_area=}")
print(f"{delta_x=}")
print(f"{f_xi=}")
print(f"{riemann_sum_simplified=}")
print(f"{limit_area=}")

Output:

exact_area=8
delta_x=2/n
f_xi=6*i/n + 1
riemann_sum_simplified=8 + 6/n
limit_area=8