$$ \newcommand{\norm}[1]{\left\Vert#1\right\Vert} \newcommand{\abs}[1]{\left\vert#1\right\vert} \newcommand{\set}[1]{\left\{#1\right\}} \newcommand{\Comp}{\mathbb C} \newcommand{\Real}{\mathbb R} \newcommand{\Nat}{\mathbb N} \newcommand{\of}[1]{\left ( #1 \right ) } \newcommand{\To}{\rightarrow} \newcommand{\eps}{\varepsilon} $$

Integral Calculus lecture notes: Simpson's rule (more advanced topic)

Eugene Kritchevski, Vanier College, last updated on August 18, 2016

Rolle's theorem and its generalization

Rolle's Theorem $\to$ Generalized Rolle's Theorem $\to$ Lagrange interpolating formula with remainder

Lagrange interpolation

If $(x_0,y_0)$ and $(x_1,y_1)$ are given points in the plane with $x_0\neq x_1$, then there exists a unique linear function $l(x)=mx+b$ such that $l(x_0)=y_0$ and $l(x_1)=y_1$. The parameters $m$ and $b$ can be found by solving the system of equations $$\begin{cases} mx_0+b =y_0 \\ mx_0+b =y_0 \\ \end{cases}$$ Taking the difference of the equations gives $m(x_1-x_0)=y_1-y_0$, which implies that $m=(y_1-y_0)/(x_1-x_0)$ and $b=y_0-mx_0=y_0-x_0(y_1-y_0)/(x_1-x_0)$. The linear function is then $$(1)\qquad l(x)=\frac{y_1-y_0}{x_1-x_0}x + y_0-x_0\frac{y_1-y_0}{x_1-x_0}.$$ The formula $(1)$ for $l(x)$ is correct but it is not particularly beautiful or easy to remember. Moreover, if we interchange the points $(x_0,y_0)$ and $(x_1,y_1)$, then $l(x)$ should not change, which is not immediately obious by looking at (1). Let us then rewrite $(1)$ in a nicer, more symmetric form.

... complete this part up to lagrange formula ...

Lagrange interpolation remainder formula
Suppose that $f$ is $n+1$ times differentiable on $[a,b]$. Let $x_0,x_1,\cdots,x_n$ be distinct nodes in $[a,b]$ and let $p$ be the interpolating polynomial for $f$ at these nodes. Is $x$ in $[a,b]$, then there is a $\xi$ in $[a,b]$ such that $$(*)\qquad f(x)=p(x)+\frac{f^{n+1}(\xi)}{(n+1)!}\prod_{k=0}^{n}(x-x_k)$$
Proof: Let $x$ be fixed. If $x$ coincides with a node, then $(*)$ holds for every $\xi$. Assuming $x$ is distinct from all the nodes, consider the function $$(**)\qquad g(t)=[f(t)-p(t)]\prod_{k=0}^{n}(x-x_k)-[f(x)-p(x)]\prod_{k=0}^{n}(t-x_k).$$ Then $g(x)=0$ and $g(x_k)=0$ for every node $x_k$. So $g$ has at least $n+2$ distinct zeros in $[a,b]$. Since $g$ is $n+1$ times differentiable, the generalized Rolles's theorem implies that $g^{n+1}(\xi)=0$ for some $\xi$ in $[a,b]$. Differentiaing $(**)$ $n+1$ times gives $$g^{n+1}(t)=[f^{n+1}(t)-0]\prod_{k=0}^{n}(x-x_k)-[f(x)-p(x)](n+1)!$$ and $(*)$ follows after replacing $t$ with $\xi$ and isolating $f(x)$.
Note that when $n=0$, we get the mean value thorem $$f(x)=f(x_0)+f'(\xi)(x-x_0)$$ Collide nodes to get Taylor... The other approach is to use (a clever) integration by parts to get Taylor with integral remainder form.

Simpson's Rule with optimal error bound

Simpson's rule
If $f^{(4)}$ is continuous on $[c-h,c+h]$, then $$(*)\qquad \int_{c-h}^{c+h}f(x)\,dx=\frac{h}{3}\left[f(c-h)+4f(c)+f(c+h)\right]-\frac{f^{(4)}(\xi)}{90}h^5$$ for some $\xi$ in [c-h,c+h].
Proof: We will first prove the formula in the special case when $c=0$ and $h=1$, as the general case will follow by a simple substitution. In this special case, the Simpson's rule is $$(1)\qquad \int_{-1}^{1}f(x)\,dx=\frac{1}{3}\left[f(-1)+4f(0)+f(1)\right]-\frac{f^{(4)}(\tau)}{90},$$ for some $\tau$ in $[-1,1]$. Given a small $\eps>0$, the Lagrange interpolation formula for the function $f$ at the four nodes $\set{-1,-\eps,\eps,1}$ gives $$(2)\qquad f(x)=P_\eps(x)+R_\eps(x), \qquad -1\leq x\leq 1,$$ where $P_\eps(x)$ is the interpolating polynomial $$(3)\qquad P_\eps(x)=f(-1)\frac{(x+\eps)(x-\eps)(x-1)}{(-1+\eps)(-1-\eps)(-1-1)}+f(1)\frac{(x+\eps)(x-\eps)(x+1)}{(1+\eps)(1-\eps)(1+1)}+ f(-\eps)\frac{(x+1)(x-1)(x-\eps)}{(-\eps+1)(-\eps-1)(-\eps-\eps)}+f(\eps)\frac{(x+1)(x-1)(x+\eps)}{(\eps+1)(\eps-1)(\eps+\eps)} $$ and the remainder term $R_\eps(x)$ is given by $$(4)\qquad R_\eps(x)=\frac{f^{(4)}(\xi)}{4!}(x-1)(x+1)(x-\eps)(x+\eps)$$ for some $\xi$ in $[-1,1]$.

Integrating both sides of $(2)$ from $-1$ to $1$ gives $$(5)\qquad\int_{-1}^{1}f(x)\,dx=\int_{-1}^{1}P_\eps(x)\,dx+\int_{-1}^{1}R_\eps(x)\,dx.$$ The integral $\int_{-1}^{1}P_\eps(x)\,dx$ can be worked out using $(3)$, and the result is $$\int_{-1}^{1}P_\eps(x)\,dx=[f(-1)+f(1)]\frac{\frac{1}{3}-\eps^2}{1-\eps^2}+[f(-\eps)+f(\eps)]\frac{2/3}{(1-\eps^2)}.$$ If we take the limit $\eps\to 0$, we get $$(6)\qquad\lim_{\eps\to 0}\int_{-1}^{1}P_\eps(x)\,dx= \frac{1}{3}[f(-1)+4f(0)+f(1)]$$ It follows from $(5)$ and $(6)$ that $$(7)\qquad\lim_{\eps\to 0}\int_{-1}^{1}R_\eps(x)\,dx=\int_{-1}^{1}f(x)\,dx - \frac{1}{3}[f(-1)+4f(0)+f(1)].$$ To complete the proof of $(1)$, we need to show that $$(8)\qquad\lim_{\eps\to 0}\int_{-1}^{1}-R_\eps(x)=\frac{f^{(4)}(\tau)}{90},$$ for some $\tau$ in $[-1,1]$. We denote by $m$ and $M$ the minimum and the maximum value of $f^{(4)}$ on $[-1,1]$. We have from $(4)$ that $$\int_{-1}^{1}-R_\eps(x)\,dx=\int_{1}^{1}\frac{f^{(4)}(\xi)}{4!}(1-x^2)(x^2-\eps^2)\,dx$$ and note that $(1-x^2)(x^2-\eps^2)\leq 0$ on the small interval $[-\eps,\eps]$ and $(1-x^2)(x^2-\eps^2)\geq 0$ on the rest of $[-1,1]$. As $\eps\to 0$ the contribution to the integral from $[-\eps,\eps]$ become negligible and therefore $$ (9)\qquad mI\leq\lim_{\eps\to 0}\int_{-1}^{1}-R_\eps(x) \leq M I,$$ where $$I=\int_{-1}^{1}\frac{1}{4!}(1-x^2)(x^2)\,dx=\frac{1}{90}.$$ Using $(9)$ and applying the intermediate value theorem to $f^{(4)}$, we get $$\lim_{\eps\to 0}\int_{-1}^{1}-R_\eps(x)=f^{(4)}(\xi)\frac{1}{90}$$ for some $\xi$ in $[-1,1]$, which proves $(8)$ and completes the proof of the Simpson's rule in the special case when $c=0$ and $h=1$. For the genral case, let $g(x)=f(c+hx)$. Applying $(1)$ to $g$ and and using the Chain rule proves $(*)$.