Romberg from class

The following does Romberg integration one step at a time. The next section gives code to automate this.

f[x_] = Sin[x];
a = 1. ;
b = 3. ;
n = 2;
h = (b - a)/n;
approx = 
      h/2*(f[a] + 2  Sum[f[a + k*h], {k, 1, n - 1, 1}] + 
            f[b])

exact = -Cos[3. ] + Cos[1. ]

exact - approx

trap1 = 
      Table[n = 
          2^m; 
h = (b - a)/n; 
h/
            2*(f[a] + 2  Sum[f[a + k*h], {k, 1, n - 1, 1}] + f[b]), {m, 
          1, 4}]

exact - trap1[[-1]]

trap2 = 
      Table[trap1[[k]] + (trap1[[k]] - trap1[[k - 1]])/3, {k, 
          2, Length[trap1]}]

exact - trap2[[-1]]

trap3 = 
      Table[trap2[[k]] + (trap2[[k]] - trap2[[k - 1]])/
            15, {k, 2, Length[trap2]}]

exact - trap3[[-1]]

trap4 = 
      Table[trap3[[k]] + (trap3[[k]] - 
                trap3[[k - 1]])/(4^3 - 1), {k, 2, 
          Length[trap3]}]

exact - trap4[[-1]]

Automated version of Romberg

f[x_] = Sin[x];
a = 1. ;
b = 3;
n = 1;
levels = 5;
h = (b - a)/n;
romberg = {};
romberg = 
        Append[romberg, 
          Table[n = 
              2^m; 
h = (b - a)/n; 
h/
                2*(f[a] + 2  Sum[f[a + k*h], {k, 1, n - 1, 1}] + 
                  f[b]), {m, 0, levels}]];
Do[
        romberg = 
          Append[romberg, 
Table[
              romberg[[m, 
                    k]] + (romberg[[m, k]] - 
                      romberg[[m, k - 1]])/(4^m - 1), {k, 2, 
                Length[romberg[[m]]]}]], 
{m, 1, 
          levels}];
romberg

TableForm[romberg]

Note that the table above interchanges rows and columns from what we usually write down by hand.  That is, the first row here gives the initial approximations, the second row gives the first step of the extrapolation, and so on.

exact = -Cos[3. ] + Cos[1. ]

approx = romberg[[-1, -1]]

approx - exact