How to prove that the error in Simpson's rule is $- \frac{(b-a)^5}{2880}{{f^{(4)}(\zeta)}}$? [duplicate]
Solution 1:
Let $F$ be an anti-derivative of $f$, $F'=f$. W.l.o.g. $x_1=0$, set $x=h$ then we are interested in the error expression $$ g(x)=F(x)-F(-x)-\frac{x}3(f(x)+4f(0)+f(-x)). $$ This has derivatives \begin{alignat}{2} g'(x)&=\frac23(f(x)-2f(0)+f(-x))&&-\frac x3 (f'(x)-f'(-x)) \\ g''(x)&=\frac13(f'(x)-f'(-x))&&-\frac x3(f''(x)+f''(-x)) \\ g'''(x)&=&&-\frac x3(f'''(x)-f'''(-x)). \end{alignat} Consequently, by the extended mean value theorem $$ \frac{g(x)}{x^m} =\frac{g'(x_1)}{mx_1^{m-1}} =\frac{g''(x_2)}{m(m-1)x_2^{m-2}} =\frac {g'''(x_3)}{m(m-1)(m-2)x_3^{m-3}} $$ with $0<x_3<x_2<x_1<x$. Using $m=5$ this gives $$ \frac{g(x)}{x^5}=\frac{g'''(x_3)}{60x_3^2}=-\frac1{90}·\frac{f'''(x_3)-f'''(-x_3)}{2x_3}=-\frac1{90}·f^{(4)}(x_4). $$ with $|x_4|<x_3<x$. This results in the error formula $$ g(x)=-\frac1{90}·f^{(4)}(x_4)·x^5. $$ or after translating the initial simplifications back, \begin{align} \int_{x_0}^{x_2}f(x)\,dx &=\frac{h}{3}[f(x_0)+4f(x_1)+f(x_2)]-\frac{h^5}{90}f^{(4)}(\xi) \end{align} with $\xi\in(x_0,x_2)$
Solution 2:
Oh well, it looks like my brilliant counterexample is in the shop for repairs. But here's the idea of how we try to make this sort of error estimate, and it's the same idea as goes into proving the error estimate for the Taylor series. Simpson's rule is based on a polynomial $p(x)$ that matches the target function $f(x)$ at three points. Because we can always make the interval of integration $[-1,1]$ by a linear transformation, it is conventional to assume this has been done so we have $$y(x_0)=y(-1)=p(x_0)$$ $$y(x_1)=y(0)=p(x_1)$$ $$y(x_2)=y(1)=p(x_2)$$ We can construct a polynomial $q(x)$ of degree at most $3$ which also matches $y(x)$ at these three points and also at $x_3$ where $x_3\in(x_0,x_1)\cup(x_1,x_2)$ as follows: $$q(x)=p(x)+(y(x_3)-p(x_3))\prod_{k=0}^2\frac{(x-x_k)}{(x_3-x_k)}$$ Just evaluate $q(x)$ at each of the four points to assure yourself that this is so. Now consider the error in that approximation $$e(x)=y(x)-q(x)$$ We just now checked that $e(x)=0$ at $x_0$, $x_2$, and at two more points in $(x_0,x_2)$. By $3$ applications of Rolle's theorem we find that there is at least one point in $(x_0,x_2)$ such that $e^{\prime\prime\prime}(\xi)=y^{\prime\prime\prime}(\xi)-q^{\prime\prime\prime}(\xi)=0$. That is $$0=y^{\prime\prime\prime}(\xi)-3!\frac{y(x_3)-p(x_3)}{\prod_{k=0}^2(x_3-x_k)}$$ Since $x_3$ could have been any point in $[x_0,x_2]$ we have shown that $$y(x)=p(x)+\frac1{3!}y^{\prime\prime\prime}(\xi(x))\prod_{k=0}^2(x-x_k)$$ although we have to double check for $x\in\{x_0,x_1,x_2\}$. Notice that in general the point $\xi$ varies depending on $x$. But if we try to use this estimate to also estimate the error in Simpson's rule there is a complication. $$\int_{-1}^1y(x)dx=\int_{-1}^1p(x)dx+\frac1{3!}\int_{-1}^1y^{\prime\prime\prime}(\xi(x))\prod_{k=0}^2(x-x_k)dx$$ And at this point we would like to say something like 'the smallest the error term could be is' $$\frac1{3!}\max\left(y^{\prime\prime\prime}(\xi)\right)\int_{-1}^1\prod_{k=0}^2(x-x_k)dx$$ But we can't do that because the integrand changes sign in the interval so that if the wiggles in $y^{\prime\prime\prime}(\xi(x))$ line up with the wiggles in the integrand that could make the integral even smaller. The problem of the integrand in the error estimate changing sign in the interval of integration can be fixed by revamping our derivation of Simpson's rule from the start. We will attempt $4$-point osculation with a cubic polynomial $q(x)$ such that $$y(-1)=y(x_0)=q(-1)$$ $$y(0)=y(x_1)=q(0)$$ $$y^{\prime}(0)=y^{\prime}(x_1)=q^{\prime}(0)$$ $$y(1)=y(x_2)=q(1)$$ This time we construct a polynomial $r(x)$ of degree at most $4$ such that it satisfies the above conditions and also $r(x_3)=y(x_3)$ at a point $x_3\in(x_0,x_1)\cup(x_1,x_2)$. $$r(x)=q(x)+(y(x_3)-q(x_3))\frac{(x-x_0)(x-x_1)^2(x-x_2)}{(x_3-x_x)(x_3-x_1)^2(x_3-x_2)}$$ Again it can be checked to satisfy the $5$ conditions. Now the error in this approximation is $e(x)=y(x)-r(x)$ and has $4$ zeros in $[x_0,x_2]$. but this time $e^{\prime}(x)$ also has $4$ zeros in $(x_0,x_2)$: $3$ it gets from Rolle's theorem, and one more for the double root at $x=x_1$. So after $3$ more applications of Rolle's theorem we arrive at $$e^{(4)}(\xi)=0=y^{(4)}(\xi)-\left(y(x_3)-q(x_3)\right)\frac{4!}{(x_3-x_0)(x_3-x_1)^2(x_3-x_2)}$$ For some $\xi\in(x_0,x_2)$. So we have established $$y(x)=q(x)+\frac1{4!}(x-x_0)(x-x_1)^2(x-x_2)y^{(4)}(\xi(x))$$ And when we progress to $$\int_{-1}^1y(x)dx=\int_{-1}^1q(x)dx+\frac1{4!}\int_{-1}^1(x-x_0)(x-x_1)^2(x-x_2)y^{(4)}(\xi(x))dx$$ We really can say that $$\begin{align}-\frac4{15}\frac1{4!}\max\left(y^{(4)}(\xi)\right)&\le\frac1{4!}\int_{-1}^1(x-x_0)(x-x_1)^2(x-x_2)y^{(4)}(\xi(x))dx\\ &\le-\frac4{15}\frac1{4!}\min\left(y^{(4)}(\xi)\right)\end{align}$$ Because the integrand now has the same sign in $(x_0,x_2)$ and $$\int_{-1}^1(x-x_0)(x-x_1)^2(x-x_2)dx=\int_{-1}^1(x^4-x^2)dx=\frac25-\frac23=-\frac4{15}$$ By the intermediate value theorem we can conclude that $$\frac1{4!}\int_{-1}^1(x-x_0)(x-x_1)^2(x-x_2)y^{(4)}(\xi(x))dx= -\frac4{15}\frac1{4!}y^{(4)}(\xi)$$ for some $\xi\in(x_0,x_2)$ and then conclude that $$\int_{-1}^1y(x)dx=\int_{-1}^1q(x)dx-\frac1{90}y^{(4)}(\xi)$$ We just have to find $q(x)$ and integrate. Let $q(x)=a+bx+c(1-x^2)+d(x-x^3)$. Then $$q(0)=y(0)=a+c$$ $$q(-1)=y(-1)=a-b$$ $$q(1)=y(1)=a+b$$ $$q^{\prime}(0)=y^{\prime}(0)=b+d$$ It's easy to solve these for $$a=\frac{y(1)+y(-1)}2$$ $$b=\frac{y(1)-y(-1)}2$$ $$c=\frac{-y(1)+2y(0)-y(-1)}2$$ $$d=\frac{-y(1)+2y^{\prime}(0)+y(-1)}2$$ And $$\int_{-1}^1q(x)dx=2a+\frac43c=\frac13y(-1)+\frac43y(0)+\frac13y(1)$$ So even though we included $y^{\prime}(0)$ in our derivation its value was ignored in the eventual formula. This made it easier to get the error estimate for Simpson's rule and relied on the fact that the rule is really a Gauss-Lobatto formula that can actually ignore the values of the first derivatives of the function being approximated at all interior points. The fact that the error is not a bound, but actually attained for some $\xi\in(x_0,x_2)$ is useful for estimating the error in the composite Simpson's rule.
Hoping that this discussion at least gives you some insight on the leap of faith you are being asked to make as described in your question.
Solution 3:
My brilliant counterexample just got back from the shop, shiny new paint job and all, so it's time for a second answer. The idea of the exercise is to create an integration formula where it is not true that $$\int_{-1}^1y(x)dx=\int_{-1}^1q(x)dx+Cf^{(4)}(\xi)\tag{1}$$ For any $\xi\in(-1,1)$. To do this we are going to evaluate $y(x)$ for $x\in\{-1,-r,r,1\}$ and integrate the interpolation polynomial $q(x)$ from $-1$ to $1$ to get $$\int_{-1}^1q(x)dx=\frac{1-3r^2}{3(1-r^2)}(f(1)+f(-1))+\frac{2}{3(1-r^2)}(f(r)+f(-r))\tag{2}$$ If $(1)$ were true, we could just integrate $f_4(x)=x^4/24$ both ways and the difference between the exact and $4$-point methods would be $C$. If $r=\frac13$, this is just Simpson's $\frac38$ rule and I think there may be some theorem about Newton-Cotes formulas satisfying equation $(1)$, Runge phenomenon and all. But when $r=\frac1{\sqrt5}$, this is the $4$-point Gauss-Lobatto formula, exact for polynomials of degree up to $5$. So for some values of $r$ near that, we might anticipate equation $(1)$ to break down, and in fact for $0.434\le r\le0.481$ my program does demonstrate this failure. It shows that the formula is exact for polynomials of degree up to $3$, then applies the formula $(2)$ to $f_4(x)$, thus determining $C$ assuming the truth of equation $(1)$. Then it applies equation $(1)$ to polynomials $f_k(x)$, where $f_k^{(4)}(x)=T_{k-4}(x)$ for $5\le k\le10$, where $T_n(x)$ is the Chebyshev polynomial of degree $n$. We can then solve for $f^{(4)}(\xi)$ and if $\left|f^{(4)}(\xi)\right|>1$ it is indicative of a failure of equation $(1)$ because $\left|T_n(x)\right|\le1$ for $x\in[-1,1]$. This means that for $0.434\le r\le0.481$ the simple addition of $f^{(4)}(\xi_1)+f^{(4)}(\xi_2)=f^{(4)}(\xi_3)$ would be shown invalid because the error just isn't equal to $Cf^{(4)}(\xi)$ for some $C$ and any $\xi\in(-1,1)$.
module simpson
use ISO_FORTRAN_ENV, only: dp => REAL64
implicit none
type T
procedure(f0), nopass, pointer :: f
end type T
contains
function s1(f,r)
real(dp) s1, r
procedure(f0) f
s1 = (1-3*r**2)/(3*(1-r**2))*(f(1.0_dp)+f(-1.0_dp))+ &
2/(3*(1-r**2))*(f(r)+f(-r))
end function s1
function f0(x)
real(dp) f0, x
f0 = 1
end function f0
function if0(x)
real(dp) if0, x
if0 = x
end function if0
function f1(x)
real(dp) f1, x
f1 = x
end function f1
function if1(x)
real(dp) if1, x
if1 = x**2/2
end function if1
function f2(x)
real(dp) f2, x
f2 = x**2
end function f2
function if2(x)
real(dp) if2, x
if2 = x**3/3
end function if2
function f3(x)
real(dp) f3, x
f3 = x**3
end function f3
function if3(x)
real(dp) if3, x
if3 = x**4/4
end function if3
function f4(x)
real(dp) f4, x
f4 = x**4/24
end function f4
function if4(x)
real(dp) if4, x
if4 = x**5/5/24
end function if4
function f5(x)
real(dp) f5, x
f5 = x**5/120
end function f5
function if5(x)
real(dp) if5, x
if5 = x**6/6/120
end function if5
function f6(x)
real(dp) f6, x
f6 = x**6/180-x**4/24
end function f6
function if6(x)
real(dp) if6, x
if6 = x**7/7/180-x**5/5/24
end function if6
function f7(x)
real(dp) f7, x
f7 = x**7/210-x**5/40
end function f7
function if7(x)
real(dp) if7, x
if7 = x**8/8/210-x**6/6/40
end function if7
function f8(x)
real(dp) f8, x
f8 = x**8/210-x**6/45+x**4/24
end function f8
function if8(x)
real(dp) if8, x
if8 = x**9/9/210-x**7/7/45+x**5/5/24
end function if8
function f9(x)
real(dp) f9, x
f9 = x**9/189-x**7/42+x**5/24
end function f9
function if9(x)
real(dp) if9, x
if9 = x**10/10/189-x**8/8/42+x**6/6/24
end function if9
function f10(x)
real(dp) f10, x
f10 = 2*x**10/315-x**8/35+x**6/20-x**4/24
end function f10
function if10(x)
real(dp) if10, x
if10 = 2*x**11/11/315-x**9/9/35+x**7/7/20-x**5/5/24
end function if10
end module simpson
program test
use simpson
implicit none
real(dp) r
(f8),T(f9),T(f10)]
type(T) :: f(0:10)
type(T) :: jf(0:10)
integer i
real(dp) C
write(*,'(a)',advance='no') 'Enter r:> '
read(*,*) r
f = [T(f0),T(f1),T(f2),T(f3),T(f4),T(f5),T(f6),T(f7),T(f8),T(f9),T(f10)]
jf = [T(if0),T(if1),T(if2),T(if3),T(if4),T(if5), &
T(if6),T(if7),T(if8),T(if9),T(if10)]
do i = 0, 3
write(*,*) i,s1(f(i)%f,r),jf(i)%f(1.0_dp)-jf(i)%f(-1.0_dp)
end do
i = 4
write(*,*) i,s1(f(i)%f,r),jf(i)%f(1.0_dp)-jf(i)%f(-1.0_dp)
C = (jf(i)%f(1.0_dp)-jf(i)%f(-1.0_dp)-s1(f(i)%f,r))
write(*,*) C
do i = 5, 10
write(*,*) i,s1(f(i)%f,r),jf(i)%f(1.0_dp)-jf(i)%f(-1.0_dp)
write(*,*) (jf(i)%f(1.0_dp)-jf(i)%f(-1.0_dp)-s1(f(i)%f,r))/C
end do
end program test