Composite trapezoidal rule's error term representation has $k_1h^2 + k_2 h^4 +... $ to apply Romberg integration

Question

I am trying to apply the Romberg method but when trying to evaluate the Error, I can't eliminate the appearance of the term $\frac{1}{12}f^{(3)}(\psi)h^3$, When I wish to get $E = K_1h^2 +K_2h^4 +K_3 h^6...$

What I tried to do was:

For each interval $[x_i, x_{i+1}]$ the error element is of the form: $\frac{1}{12}h^{3}f^{(2)}(s_i)$ when $s_i\in[x_{i},x_{i+1}]$ and $h=[x_{i+1}-x_i]$ (the equi-distanced partition). Using Taylor polynomial at $x_i$ $f^{(2)}(s_i) = f^{(2)}(x_i) + f^{(3)}(x_i)(s_i-x_i) + \frac{1}{2}f^{(4)}(\psi_i)(s_i-x_i)^2$.

So $|E_i| \le |\frac{1}{12}h^3[ f^{(2)}(x_i) + f^{(3)}(x_i)(h) + \frac{1}{2}f^{(4)}(\psi_i)(h)^2]|$

Thus the error sum is bounded by the absolute value of: $\frac{1}{12}h^2 \cdot h\cdot \sum [ f^{(2)}(x_i) + f^{(3)}(x_i)(h) + \frac{1}{2}f^{(4)}(\psi_i)(h)^2]$ Now applying intermediate value theorem, we get (for example for $\sum f^{(2)}(x_i)$: $\min\{f(x_i)\} \le h\cdot\sum_{i=0} ^{n-1} f^{(2)}(x_i) =\frac{1}{n}\cdot\sum_{i=0} ^{n-1} f^{(2)}(x_i) \le \max\{f(x_i)\} $ .

Thus the error can be represented: $\frac{1}{12}h^2 \cdot [ f^{(2)}(\varphi) + f^{(3)}(\psi)(h) + \frac{1}{2}f^{(4)}(\chi)(h)^2]$ for $\varphi , \psi,\chi \in [x_0,x_n]$ derived from the intermediate value theorem (and the assumption $f^{(4)}$ is continuous. So the error is bounded by $\frac{1}{12} \cdot [ f^{(2)}(\varphi)h^2 + f^{(3)}(\psi)h^3 + \frac{1}{2}f^{(4)}(\chi)h^4]$.

In Burden's book it tells that there is a detailed explanation in Ralston and Rabinowitz book, at pages 136-140. But I don't have access to this book for the next few days.

Answer 1

It's about even/odd symmetry. Translate that interval $[x_i,x_{i+1}]$ to $\left[-\frac{h}{2},\frac{h}{2}\right]$ and split $f$ into even and odd parts; $f(x)=g(x)+q(x)$ on that interval with $g$ even and $q$ odd.

Then $\int_{-h/2}^{h/2}f(x)\,dx=\int_{-h/2}^{h/2}g(x)\,dx$ since the integral of an odd function is zero, and the trapezoid rule estimate is $\frac{h}{2}\left(g(\frac h2)+q(\frac h2)+g(-\frac h2)+q(-\frac h2)\right) = hg(\frac h2)$.

Now, since we've eliminated the odd part, expand $g$ by its power series $$g(x)=g(0)+\frac{g''(0)}{2}x^2+\frac{g^{(4)}(0)}{4!}x^4+\frac{g^{(6)}(0)}{6!}x^6+\cdots$$ Because we're working with an even function, there aren't any odd derivatives involved. Integrate that power series, and $$\int_{-h/2}^{h/2}g(x)\,dx = hg(0)+\frac{h^3}{24}g''(0)+\frac{h^5}{5!\cdot 2^4}g^{(4)}(0)+\frac{h^7}{7!\cdot 2^6}g^{(6)}(0)+\cdots$$ compared to a trapezoid rule estimate of $$hg\left(\frac h2\right) = hg(0)+\frac{h^3}{8}g''(0)+\frac{h^5}{4!\cdot 2^4}g^{(4)}(0)+\frac{h^7}{6!\cdot 2^6}g^{(6)}(0)+\cdots$$ All of those even derivatives of $g$ at zero are equal to the even derivatives of $f$ there, of course.

On an arbitrary interval, we would get the same effect, with a lot more clutter, by expanding the Taylor series around the center of the interval. That's the key to it. It's just clearer to move that symmetry point to the origin for illustration purposes.

Answer 2

We remark that the trapezoidal rule on the interval $[a,b]$ with $n+1$ data points is exact for $$f_{2n+1}(x)=\left(x-\frac{a+b}2\right)^{2n+1}$$ because the function is odd about $x=(a+b)/2$ so both the trapezoidal rule result and the integral are identically zero. Then if Romberg's rule is exact for polynomials of degree at most $2n$, since we can write any polynomial of degree at most $2n+1$ as $p_{2n+1}(x)=a_{2n+1}f_{2n+1}(x)+p_{2n}(x)$ where $a_{2n+1}$ is a constant and $p_{2n}(x)$ is a polynomial of degree at most $2n$, it is exact for $p_{2n+1}(x)$ because it's a linear combination of trapezoidal rules which all approximated $$\int_a^bf_{2n+1}(x)dx\approx0$$ Now this doesn't tell us how big the error in Romberg's rule is, but we can find out the error for a polynomial of lowest degree for which it isn't exact. The general idea is to develop the theory of Bernoulli polynomials and then establish the Euler-MacLaurin summation formula, hence the trapezoidal rule with correction terms, and then apply Romberg's rule.

The Bernoulli polynomials are defined such that $$B_0(x)=1$$ $$B_n^{\prime}(x)=nB_{n-1}(x),\quad n>0$$ $$B_n(0)=B_n(1),\quad n>1$$ Then we can integrate to get $$B_1(x)=1\cdot\int B_0(x)dx=\left(x-\frac12\right)+c_1$$ $$B_2(x)=2\int B_1(x)dx=\left(x-\frac12\right)^2+2c_1\left(x-\frac12\right)+c_2$$ Applying the boundary condition $$B_2(1)-B_2(0)=2c_1=0$$ Hopefully it's obvious that each $c_{2n+1}=0$, so the Bernoulli polynomials of odd degree are odd functions of $\left(x-\frac12\right)$. Then $$B_{2n+1}(x)=B_{2n+1}\left(\frac12+\left(x-\frac12\right)\right)=-B_{2n+1}\left(\frac12-\left(x-\frac12\right)\right)=B_{2n+1}(1-x)$$ so $$B_{2n+1}(0)=-B_{2n+1}(1)=B_{2n+1}(1)=0$$ For $n>1$, where the middle equality follows from the boundary condition and the last from the middle one. The Bernoulli numbers may be defined as $B_n=B_n(1)$. Now we can do some integration by parts: $$\begin{align}\int_0^1f(x)dx&=\int_0^1B_0(x)f(x)dx=\frac11\left[B_1(x)f(x)\right]_0^1-\frac1{1!}\int_0^1B_1(x)f^{\prime}(x)dx\\ &=\left(1-\frac12\right)f(1)-\left(0-\frac12\right)f(0)-\frac1{1!}\left\{\left[\frac12B_2(x)f^{\prime}(x)\right]_0^1-\frac12\int_0^1B_2(x)f^{\prime\prime}(x)dx\right\}\\ &=\frac12\left[f(0)+f(1)\right]-\frac1{2!}B_2\left[f^{\prime}(1)-f^{\prime}(0)\right]+\frac1{2!}\int_0^1B_2(x)f^{\prime\prime}(x)dx\\ &=\cdots\\ &=\frac12\left[f(0)+f(1)\right]+\sum_{j=2}^{2m+2}\frac{(-1)^{j+1}}{j!}B_j\left[f^{(j-1)}(1)-f^{(j-1)}(0)\right]\\ &+\frac{(-1)^{2m+2}}{(2m+2)!}\int_0^1B_{2m+2}(x)f^{(2m+2)}(x)dx\\ &=\frac12\left[f(0)+f(1)\right]-\sum_{j=1}^{m+1}\frac1{(2j)!}B_{2j}\left[f^{(2j-1)}(1)-f^{(2j-1)}(0)\right]\\ &+\frac1{(2m+2)!}\int_0^1B_{2m+2}(x)f^{(2m+2)}(x)dx\end{align}$$ Because the odd Bernoulli numbers are zero except for $B_1=\frac12$. Now we can establish the Euler-MacLaurin summation formula: $$\begin{align}\int_0^nf(x)dx&=\sum_{\ell=0}^{n-1}\int_{\ell}^{\ell+1}f(x)dx=\sum_{\ell=0}^{n-1}\int_0^1f(x+\ell)dx\\ =&\sum_{\ell=0}^{n-1}\left\{\frac12\left[f(\ell)+f(\ell+1)\right]-\sum_{j=1}^{m+1}\frac1{(2j)!}B_{2j}\left[f^{(2j-1)}(\ell+1)-f^{(2j-1)}(\ell)\right]\right.\\ &\left.+\frac1{(2m+2)!}\int_0^1B_{2m+2}(x)f^{(2m+2)}(x+\ell)dx\right\}\\ &=\frac12f(0)+\sum_{\ell=1}^{n-1}f(\ell)+\frac12f(n)-\sum_{j=1}^{m+1}\frac1{(2j)!}B_{2j}\left[f^{(2j-1)}(n)-f^{(2j-1)}(0)\right]\\ &+\frac1{(2m+2)!}\int_0^1B_{2m+2}(x)\sum_{\ell=0}^{n-1}f^{(2m+2)}(x+\ell)dx\end{align}$$ Where we have telescoped the middle term. Then if we let $x=a+\frac{b-a}nt$, $$\begin{align}\int_a^bf(x)dx&=\left(\frac{b-a}n\right)\int_0^nf\left(a+\frac{b-a}nt\right)dt\\ &=\left(\frac{b-a}n\right)\left\{\frac12f(a)+\sum_{\ell=0}^{n-1}f\left(a+\frac{b-a}n\ell\right)+\frac12f(b)\right.\\ &\left.-\sum_{j=1}^{m+1}\frac1{(2j)!}B_{2j}\left(\frac{b-a}n\right)^{2j-1}\left[f^{2j-1}(b)-f^{2j-1}(a)\right]\right.\\ &\left.+\frac1{(2m+2)!}\left(\frac{b-a}n\right)^{2m+2}\int_0^1B_{2m+2}(t)\sum_{\ell=0}^{n-1}f^{(2m+2)}\left(a+\frac{b-a}n(t+\ell)\right)dt\right\}\\ &=\left(\frac{b-a}n\right)\left[\frac12f(a)+\sum_{\ell=0}^{n-1}f\left(a+\frac{b-a}n\ell\right)+\frac12f(b)\right]\\ &-\sum_{j=1}^{m+1}\frac1{(2j)!}B_{2j}\left(\frac{b-a}n\right)^{2j}\left[f^{2j-1}(b)-f^{2j-1}(a)\right]\\ &+\frac1{(2m+2)!}\left(\frac{b-a}n\right)^{2m+3}\int_0^1B_{2m+2}(t)\sum_{\ell=0}^{n-1}f^{(2m+2)}\left(a+\frac{b-a}n(t+\ell)\right)dt\end{align}$$ This is the trapezoidal rule with correction terms. If $f(x)$ is a polynomial of degree at most $2m+2$, then $f^{2m+2}\left(a+\frac{b-a}n(t+\ell)\right)$ is a constant and $$\int_0^1B_{2m+2}(x)dx=\frac1{2m+3}\left.B_{2m+3}(x)\right|_0^1=\frac1{2m+3}\left[B_{2m+3}(1)-B_{2m+3}(0)\right]=0$$ So the third line above drops out and we have $$\int_a^bf(x)dx=R(f,a,b,n,0)-\sum_{j=1}^{m+1}\frac1{(2j)!}B_{2j}\left(\frac{b-a}n\right)^{2j}\left[f^{2j-1}(b)-f^{2j-1}(a)\right]$$ Where $R(f,a,b,n,k)$ is the Romberg rule applied to $f(x)$ over the interval $[a,b]$ with initially $n+1$ points after $k$ refinements. Of course $R(f,a,b,n,0)$ is the trapezoidal rule and $$R(f,a,b,n,k)=\frac{2^{2k}R(f,a,b,2n,k-1)-R(f,a,b,n,k-1)}{2^{2k}-1}$$ Applying $m$ refinements, $$\int_a^bf(x)dx=R(f,a,b,n,m)-\sum_{j=1}^{m+1}\frac1{(2j)!}B_{2j}\left(\frac{b-a}n\right)^{2j}\left[f^{2j-1}(b)-f^{2j-1}(a)\right]\prod_{k=1}^m\frac{\frac{2^{2k}}{2^{2j}}-1}{2^{2k}-1}$$ Now, if $1\le j\le m$, then $\prod_{k=1}^m\left(\frac{2^{2k}}{2^{2j}}-1\right)=0$, so the only term in the sum that survives is the $j=m+1$ term and $$\begin{align}\prod_{k=1}^m\left(\frac{2^{2k}}{2^{2m+2}}-1\right)&=\prod_{k=1}^m\left(2^{2k-2m-2}-1\right)=\prod_{\ell=1}^m\left(2^{-2\ell}-1\right)\\ &=\frac{\prod_{\ell=1}^m\left(1-2^{2\ell}\right)}{\prod_{\ell=1}^m2^{2\ell}}=\frac{(-1)^m}{2^{m^2+m}}\prod_{\ell=1}^m\left(2^{2\ell}-1\right)\end{align}$$ Where we have made the substitution $\ell=m+1-k$. So if $f(x)$ is a polynomial of degree at most $2m+2$, then $$\begin{align}\int_a^bf(x)dx&=R(f,a,b,n,m)+\frac{(-1)^{m+1}B_{2m+2}}{2^{m^2+m}(2m+2)!}\left(\frac{b-a}n\right)^{2m+2}\left[f^{2m+1}(b)-f^{2m+1}(a)\right]\\ &=R(f,a,b,n,m)+\frac{(-1)^{m+1}B_{2m+2}}{2^{m^2+m}(2m+2)!}\frac{(b-a)^{2m+3}}{n^{2m+2}}f^{2m+2}(\xi)\end{align}$$ For any $a<\xi<b$. I don't have a proof that this is the case for some $a<\xi<b$ if $f(x)$ is not a polynomial of degree at most $2m+2$, however. Not to mention that I'm too tired to check my result just now...

EDIT I wrote a program that computes $\int_a^bx^{2m+2}dx$ and $R(x^{2m+2},a,b,n,m)$ then takes their difference and multiplies by $2^{m^2+m}n^{2m+2}/(b-a)^{2m+3}$ because $f^{(2m+2)}(\xi)=(2m+2)!$ in this case and it does seem to be producing $(-1)^{m+1}B_{2m+2}$ before roundoff error invalidates its results:

program romberg
   use ISO_FORTRAN_ENV
   implicit none
   real(REAL128) a, b
   integer i, j, m, n
   real(REAL128) h
   integer, parameter :: mmax = 7
   real(REAL128) s(0:mmax)
   a = 3
   b = 8
   n = 7
   do m = 0, mmax
      h = (b-a)/n
      s(0) = h*(a**(2*m+2)+2*sum([((a+(b-a)/n*j)**(2*m+2),j=1,n-1)])+b**(2*m+2))/2
      do i = 1, m
         h = h/2
         s(i) = s(i-1)/2+h*sum([((a+(2*j-1)*h)**(2*m+2),j=1,n*2**(i-1))])
         do j = i-1, 0, -1
            s(j) = (4**(i-j)*s(j+1)-s(j))/(4**(i-j)-1)
         end do
      end do
      write(*,*) m, ((b**(2*m+3)-a**(2*m+3))/(2*m+3)-s(0))* &
         real(n,REAL128)**(2*m+2)*2.0_REAL128**(m**2+m)/(b-a)**(2*m+3)
   end do
end program romberg

Output:

   0 -0.166666666666666666666666666666664
   1 -3.333333333333333333333333333536527E-0002
   2 -2.380952380952380952380953375265230E-0002
   3 -3.333333333333333333335858504793432E-0002
   4 -7.575757575757575760082802140971125E-0002
   5 -0.253113553113710066890623420476913
   6  -1.16666671258131681811523437500000
   7  -7.06133308742882048000000000000000