Incorrect result during evaluating Riemann sum for $\frac{1}{\sqrt{1 + x^2}}$ using Python

Question

While solving an exercise problem from here, I stumbled upon a problem in the evaluation of:

$$ \int_0^u \frac{1}{\sqrt{1 + x^p}} dx \approx \sum_{k=1}^N {\frac{1}{1 + (x_k^*)^p} \Delta x} $$

where $ \Delta x = \frac{u}{N} $ and $ x_k^* = \frac{x_k + x_{k - 1}}{2} $ for endpoints $ x_k = k \Delta x $ using Python.

Note: $x_k^* = \frac{x_k + x_{k - 1}}{2} = \frac{\Delta x(2k - 1)}{2}$

For a test case, plugging $u = 1, p = 2$ yields $\int_0^1 \frac{1}{\sqrt{1 + x^2}} dx \approx 0.8813$

---

Have a look at my code:

def sqrt_integral(N, upper, p):
    "Evaluates the integral of 1/sqrt(1 + x^p) from 0 to u."

    dx = upper / N
    term = [dx / (1 + ((2*k - 1)*dx / 2)**p) for k in range(1, N + 1)]

    summation = sum(term)
    return summation


print(sqrt_integral(10000, 1, 2))

The term is a list comprehension that evaluates each term in the form:

$$ \sum_{k=1}^N \frac{\Delta x}{1 + \left(\frac{\Delta x(2k - 1)}{2}\right)^p} $$

Now, I don't understand why it yields $\approx 0.7853$.

Answer 1

The formula for the Riemann sum should be

$$\int_{x=a}^b f(x) \, dx = \lim_{N \to \infty} \sum_{k=1}^N f(x_k^*) \Delta x, \tag{1}$$ where $\Delta x = \frac{b-a}{N}$, $x_k^*$ is defined as you have done, and $x_k = a + k \Delta x$.

In your case, $a = 0$, $b = 1$, and $$f(x) = \frac{1}{\sqrt{1+x^2}}.$$ This makes the sum $(1)$

$$\lim_{N \to \infty} \sum_{k=1}^N \frac{1}{\color{red}{\sqrt{1+(x_k^*)^2}}} \Delta x. \tag{2}$$

Try it out. You will find that the result numerically matches.