{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Polynomial Collocation (Interpolation/Extrapolation)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**References:**\n",
"\n",
"- Section 3.1 *Data and Interpolating Functions* in {cite}`Sauer`.\n",
"- Section 3.1 *Interpolation and the Lagrange Polynomial* in {cite}`Burden-Faires`.\n",
"- Section 4.1 in {cite}`Chenney-Kincaid`."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Introduction\n",
"\n",
"Numerical methods for dealing with functions require describing them, at least approximately, using a finite list of numbers, and the most basic approach is to approximate by a polynomial.\n",
"(Other important choices are rational functions and \"trigonometric polynomials\": sums of multiples of sines and cosines.)\n",
"Such polynomials can then be used to aproximate derivatives and integrals.\n",
"\n",
"The simplest idea for approximating $f(x)$ on domain $[a, b]$ is to start with a finite collection of **node** values $x_i \\in [a, b]$, $0 \\leq i \\leq n$ and then seek a polynomial $p$ which **collocates** with $f$ at those values: $p(x_i) = f(x_i)$ for $0 \\leq i \\leq n$.\n",
"Actually, we can put the function aside for now, and simply seek a polynomial that passes through a list of points $(x_i, y_i)$; later we will achieve collocation with $f$ by choosing $y_i = f(x_i)$.\n",
"\n",
"In fact there are infinitely many such polynomials: given one, add to it any polynomial with zeros at all of the $n+1$ notes.\n",
"So to make the problem well-posed, we seek the collocating polynomial of *lowest degree*."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```{prf:theorem} Existence and uniqueness of a collocating polynomial\n",
":label: theorem-collocation\n",
"\n",
"Given $n+1$ distinct values $x_i$, $0 \\leq i \\leq n$, and corresponding $y$-values $y_i$,\n",
"there is a unique polynomial $P$ of degree at most $n$ with $P(x_i) = y_i$ for $0 \\leq i \\leq n$.\n",
"\n",
"(*Note:* although the degree is typically $n$, it can be less; as an extreme example, if all $y_i$ are equal to $c$, then $P(x)$ is that constant $c$.)\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Historically there are several methods for finding $P_n$ and proving its uniqueness, in particular, the *divided difference* method introduced by Newton and the *Lagrange polynomial* method.\n",
"However for our purposes, and for most modern needs, a different method is easiest, and it also introduces a strategy that will be of repeated use later in this course: the **Method of Undertermined Coefficients** or **MUC**.\n",
"\n",
"In general, this method starts by assuming that the function wanted is a sum of unknown multiples of a collection of known functions.\n",
"Here, $P(x) = c_n x^n + c_{n-1} x^{n-1} + \\cdots+ c_1 x + c_0 = \\sum_{j=0}^n c_j x^j$.\n",
"
\n",
"(*Note:* any of the $c_i$ could be zero, including $c_n$, in which case the degree is less than $n$.)\n",
"
\n",
"The unknown factors ($c_0 \\cdots c_n$) are the *undetermined coefficients.*\n",
"\n",
"Next one states the problem as a system of equations for these undetermined coefficients, and solves them.\n",
"
\n",
"Here, we have $n+1$ conditions to be met:\n",
"\n",
"$$P(x_i) = \\sum_{j=0}^n c_j x_i^j = y_i, \\quad 0 \\leq i \\leq n$$\n",
"\n",
"This is a system if $n+1$ simultaneous linear equations in $n+1$ iunknowns, so the question of existence and uniqueness is exactly the question of whether the corresponding matrix is singular,\n",
"and so is equivalent to the case of all $y_i = 0$ having *only* the solution with all\n",
"$c_i = 0$.\n",
"\n",
"Back in terms of polynomials, this is the claim that the only polynomial of degree at most $n$ with zeros $x_0 \\dots x_n$.\n",
"And this is true, because any non-trivial polynomial with those $n+1$ distinct roots is of degree at least $n+1$, so the only \"degree n\" polynomial fitting this data is $P(x) \\equiv 0$.\n",
"The theorem is proven."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The proof of this theorem is completely constructive; it gives the only numerical method we need, and which is the one implemented in Numpy through the pair of functions `numpy.polyfit` and `numpy.polyval`.\n",
"(Aside: here as in many places, Numpy mimics the names and functionality of corresponding Matlab tools.)\n",
"\n",
"Briefly, the algorithm is this (indexing from 0 !)\n",
"- Create the $n+1 \\times n+1$ matrix $V$ with elements\n",
"\n",
"$$v_{i,j} = x_i^j,\\; 0 \\leq i \\leq n, \\, 0 \\leq j \\leq n$$\n",
"\n",
"and the $n+1$-element column vector $y$ with elements $y_i$ as above.\n",
"- Solve $V c = y$ for the vector of coefficients $c_j$ as above.\n",
"\n",
"I use the name $V$ because this is called the **Vandermonde Matrix.**"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"from numpy import array, linspace, zeros, zeros_like, exp\n",
"from matplotlib.pyplot import figure, plot, title, grid, legend"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```{prf:example}\n",
":label: interpolation-example-1\n",
"\n",
"As usual, I concoct a first example with known correct answer, by using a polynomial as $f$:\n",
"\n",
"$$ f(x) = 4 + 7x -2x^2 - 5x^3 + 2x^4 $$\n",
"\n",
"using the nodes $x_0 = 1$, $x_1 = 2$, $x_2= 0$, $x_3 = 3.3$ and $x_4 = 4$ (They do not need to be in order.)\n",
"```"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"def f(x):\n",
" return 4 + 7*x - 2*x**2 -5*x**3 + 3*x**4"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The x nodes 'x_i' are [1. 2. 7. 5. 4.]\n",
"The y values at the nodes are [ 7. 18. 5443. 1239. 448.]\n"
]
}
],
"source": [
"xnodes = array([1., 2., 7., 5., 4.]) # They do not need to be in order\n",
"nnodes = len(xnodes)\n",
"n = nnodes-1\n",
"print(f\"The x nodes 'x_i' are {xnodes}\")\n",
"ynodes = zeros_like(xnodes)\n",
"for i in range(nnodes):\n",
" ynodes[i] = f(xnodes[i])\n",
"print(f\"The y values at the nodes are {ynodes}\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The Vandermonde matrix:"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"V = zeros([nnodes, nnodes])\n",
"for i in range(nnodes):\n",
" for j in range(nnodes):\n",
" V[i,j] = xnodes[i]**j"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Solve, using our functions seen in earlier sections and gathered in {doc}`numericalMethods`"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"ename": "ImportError",
"evalue": "cannot import name 'rowReduce' from 'numerical_methods' (/Users/brenton/Library/CloudStorage/OneDrive-CollegeofCharleston/numerical-methods-and-analysis/python-edition/main/numerical_methods.py)",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mImportError\u001b[0m Traceback (most recent call last)",
"Cell \u001b[0;32mIn[7], line 1\u001b[0m\n\u001b[0;32m----> 1\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mnumerical_methods\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m rowReduce, backwardSubstitution\n",
"\u001b[0;31mImportError\u001b[0m: cannot import name 'rowReduce' from 'numerical_methods' (/Users/brenton/Library/CloudStorage/OneDrive-CollegeofCharleston/numerical-methods-and-analysis/python-edition/main/numerical_methods.py)"
]
}
],
"source": [
"from numericalMethods import rowReduce, backwardSubstitution"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"ename": "NameError",
"evalue": "name 'rowReduce' is not defined",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)",
"Cell \u001b[0;32mIn[8], line 1\u001b[0m\n\u001b[0;32m----> 1\u001b[0m (U, z) \u001b[38;5;241m=\u001b[39m \u001b[43mrowReduce\u001b[49m(V, ynodes)\n\u001b[1;32m 2\u001b[0m c \u001b[38;5;241m=\u001b[39m backwardSubstitution(U, z)\n\u001b[1;32m 3\u001b[0m \u001b[38;5;28mprint\u001b[39m(\u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mThe coefficients of P are \u001b[39m\u001b[38;5;132;01m{\u001b[39;00mc\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m\"\u001b[39m)\n",
"\u001b[0;31mNameError\u001b[0m: name 'rowReduce' is not defined"
]
}
],
"source": [
"(U, z) = rowReduce(V, ynodes)\n",
"c = backwardSubstitution(U, z)\n",
"print(f\"The coefficients of P are {c}\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can also check the resulting values of the polynomial:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"P = c[0] + c[1]*xnodes + c[2]*xnodes**2 + c[3]*xnodes**3 + c[4]*xnodes**4"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"for (x, y, P_i) in zip(xnodes, ynodes, P):\n",
" print(f\"P({x}) should be {y}; it is {P_i}\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Functions for computing the coefficients and evaluating the polynomials\n",
"\n",
"We will use this procedure several times, so it time to put it into a functions —\n",
"and add a pretty printer for polynomials."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def fitpolynomial(x, y):\n",
" \"\"\"Compute the coefficients c_i of the polynomial of lowest degree that collocates the points (x[i], y[i]).\n",
" These are returned in an array c of the same length as x and y, even if the degree is less than the normal length(x)-1,\n",
" in which case the array has some trailing zeroes.\n",
" The polynomial is thus p(x) = c[0] + c[1]x + ... c[d] x^d where n =length(x)-1, the nominal degree.\n",
" \"\"\"\n",
" nnodes = len(x)\n",
" n = nnodes - 1\n",
" V = zeros([nnodes, nnodes])\n",
" for i in range(nnodes):\n",
" for j in range(nnodes):\n",
" V[i,j] = x[i]**j\n",
" (U, z) = rowReduce(V, y)\n",
" c = backwardSubstitution(U, z)\n",
" return c"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def evaluatepolynomial(x, coeffs):\n",
" # Evaluate the polynomial with coefficients in coeffs at the points in x.\n",
" npoints = len(x)\n",
" ncoeffs = len(coeffs)\n",
" n = ncoeffs - 1\n",
" powers = linspace(0, n, n+1)\n",
" y = zeros_like(x)\n",
" for i in range(npoints):\n",
" y[i] = sum(coeffs * x[i]**powers)\n",
" return y"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def showpolynomial(c):\n",
" print(\"P(x) = \", end=\"\")\n",
" n = len(c)-1\n",
" print(f\"{c[0]:.4}\", end=\"\")\n",
" if n > 0:\n",
" coeff = c[1]\n",
" if coeff > 0:\n",
" print(f\" + {coeff:.4}x\", end=\"\")\n",
" elif coeff < 0:\n",
" print(f\" - {-coeff:.4}x\", end=\"\")\n",
" if n > 1:\n",
" for j in range(2, len(c)):\n",
" coeff = c[j]\n",
" if coeff > 0:\n",
" print(f\" + {coeff:.4}x^{j}\", end=\"\")\n",
" elif coeff < 0:\n",
" print(f\" - {-coeff:.4}x^{j}\", end=\"\")\n",
" print()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"While debugging, redo the first example:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"c_new = fitpolynomial(xnodes, ynodes)\n",
"print(c_new)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"P_i_new = evaluatepolynomial(xnodes, c_new)\n",
"\n",
"print(P_i_new)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print(f\"The values of P(x_i) are {P_i_new}\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"showpolynomial(c_new)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```{prf:example} $f(x)$ not a polynomial of degree $\\leq n$\n",
":label: interpolation-example-2\n",
"\n",
"Make an exact fit impossible by using the same function but using only four nodes and reducing the degree of the interpolating $P$ to three:\n",
"$x_0 = 1$, $x_1 = 2$, $x_2 = 3$ and $x_3 = 4$\n",
"```"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Reduce the degree of $P$ to at most 3:\n",
"n = 3\n",
"xnodes = array([-2.0, 0., 1.0, 2.])\n",
"ynodes = zeros_like(xnodes)\n",
"for i in range(len(xnodes)):\n",
" ynodes[i] = f(xnodes[i])\n",
"print(f\"n is now {n}, the nodes are now {xnodes}, with f(x_i) values {ynodes}\")\n",
"c = fitpolynomial(xnodes, ynodes)\n",
"print(f\"The coefficients of P are now {c}\")\n",
"showpolynomial(c)\n",
"print(f\"The values of P at the nodes are now {evaluatepolynomial(xnodes, c)}\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"There are several ways to assess the accuracy of this fit; we start graphically, and later consider the maximum and root-mean-square (RMS) errors."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"x = linspace(min(xnodes), max(xnodes)) # defaulting to 50 points, for graphing"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"figure(figsize=[12,6])\n",
"plot(x, f(x), label=\"y=f(x)\")\n",
"plot(xnodes, ynodes, \"*\", label=\"nodes\")\n",
"P_n_x = evaluatepolynomial(x, c)\n",
"plot(x, P_n_x, label=\"y = P_n(x)\")\n",
"legend()\n",
"grid(True);"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"P_error = f(x) - P_n_x\n",
"figure(figsize=[12,6])\n",
"title(\"Error in P_n(x)\")\n",
"plot(x, P_error, label=\"y=f(x)\")\n",
"plot(xnodes, zeros_like(xnodes), \"*\")\n",
"grid(True);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```{prf:example} $f(x)$ not a polynomial at all\n",
":label: interpolation-example-3\n",
"\n",
"$f(x) = e^x$ with five nodes, equally spaced from $-1$ to $1$\n",
"```"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def g(x): return exp(x)\n",
"a_g = -1.0\n",
"b_g = 1.0"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"n = 3\n",
"xnodes_g = linspace(a_g, b_g, n+1)\n",
"ynodes_g = zeros_like(xnodes_g)\n",
"for i in range(len(xnodes_g)):\n",
" ynodes_g[i] = g(xnodes_g[i])\n",
"print(f\"{n=}\")\n",
"print(f\"node x values {xnodes_g}\")\n",
"print(f\"node y values {ynodes_g}\")\n",
"c_g = fitpolynomial(xnodes_g, ynodes_g)\n",
"print(f\"The coefficients of P are {c_g}\")\n",
"showpolynomial(c_g)\n",
"P_values = evaluatepolynomial(c_g, xnodes_g)\n",
"print(f\"The values of P(x_i) are {P_values}\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"There are several ways to assess the accuracy of this fit.\n",
"We start graphically, and later consider the maximum and root-mean-square (RMS) errors."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"x_g = linspace(a_g - 0.25, b_g + 0.25) # Go a bit beyond the nodes in each direction"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"figure(figsize=[14,10])\n",
"title(\"With $g(x) = e^x$\")\n",
"plot(x_g, g(x_g), label=\"y = $g(x)$\")\n",
"plot(xnodes_g, ynodes_g, \"*\", label=\"nodes\")\n",
"P_g = evaluatepolynomial(x_g, c_g)\n",
"plot(x_g, P_g, label=f\"y = $P_{n}(x)$\")\n",
"legend()\n",
"grid(True);"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"P_error = g(x_g) - P_g\n",
"figure(figsize=[14,10])\n",
"title(f\"Error in $P_{n}(x)$ for $g(x) = e^x$\")\n",
"plot(x_g, P_error)\n",
"plot(xnodes_g, zeros_like(xnodes_g), \"*\")\n",
"grid(True);"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"P_error = f(x) - P_m_x\n",
"figure(figsize=[14,10])\n",
"title(\"Error in $P_1(x)$ for $f(x)$\")\n",
"plot(x, P_error, label=\"y=f(x)\")\n",
"plot(xnodes, zeros_like(xnodes), \"*\")\n",
"grid(True);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"---\n",
"This work is licensed under [Creative Commons Attribution-ShareAlike 4.0 International](https://creativecommons.org/licenses/by-sa/4.0/)"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"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.9.16"
}
},
"nbformat": 4,
"nbformat_minor": 4
}