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.

  • $\begingroup$ You are absolutely correct that this is a tricky point and in fact doesn't work in general. I have to work in a few minutes so I don't have time just now to post a really good answer, but if you are patient I might get time later to both show how difficult a question this is and to answer it affirmatively. $\endgroup$ – user5713492 May 2 '16 at 22:50
  • $\begingroup$ @user5713492 I would love an answer or some insight into the question if you find some time! When I first read the textbook, I guessed that it just worked out because $f^{(4)}(x)$ is assumed to be continuous - but playing with an example immediately showed that was not good enough. $\endgroup$ – Peter Woolfitt May 3 '16 at 1:03

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

  • $\begingroup$ I think I've seen that proof before, but it only seems to work for Simpson's rule and the trapezoidal rule. It doesn't seem to be applicable to Simpson's $\frac38$ rule or higher order Newton-Cotes formulas. $\endgroup$ – user5713492 May 3 '16 at 15:47
  • $\begingroup$ Thank you! Without much understanding of how all these formulas work, the use of the extended mean value theorem seems magical - though I suppose if the error is proportional to $f^{(4)}(\xi)$, it has to work out that way. $\endgroup$ – Peter Woolfitt May 5 '16 at 16:46
  • $\begingroup$ It doesn't have to work out that way. Try the same proof for Simpson's $\frac38$ rule and it goes nowhere. The proof really is magical in that sense. Of course my proof doesn't work for Simpson's $\frac38$ rule either, but it does for all Gaussian formulas inclusive of Gauss-Lobatto and Gauss-Radau. $\endgroup$ – user5713492 May 5 '16 at 20:06
  • $\begingroup$ Can you please clarify how the $ g(x)=F(x)-F(-x)-\frac{x}3(f(x)+4f(0)+f(-x)) $ is motivated? $\endgroup$ – sequence Jan 30 at 17:28
  • $\begingroup$ @sequence : It is the error of the numerical method. Replace $x$ by $h$ if you need it to look more familiar. $\endgroup$ – LutzL Jan 30 at 18:54

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.

  • $\begingroup$ Wow! Thank you so much. It has taken me a while to work through, but it was worth it. Your commentary with your work was extremely helpful in my understanding (and the derivation using $y'(0)$ is super sweet). $\endgroup$ – Peter Woolfitt May 5 '16 at 16:34
  • $\begingroup$ @Peter Woolfitt Thanks for the bounty! I am happy to see that you have gotten something of value from my posts. I always try to put my own spin on a problem, and not everyone appreciates that. $\endgroup$ – user5713492 May 22 '16 at 5:14

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
      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))+ &
      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
   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), &
   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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.