Establish the composite Simpson's 3/8 rule from basic Simpson's 3/8 rule over n subintervals and the error

Question

Establish the composite Simpson's 3/8 rule from basic Simpson's 3/8 rule over n subintervals and the error formula.

$$I(f) = \int_a^bf(x) dx \ \approxeq \ \frac{3h}{8}\left(f(a) \ + \ 3f\left(\frac{2a+b}{3}\right) \ + \ 3f\left(\frac{a+2b}{3}\right) \ + f(b)\right)$$ where $3h=b-a$, and the error of approximation $-3/80*f''''(c)*h^5$

I've spent a long time trying to understand how to derive the following result: $$ \int_a^b{f(x) dx} \approx \frac{3h}{8} \left[ f(x_0) + 3 \sum_{i=1}^{m}{\left(f(x_{3i-2})+f(x_{3i-1})\right)} + \ 2 \sum_{i=1}^{m-1}f(x_{3i}) + f(x_{3m}) \right] $$

I would appreciate it a lot if someone could give me an explanation and how to derive its error term.

Answer 1

We want $$\int_0^3f(x)dx\approx\sum_{i=0}^3w_if(i)$$ To be exact for polynomials of degree at most $3$ Then $$\begin{align}\int_0^3\frac{(x-1)(x-2)(x-3)}{(0-1)(0-2)(0-3)}dx&=\frac38=w_0\\ \int_0^3\frac{(x-0)(x-2)(x-3)}{(1-0)(1-2)(1-3)}dx&=\frac98=w_1\\ \int_0^3\frac{(x-0)(x-1)(x-3)}{(2-0)(2-1)(2-3)}dx&=\frac98=w_2\\ \int_0^3\frac{(x-0)(x-1)(x-2)}{(3-0)(3-1)(3-2)}dx&=\frac38=w_2\end{align}$$ Given a polynomial $p_3(i)=f(i)$ for $i\in\{0,1,2,3\}$ for any $x\in(0,3)\setminus\{0,1,2,3\}$ we can create $$e(y)=f(y)-p_3(y)-\frac{(y-0)(y-1)(y-2)(y-3)}{(x-0)(x-1)(x-2)(x-3)}\left(f(x)-p_3(x)\right)$$ So that $e(0)=e(1)=e(2)=e(3)=e(x)=0$ so after $4$ applications of Rolle's theorem we see that there is some $v(x)\in(0,3)$ such that $$e^{(4)}(v(x))=f^{(4)}(v(x))-\frac{4!}{(x-0)(x-1)(x-2)(x-3)}\left(f(x)-p_3(x)\right)=0$$ Thus for any $x\in(0,3)\setminus\{0,1,2,3\}$ there is some $v(x)\in(0,3)$ such that $$f(x)=p_3(x)+\frac{(x-0)(x-1)(x-2)(x-3)}{4!}f^{(4)}(v(x))$$ And also for any $x\in\{0,1,2,3\}$ by inspection. Thus $$\int_0^2f(x)dx=\int_0^3p_3(x)dx+\frac1{4!}\int_0^3(x-0)(x-1)(x-2)(x-3)f^{(4)}(v(x))dx$$ $$\int_0^2f(x)dx-\sum_{i=0}^3w_if(i)=\frac1{4!}\int_0^3(x-0)(x-1)(x-2)(x-3)f^{(4)}(v(x))dx$$ Now here we hit a snag. If $(x-0)(x-1)(x-2)(x-3)$ were negative throughout the interval of integration, we could say that $$\begin{align}&\max_{0\le x\le3}f^{(4)}(x)\int_0^3(x-0)(x-1)(x-2)(x-3)dx\\ &\quad\le \int_0^3(x-0)(x-1)(x-2)(x-3)f^{(4)}(x)dx\\ &\quad\le\min_{0\le x\le3}f^{(4)}(x)\int_0^3(x-0)(x-1)(x-2)(x-3)dx\\ &\quad=-\frac9{10}\min_{0\le x\le3}f^{(4)}(x)\end{align}$$ And then using the intermediate value theorem assuming continuous $f^{(4)}(x)$ we would arrive at $$\int_0^2f(x)dx-\sum_{i=0}^3w_if(i)=-\frac9{10}\frac1{4!}f^{(4)}(\xi)$$ For some $\xi\in(0,3)$, and this is the error quoted in the problem statement. Unfortunately $(x-0)(x-1)(x-2)(x-3)$ changes sign in the interval of integration, so the only easy result is that $$\begin{align}&\left|\int_0^3(x-0)(x-1)(x-2)(x-3)f^{(4)}(v(x))dx\right|\\ &\quad\le\int_0^3\left|(x-0)(x-1)(x-2)(x-3)\right|\max_{0\le y\le3}\left|f^{(4)}(y)\right|dx\\ &\quad=\max_{0\le x\le3}\left|f^{(4)}(x)\right|\int_0^3\left|(x-0)(x-1)(x-2)(x-3)\right|dx\\ &\quad=\frac{49}{30}\max_{0\le x\le3}\left|f^{(4)}(x)\right|\end{align}$$ So the best result you can get by this method is $$\left|\int_0^2f(x)dx-\sum_{i=0}^3w_if(i)\right|\le\frac1{4!}\frac{49}{30}\max_{0\le x\le3}\left|f^{(4)}(x)\right|=\frac{49}{720}\max_{0\le x\le3}\left|f^{(4)}(x)\right|$$ I have looked at several web sites and none of them establish the sharper number $\frac3{80}$ to my satisfaction or even quote the easier to establish $\frac{49}{720}$ so they are all $\text{IMHO}$ just talking out of their collective ass. There really are cases where assuming that error is $Ch^{n+1}f^{(n)}(\xi)$ isn't valid. So I am afraid that I can't explain how it comes about that everybody says that error term is correct either, nor have I seen anybody explain it.

Answer 2

Having establish the the quadrature formula by cubic interpolation of $f$ based at nodes $(a,f(a))$, $(a+h,f(a+h))$, $(a+2h,f(a+2+2h))$ and $(b,f(b))$ with $h=\frac{b-a}{3}$ as well as its write this quadrature formula in the form $$\int^b_af\approx\frac{3h}{8}\Big(f(a)+3f(a+h)+3f(a+2h)+f(b)\Big)$$ with error $e_h=-\frac{(b-a)^5}{6480}f^{(iv)}(c)$ for some $c\in (a, b)$.

Divide the interval $[a,b]$ in $3m$ pieces of same length $h=\frac{b-a}{3m}$. This defines a partition $x_k=a+\frac{b-a}{3m}$. On the each subinterval of the form $[x_{3k},x_{3(k+1)}]$, $k=0,\ldots ,m-1$ apply Simpson's rules based on cubic interpolation (the fist formula in the OP). That is $$\begin{align} \int^b_af=\sum^{m-1}_{k=0}\int^{x_{3(k+1)}}_{x_{3k}}f\tag{1}\label{one} \end{align}$$ and $$I_k:=\int^{x_{3(k+1)}}_{x_{3k}}f\approx\frac{3h}{8}\Big(f(x_{3k})+3f(x_{3k+1})+3f(x_{3k+2})+f(x_{3(k+1)})\Big)=A_k$$ Denoting the righthand side of \eqref{one} by $A$ we have that \begin{align} \frac{8}{3h}A&=\sum^{m-1}_{k=0}\Big(f(x_{3k})+3f(x_{3k+1})+3f(x_{3k+2})+f(x_{3(k+1)})\Big)\\ &=f(x_0)+\sum^{m-1}_{k=0}3\big(f(x_{3k+1})+f(x_{3k+2})\big) +2\sum^{m-1}_{k=1}f(x_{3k}) + f(x_{3m}) \end{align} The composite formula follows.

As for the error, each approximation $A_k$ to $I_k$ has error $-\frac{3}{80} h^5 f^{(iv)}(c_k)$ for some $c_k\in (x_{3k},x_{3(k+1)})$. The sum of all errors is $$E:=-\frac{3}{80}h^5\sum^{m-1}_{j=0}f^{(iv)}(c_k)$$ If we assume that $f\in\mathcal{C}^4([a,b])$, there there is $c\in(a, b)$ such that $f^{(iv)}(c)=\frac1m\sum^{m-1}_{k=0}f^{(iv)}(c_k)$ and so $$ E=-\frac{1}{80}(b-a)h^4f^{(iv)}(c) $$

---

With regards to the error in this quadrature, there are somewhat explicit formulas for the errors resulting from interpolation methods on evenly distributed nodes:

For $n\in\mathbb{N}$, $n\geq2$, the error of approximation to $\int^b_af$ in the quadrature method based on interpolating over the uniform partition $a_k=\frac{k}{n}(b-a)$, $0\leq k\leq n$, is

1. If $n$ is even and $f\in C^{(n+2)}[a,b]$, $$E_n=\frac{K_n}{(n+2)!}f^{(n+2)}(\eta),\quad K_n=\int^b_a x\prod^n_{k=0}(x-x_k)\,dx$$ for some $a<\eta<b$,
2. If $n$ is odd and $f\in C^{(n+1)}[a,b]$, $$E_n=\frac{K_n}{(n+1)!}f^{(n+1)}(\eta),\qquad K_n=\int^b_a\prod^n_{j=1}(x-x_j)\,dx$$

For Simpson's 3/8 rule, $n=3$ and so $K_3=\int^b_a (x-a)(x-\tfrac{2a+ b}{3})(x-\tfrac{a+2b}{3})(x-b)\,dx$

The simplification, though messy with pencil and paper, yields $E_3=-\frac{(b-a)^5}{6480}f^{(4)}(\eta)$ for some $a<\eta<b$.

Reference: Krylov, V. I. (translated by Stroud, A.), Approximate calculation of integrals, Dover Publications. 2005 pp 89-91