Combining error terms in Simpson's rule

Question

My numerical analysis textbook (Burden and Faires) derives Simpson's rule as

$$\begin{align} \int_{x_0}^{x_2}f(x)\,dx&=2hf(x_1)+\frac{h^3}{3}f''(x_1)+\frac{h^5}{60}f^{(4)}(\xi_1) \\&=\frac{h}{3}[f(x_0)+4f(x_1)+f(x_2)]-\frac{h^5}{12}\left[\frac{1}{3}f^{(4)}(\xi_2)-\frac{1}{5}f^{(4)}(\xi_1)\right] \end{align}$$ for some $\xi_1,\xi_2\in(x_0,x_2)$. Here $x_0=x_1-h$ and $x_2=x_1+h$.

Note that $\xi_1$ comes from the error term when integrating the Taylor series of $f(x)$ while $\xi_2$ comes from the error term in rewriting $f''(x_1)$.

How can we show that we can replace $\xi_1$ and $\xi_2$ in the above formula with some $\xi\in(x_0,x_2)$?

My textbook leaves this part as an exercise, but I am unconvinced the exercise demonstrates what we want. The following is the exercise verbatim:

"Derive Simpson's rule with error term by using $$\int_{x_0}^{x_2}f(x)\,dx=a_0f(x_0)+a_1f(x_1)+a_2f(x_2)+kf^{(4)}(\xi)$$ Find $a_0$, $a_1$, and $a_2$ from the fact that Simpson's rule is exact for $x^n$ when $n=1$, $2$, and $3$. Then find $k$ by applying the integration formula with $f(x)=x^4$."

Performing this exercise, we certainly get the correct coefficients, but I don't see how this tells us we can combine $\xi_1$ and $\xi_2$ - it seems that we just assumed it.

Answer 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)$

Answer 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.

Answer 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

Answer 4

The question came up if the mean value strategy also works for other methods, on the example of the 3/8 method. Thus take $$ g(x)=F(x)-F(-x)-\frac{x}{4}(f(x)+3f(x/3)+3f(-x/3)+f(-x)) $$ as difference of the integral of $f$ over $[-x,x]$, expressed via a primitive function $F'=f$, and its approximation via a finite number of values.

Compute derivatives \begin{align} g'(x)&=\frac{3}{4}(f(x)-f(x/3)-f(-x/3)+f(-x)) -\frac{x}{4}(f'(x)+f'(x/3)-f'(-x/3)-f'(-x)) \\ g''(x)&=\frac12(f'(x)-f'(x/3)+f'(-x/3)-f'(-x)) -\frac{x}{12}(3f''(x)+f''(x/3)+f''(-x/3)+3f''(-x)) \\ g'''(x)&=\frac14(f''(x)-f''(x/3)-f''(-x/3)+f''(-x)) -\frac{x}{36}(9f'''(x)+f'''(x/3)-f'''(-x/3)-9f'''(-x)) \\ g^{(4)}(x)&=-\frac29(f'''(x/3)-f'''(-x/3)) -\frac{x}{108}(27f^{(4)}(x)+f^{(4)}(x/3)+f^{(4)}(-x/3)+27f^{(4)}(-x)) \end{align} Applying the chained mean value theorem gives $$ \frac{g(x)}{x^5}=...=\frac{g^{(4)}(x_3)}{5!x_3}=-\frac{4}{5!\cdot 27}f^{(4)}(x_4)-\frac1{5!\cdot4}(f^{(4)}(x_3)+f^{(4)}(-x_3))-\frac1{5!\cdot108}(f^{(4)}(x_3/3)+f^{(4)}(-x_3/3)) $$ This is now a linear combination of 4th derivatives inside the interval, $|x_4|<x_3/3$, with all coefficients having the same negative sign. Thus we can apply the mean value theorem of summation to replace all the points with the same intermediate point $$ \frac{g(x)}{x^5} = -\frac1{5!\cdot54} \left(8+27+1\right) f^{(4)}(\xi) = -\frac1{180} f^{(4)}(\xi) $$ Remember that the length of the integration interval was $2x$, so that if that is replaced by $h$, an appropriate power of $2$ will occur in the error term.