{
"cells": [
{
"cell_type": "markdown",
"id": "configured-briefs",
"metadata": {},
"source": [
"# The Romberg Method\n",
"\n",
"**Revised on Monday April 7**, adding some notes on presenting the results and errors nicely."
]
},
{
"cell_type": "markdown",
"id": "peaceful-shipping",
"metadata": {},
"source": [
"Write a function to implement the Romberg method.\n",
"\n",
"- Again, the input should include an error tolerance, and the output should include estimates of absolute error and of cost.\n",
"\n",
"- We will discuss this in class; there are lots of details to address.\n",
"\n",
"- Note that some ideas for this are first explored in the part of\n",
"[Centered Difference Approximation with Richardson Extrapolation](task-1-and-2-centered-difference-derivative-approximation-and-richardson.ipynb)."
]
},
{
"cell_type": "markdown",
"id": "cosmetic-correspondence",
"metadata": {},
"source": [
"**Additions on April 5**"
]
},
{
"cell_type": "markdown",
"id": "saving-vanilla",
"metadata": {},
"source": [
"Recall that the basic algorithm is\n",
"\n",
"$n \\leftarrow 1$\n",
"
\n",
"$h \\leftarrow b-a$\n",
"
\n",
"$\\displaystyle R_{0,0} = \\frac{f(a) + f(b)}{2} h$\n",
"
\n",
"`for i from 1 to M:`\n",
"
\n",
"$\\quad$ $R_{i,0} = \\left( R_{i-1,0} + h \\sum_{k=1}^n f(a + (i - 1/2)h) \\right)/2$\n",
"
\n",
"$\\quad$ `for j from 1 to i:`\n",
"
\n",
"$\\quad\\quad \\displaystyle R_{i,j} = \\frac{4^j R_{i,j-1} - R_{i-1,j-1}}{4^j - 1}$\n",
"
\n",
"$\\quad$ `end for`\n",
"
\n",
"$\\quad n \\leftarrow 2 n$\n",
"
\n",
"$\\quad h \\leftarrow h/2$\n",
"
\n",
"`end for`"
]
},
{
"cell_type": "markdown",
"id": "invisible-consequence",
"metadata": {},
"source": [
"**Python note:**\n",
"The needed 2D array for $R$ can be created with the numpy function `R = numpy.array([M,M])`.\n",
"It has the syntactic peculiarity that the dimensions of the array are given in a list, not as successive input parameters."
]
},
{
"cell_type": "markdown",
"id": "welsh-juvenile",
"metadata": {},
"source": [
"Once the above basic verion is working, produce a second version that adds an error estimate,\n",
"and ends once the final approximation $R_{i,i}$ at stage $i$ is good enough.\n",
"\n",
"As usual I recommend using a `for` loop organization, adding an `if ... break`"
]
},
{
"cell_type": "markdown",
"id": "about-temple",
"metadata": {},
"source": [
"**Additions on April 7:**\n",
"some tips on printing, and formatting numbers.\n",
"\n",
"The final array only has relevant value at positions `R[i,j]` with $ j \\leq i$, so here is one way to print it:"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "transsexual-pennsylvania",
"metadata": {},
"outputs": [],
"source": [
"import numpy as np"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "nuclear-practice",
"metadata": {},
"outputs": [],
"source": [
"M = 4\n",
"R = np.zeros([M,M])\n",
"# Set up a fake array R;\n",
"# the entries are from an important example: the Hilbert matrix. We will see it again soon.\n",
"for i in range(M):\n",
" for j in range(i+1):\n",
" R[i, j] = 1.0/(1+i+j)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "korean-inspector",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"R is:\n",
"[[1. 0. 0. 0. ]\n",
" [0.5 0.33333333 0. 0. ]\n",
" [0.33333333 0.25 0.2 0. ]\n",
" [0.25 0.2 0.16666667 0.14285714]]\n",
"The relevant parts of R are:\n",
"[1.]\n",
"[0.5 0.33333333]\n",
"[0.33333333 0.25 0.2 ]\n",
"[0.25 0.2 0.16666667 0.14285714]\n"
]
}
],
"source": [
"print(f\"R is:\")\n",
"print(R)\n",
"print(f\"The relevant parts of R are:\")\n",
"for i in range(M):\n",
" print(R[i,0:i+1])"
]
},
{
"cell_type": "markdown",
"id": "binary-senior",
"metadata": {},
"source": [
"This strategy can also be used to print the errors.\n",
"One can subtract a number from an array, which subtracts it from each element if the array.\n",
"So pretending that the exact value is 1/8:"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "immune-freeware",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The fake errors are:\n",
"[0.875]\n",
"[0.375 0.20833333]\n",
"[0.20833333 0.125 0.075 ]\n",
"[0.125 0.075 0.04166667 0.01785714]\n"
]
}
],
"source": [
"exact = 1.0/8.0\n",
"print(f\"The fake errors are:\")\n",
"for i in range(M):\n",
" print(R[i,0:i+1] - exact)"
]
},
{
"cell_type": "markdown",
"id": "internal-likelihood",
"metadata": {},
"source": [
"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:"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "grave-mauritius",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The fake errors are:\n",
"0.875 \n",
"0.375 0.20833333333333331 \n",
"0.20833333333333331 0.125 0.07500000000000001 \n",
"0.125 0.07500000000000001 0.04166666666666666 0.01785714285714285 \n"
]
}
],
"source": [
"exact = 1.0/8.0\n",
"print(f\"The fake errors are:\")\n",
"for i in range(M):\n",
" for j in range(i+1):\n",
" print(R[i,j] - exact, end=\" \") # This puts the \"end\" text (a space) at the end, instead of a new line.\n",
" print() # Start a new line only at the end of each row"
]
},
{
"cell_type": "markdown",
"id": "formed-vacation",
"metadata": {},
"source": [
"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:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "foreign-craps",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The fake errors are:\n",
" 0.875 \n",
" 0.375 0.2083 \n",
" 0.2083 0.125 0.075 \n",
" 0.125 0.075 0.04167 0.01786 \n"
]
}
],
"source": [
"exact = 1.0/8.0\n",
"print(f\"The fake errors are:\")\n",
"for i in range(M):\n",
" for j in range(i+1):\n",
" print(f\"{R[i,j] - exact:7.4}\", end=\" \") # This puts the \"end\" text (a space) at the end, instead of a new line.\n",
" print() # Start a new line only at the end of each row"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.5"
}
},
"nbformat": 4,
"nbformat_minor": 5
}