Order of accuracy for trapezoidal integration

Question

I have a question regarding the order of accuracy for certain function using the trapezoidal formula.

I know that from theory the formula is second order accurate, but when working with matlab I get different answers.

I calculated the order of accuracy by plotting the logarithm of the error vs the logarithm of the step size $H$.

When using functions which are infinitely differentiable, for example $\sin(x)$ or $e^x$, I get a slope of $-2$, as expected.

But when using the formula $\arctan(\sqrt{x})$ I get a slope of about -1.5.

Why does the order of accuracy depend on the function that is integrated?

Answer 1

Possible solution: The error term involves a second derivative of the function being integrated. And your example is not differentiable at 0. (So I suggested testing it on an interval where the second derivative exists).

Answer 2

If $f$ has at least two bounded derivative in the interval, the trapezium rule has error bounded by $1/N^2$ times the supremum of the second derivative. However, for $\sqrt{x}$ this is not the case on $(0,1)$, as reflected in your data.

It is possible to derive a precise estimate for the error, but it is not an easy calculation. Instead, we can get away with some simple inequalities to give the order of the error bound exactly, while not worrying too much about the constant. We work over the interval $(0,1)$, because the calculations are easier, but the same will occur for any interval with $0$ as its lower endpoint. We also write $1/H=N$, and take this to be an integer, since it makes the sums easier to work with.

The approximation with $N$ intervals can be written as $$ \frac{1}{N}\sum_{k=1}^{N} \sqrt{\frac{k}{N}} - \frac{1}{2N} = N^{-3/2} \sum_{k=1}^N \sqrt{k} - \frac{1}{2N}. $$

We now examine the sum, $S_N = \sum_{k=1}^N \sqrt{k}$. We want an estimate for this that is sufficiently good that we can show that there is an $O(N^{-3/2})$ term in the difference between the integral and the approximation. We therefore need something with leading term $\frac{2}{3} N^{3/2}$, and a term that cancels the $1/(2N)$, and want to prove tha the difference between this and $S_N$ is bounded above and below by constants with the same sign. So consider $$ T_N = \frac{2}{3}N^{3/2} + \frac{1}{2} N^{1/2}. $$ In particular, $S_0=T_0=0$, and $$ (S_{N+1}-T_{N+1})-(S_N-T_N) = \frac{1}{2}(\sqrt{N}+\sqrt{N+1})+\frac{2}{3}(N^{3/2}-(N+1)^{3/2}). $$ Expanding the right-hand side in a series for large $N$, we find the right-hand side is $$ -\frac{1}{48}N^{-3/2}(1+O(N^{-2})), $$ which is bounded away from zero for sufficiently large $N$, so in particular, the difference is negative and decreasing for $N>0$. Also, $\sum_{n=1}^{\infty} n^{-3/2}$ converges, so the difference is bounded ($ \sum_{k=1}^N [(S_k-T_k)-(S_{k-1}-T_{k-1})] = S_N-T_N$ by telescoping). Thus there are $0<A<B$ so that $$ -A < S_N -T_N < -B $$ for sufficiently large $N$.

From this, we find that the error is $$ \int_0^1 \sqrt{x} \, dx - \left( N^{-3/2} \sum_{k=1}^N \sqrt{k} - \frac{1}{2N} \right) = \frac{2}{3} + \frac{1}{2N} - N^{-3/2} S_N = N^{-3/2} ( T_N - S_N ), $$ which we have just determined is bounded between $AN^{-3/2}$ and $BN^{-3/2}$. Hence the error is definitely $O(N^{-3/2})$, or $O(H^{3/2})$ in terms of the step size.

A similar analysis will work on $x^{\alpha}$ for any $0<\alpha<1$.

Answer 3

You can compute the error in the trapeze method by using integration by parts, let $f:(0,1) \to \Bbb R$.

Then the error is the sum of local errors: $$E=\sum_{k=0}^{n-1} \left(\int_\frac{k}{n}^\frac{k+1}{n}\left(f(t)-\frac{f\left(\frac{k}{n}\right)+f\left(\frac{k+1}{n}\right)}{2} \right) dt \right)$$ integrating by parts: $$\int_\frac{k}{n}^\frac{k+1}{n}\left(f(t)-\frac{f\left(\frac{k}{n}\right)+f\left(\frac{k+1}{n}\right)}{2} \right) =0-\int_\frac{k}{n}^\frac{k+1}{n}f'(t) \left(t-\frac{k+\frac{1}{2}}{n} \right) dt$$ and integrating by parts again: $$\int_\frac{k}{n}^\frac{k+1}{n}\left(f(t)-\frac{f\left(\frac{k}{n}\right)+f\left(\frac{k+1}{n}\right)}{2} \right) =0+\int_\frac{k}{n}^\frac{k+1}{n}f''(t) \left(t-\frac{k}{n} \right)\left(t-\frac{k+1}{n} \right) dt$$

---

If $f''$ is bounded the idea to recover the standard bound is to write: $$|E|=\left|\sum_{k=0}^{n-1}\int_\frac{k}{n}^\frac{k+1}{n}f''(t) \left(t-\frac{k}{n} \right)\left(t-\frac{k+1}{n} \right) dt \right|\leq \|f''\|_{\infty} \left|\sum_{k=0}^{n-1}\int_\frac{k}{n}^\frac{k+1}{n} \left(t-\frac{k}{n} \right)\left(t-\frac{k+1}{n} \right) dt \right|=\|f''\|_{\infty} n \cdot \frac{1}{6 n^3}$$

---

Here we can write, as $f''$ is not bounded: $$E=\int_0^\frac{1}{n}f''(t)\left(t-0\right)\left(t-\frac{1}{n} \right)dt+\sum_{k=1}^{n-1}\int_\frac{k}{n}^\frac{k+1}{n}f''(t) \left(t-\frac{k}{n} \right)\left(t-\frac{k+1}{n} \right) dt$$ so for $f(t)=\sqrt{t}$ (in fact or any function that looks like $\sqrt{t}$ near $0$ and is smooth far from $0$): $$|E|\leq\int_0^\frac{1}{n} \left|\frac{1}{4t^\frac{3}{2}}t\left(\frac{1}{n} -t\right)\right|dt+\sum_{k=1}^{n-1} \frac{1}{4n^2}\int_\frac{k}{n}^\frac{k+1}{n}|f''(t)| dt$$ so as $f''$ is of constant sign: $$|E|\leq\frac{1}{3n^\frac{3}{2}}+\frac{1}{4n^2}\int_\frac{1}{n}^1 |f''(t) |dt=\frac{1}{3n^\frac{3}{2}}+\frac{1}{4n^2}\left|f'(1)-f'\left(\frac{1}{n} \right) \right|\leq\frac{1}{3n^\frac{3}{2}}+\frac{1}{4n^2}+\frac{1}{2n^\frac{3}{2}}$$

so the $\frac{3}{2}$ rate of convergence. Note that the problem lies near $t=0$, as for example the part between $0$ and $\frac{1}{n}$ contribute with a $n^\frac{-3}{2}$ in itself.