The Romberg Method¶
Revised on Monday April 7, adding some notes on presenting the results and errors nicely.
Write a function to implement the Romberg method.
Again, the input should include an error tolerance, and the output should include estimates of absolute error and of cost.
We will discuss this in class; there are lots of details to address.
Note that some ideas for this are first explored in the part of Centered Difference Approximation with Richardson Extrapolation.
Additions on April 5
Recall that the basic algorithm is
\(n \leftarrow 1\)
\(h \leftarrow b-a\)
\(\displaystyle R_{0,0} = \frac{f(a) + f(b)}{2} h\)
for i from 1 to M:
\(\quad\) \(R_{i,0} = \left( R_{i-1,0} + h \sum_{k=1}^n f(a + (i - 1/2)h) \right)/2\)
\(\quad\) for j from 1 to i:
\(\quad\quad \displaystyle R_{i,j} = \frac{4^j R_{i,j-1} - R_{i-1,j-1}}{4^j - 1}\)
\(\quad\) end for
\(\quad n \leftarrow 2 n\)
\(\quad h \leftarrow h/2\)
end for
Python note:
The needed 2D array for \(R\) can be created with the numpy function R = numpy.array([M,M]).
It has the syntactic peculiarity that the dimensions of the array are given in a list, not as successive input parameters.
Once the above basic verion is working, produce a second version that adds an error estimate, and ends once the final approximation \(R_{i,i}\) at stage \(i\) is good enough.
As usual I recommend using a for loop organization, adding an if ... break
Additions on April 7: some tips on printing, and formatting numbers.
The final array only has relevant value at positions R[i,j] with \( j \leq i\), so here is one way to print it:
import numpy as np
M = 4
R = np.zeros([M,M])
# Set up a fake array R;
# the entries are from an important example: the Hilbert matrix. We will see it again soon.
for i in range(M):
for j in range(i+1):
R[i, j] = 1.0/(1+i+j)
print(f"R is:")
print(R)
print(f"The relevant parts of R are:")
for i in range(M):
print(R[i,0:i+1])
R is:
[[1. 0. 0. 0. ]
[0.5 0.33333333 0. 0. ]
[0.33333333 0.25 0.2 0. ]
[0.25 0.2 0.16666667 0.14285714]]
The relevant parts of R are:
[1.]
[0.5 0.33333333]
[0.33333333 0.25 0.2 ]
[0.25 0.2 0.16666667 0.14285714]
This strategy can also be used to print the errors. One can subtract a number from an array, which subtracts it from each element if the array. So pretending that the exact value is 1/8:
exact = 1.0/8.0
print(f"The fake errors are:")
for i in range(M):
print(R[i,0:i+1] - exact)
The fake errors are:
[0.875]
[0.375 0.20833333]
[0.20833333 0.125 0.075 ]
[0.125 0.075 0.04166667 0.01785714]
An alternative is to print one entry at a time, telling the print command to not start a new line except at the end of a row:
exact = 1.0/8.0
print(f"The fake errors are:")
for i in range(M):
for j in range(i+1):
print(R[i,j] - exact, end=" ") # This puts the "end" text (a space) at the end, instead of a new line.
print() # Start a new line only at the end of each row
The fake errors are:
0.875
0.375 0.20833333333333331
0.20833333333333331 0.125 0.07500000000000001
0.125 0.07500000000000001 0.04166666666666666 0.01785714285714285
One advantage of this “one at a time” approach is that one can specify the number of significant digits to display, using f-string formatting. Here, the “:7.4” part at the end of the braces means to use y spaces for each value (so they line up in columns) wth 4 sigificant digits:
exact = 1.0/8.0
print(f"The fake errors are:")
for i in range(M):
for j in range(i+1):
print(f"{R[i,j] - exact:7.4}", end=" ") # This puts the "end" text (a space) at the end, instead of a new line.
print() # Start a new line only at the end of each row
The fake errors are:
0.875
0.375 0.2083
0.2083 0.125 0.075
0.125 0.075 0.04167 0.01786