{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"\n",
"# Polynomial interpolation: Newton interpolation\n",
"\n",
" \n",
"**Anne Kværnø (modified by André Massing)**\n",
"\n",
"Date: **Jan 18, 2021**"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"If you want to have a nicer theme for your jupyter notebook,\n",
"download the [cascade stylesheet file tma4125.css](https://www.math.ntnu.no/emner/TMA4125/2020v/part_II/notebooks/tma4125.css)\n",
"and execute the next cell:"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"text/html": [
" \n",
"\n"
],
"text/plain": [
""
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"from IPython.core.display import HTML\n",
"from numpy import pi\n",
"from numpy.linalg import norm, solve # Solve linear systems and compute norms\n",
"\n",
"\n",
"def css_styling():\n",
" try:\n",
" with open(\"tma4125.css\") as f:\n",
" styles = f.read()\n",
" return HTML(styles)\n",
" except FileNotFoundError:\n",
" pass # Do nothing\n",
"\n",
"\n",
"# Comment out next line and execute this cell to restore the default notebook style\n",
"css_styling()"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"And of course we want to import the required modules."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [],
"source": [
"%matplotlib inline\n",
"\n",
"\n",
"newparams = {\n",
" \"figure.figsize\": (6.0, 6.0),\n",
" \"axes.grid\": True,\n",
" \"lines.markersize\": 8,\n",
" \"lines.linewidth\": 2,\n",
" \"font.size\": 14,\n",
"}\n",
"plt.rcParams.update(newparams)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Newton interpolation"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"This is an alternative approach to find the interpolation polynomial.\n",
"Let $x_0,x_1,\\ldots,x_n$ be $n+1$ distinct real numbers.\n",
"Instead of using the Lagrange polynomials to write the\n",
"interpolation polynomial in Lagrange form, we will now employ\n",
"the **Newton** polynomials $\\omega_i$, $i=0,\\ldots, n$.\n",
"The Newton polynomials are defined by\n",
"as follows:"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"$$\n",
"\\begin{align*}\n",
"\\omega_0(x) &= 1,\n",
"\\\\ \n",
"\\omega_1(x) &= (x-x_0),\n",
"\\\\ \n",
"\\omega_2(x) &= (x-x_0)(x-x_1),\n",
"\\\\ \n",
"\\ldots\n",
"\\\\ \n",
"\\omega_n(x) &= (x-x_0)(x-x_1)\\cdots(x-x_{n-1}),\n",
"\\end{align*}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"or in more compact notation\n",
"\n",
"\n",
"\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\omega_i(x)= \\prod_{k=0}^{i-1}(x-x_k).\n",
"\\label{eq:newton_poly} \\tag{1}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"The so-called\n",
"**Newton form** of a polynomial of degree $n$ is an expansion of the form\n",
"\n",
"$$\n",
"p(x)=\\sum_{i=0}^{n} c_i \\omega_i(x)\n",
"$$\n",
"\n",
"or more explicitly\n",
"\n",
"$$\n",
"p(x)=c_n (x-x_0)(x-x_1)\\cdots(x-x_{n-1}) + c_{n-1}(x-x_0)(x-x_1)\\cdots(x-x_{n-2}) + \\cdots + c_1(x-x_0) + c_0.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"In the light of this form of writing a polynomial, the polynomial\n",
"interpolation problem leads to the following observations. Let us\n",
"start with a single node $x_0$, then $f(x_0)=p(x_0)=c_0$. "
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"Going one\n",
"step further and consider two nodes $x_0,x_1$. Then we see that\n",
"$f(x_0)=p(x_0)=c_0$ and $f(x_1)=p(x_1)=c_0 + c_1(x_1-x_0)$. The latter\n",
"implies that the coefficient"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"$$\n",
"c_1=\\frac{f(x_1)-f(x_0)}{x_1-x_0}.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"Given three nodes $x_0,x_1,x_2$ yields the coefficients $c_0,c_1$ as defined above, and from\n",
"\n",
"$$\n",
"f(x_2)=p(x_2)=c_0 + c_1(x_2-x_0) + c_2(x_2-x_0)(x_2-x_1)\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"we deduce the coefficient\n",
"\n",
"$$\n",
"c_2=\\frac{f(x_2) - f(x_0) - \\frac{f(x_1)-f(x_0)}{x_1-x_0} (x_2-x_0)}{(x_2-x_0)(x_2-x_1)}.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"Playing with this quotient gives the much more structured expression\n",
"\n",
"$$\n",
"c_2=\\frac{\\frac{f(x_2)-f(x_1)}{x_2-x_1} - \\frac{f(x_1)-f(x_0)}{x_1-x_0}}{(x_2-x_0)}.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"This procedure can be continued and yields a so-called triangular\n",
"systems that permits to define the remaining coefficients\n",
"$c_3,\\ldots,c_n$. One sees quickly that the coefficient $c_k$ only\n",
"depends on the interpolation points $(x_0,y_0),\\ldots,(x_k,y_k)$,\n",
"where $y_i:=f(x_i)$, $i=0,\\ldots,n$."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"We introduce the folllwing so-called **finite difference** notation for a\n",
"function $f$. The 0th order finite difference is defined to be\n",
"$f[x_0]:=f(x_0)$. The 1st order finite difference is\n",
"\n",
"$$\n",
"f[x_0,x_1]:=\\frac{f(x_1)-f(x_0)}{x_1-x_0}.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"The second order finite difference is defined by\n",
"\n",
"$$\n",
"f[x_0,x_1,x_2]:= \\frac{f[x_1,x_2] - f[x_0,x_1]}{x_2-x_0}.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"In general, the **nth order finite difference** of the function $f$,\n",
"also called the **nth Newton divided difference**, is defined recursively by\n",
"\n",
"$$\n",
"f[x_0,\\ldots,x_n]:= \\frac{f[x_1,\\ldots,x_n] - f[x_0,\\ldots,x_{n-1}]}{x_n-x_0}.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"Newton's method to solve the polynomial interpolation problem can be\n",
"summarized as follows. Given $n+1$ interpolation points\n",
"$(x_0,y_0),\\ldots,(x_n,y_n)$, $y_i:=f(x_i)$. If the order $n$\n",
"interpolation polynomial is expressed in Newton's form\n",
"\n",
"$$\n",
"p_n(x)=c_n (x-x_0)(x-x_1)\\cdots(x-x_{n-1}) + c_{n-1}(x-x_0)(x-x_1)\\cdots(x-x_{n-2}) + \\cdots + c_1(x-x_0) + c_0,\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"then the coefficients\n",
"\n",
"$$\n",
"c_k = f[x_0,\\ldots,x_k]\n",
"$$\n",
"\n",
"for $k=0,\\ldots,n$. "
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"In fact, a recursion is in place\n",
"\n",
"$$\n",
"p_n(x)=p_{n-1}(x) + f[x_0,\\ldots,x_n](x-x_0)(x-x_1)\\cdots(x-x_{n-1})\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"It is common to write the finite differences in a table, which for $n=3$ will\n",
"look like:\n",
"\n",
"$$\n",
"\\begin{array}{c|cccc}\n",
"x_0 & f[x_0] & \\\\ \n",
" & & f[x_0,x_1] & \\\\ \n",
"x_1 & f[x_1] & & f[x_0,x_1,x_2] \\\\ \n",
" & & f[x_1,x_2] & & f[x_0,x_1,x_2, x_3] \\\\ \n",
"x_2 & f[x_2] & & f[x_1,x_2,x_3] \\\\ \n",
" & & f[x_2,x_3] & \\\\ \n",
"x_3 & f[x_3] \\\\ \n",
"\\end{array}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"**Example 1 again:**\n",
"Given the points in Example 1. The corresponding table of divided differences\n",
"becomes:\n",
"\n",
"$$\n",
"\\begin{array}{c|cccc}\n",
"0 & 1 & \\\\ \n",
" & & -3/4 & \\\\ \n",
"2/3 & 1/2 & & -3/4 \\\\ \n",
" & & -3/2 & \\\\ \n",
"1 & 0 & \n",
"\\end{array}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"and the interpolation polynomial becomes\n",
"\n",
"$$\n",
"p_2(x) = 1 - \\frac{3}{4}(x-0)-\\frac{3}{4}(x-0)(x-\\frac23) = 1 - \\frac{1}{4}x -\n",
"\\frac{3}{4} x^2.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Implementation\n",
"The method above is implemented as two functions:\n",
"* `divdiff(xdata, ydata)`: Create the table of divided differences\n",
"\n",
"* `newtonInterpolation(F, xdata, x)`: Evaluate the interpolation polynomial.\n",
"\n",
"Here, `xdata` and `ydata` are arrays with the interpolation points, and `x` is an \n",
"array of values in which the polynomial is evaluated."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [],
"source": [
"def divdiff(xdata, ydata):\n",
" # Create the table of divided differences based\n",
" # on the data in the arrays x_data and y_data.\n",
" n = len(xdata)\n",
" F = np.zeros((n, n))\n",
" F[:, 0] = ydata # Array for the divided differences\n",
" for j in range(n):\n",
" for i in range(n - j - 1):\n",
" F[i, j + 1] = (F[i + 1, j] - F[i, j]) / (xdata[i + j + 1] - xdata[i])\n",
" return F # Return all of F for inspection.\n",
" # Only the first row is necessary for the\n",
" # polynomial.\n",
"\n",
"\n",
"def newton_interpolation(F, xdata, x):\n",
" # The Newton interpolation polynomial evaluated in x.\n",
" n, m = np.shape(F)\n",
" xpoly = np.ones(len(x)) # (x-x[0])(x-x[1])...\n",
" newton_poly = F[0, 0] * np.ones(len(x)) # The Newton polynomial\n",
" for j in range(n - 1):\n",
" xpoly = xpoly * (x - xdata[j])\n",
" newton_poly = newton_poly + F[0, j + 1] * xpoly\n",
" return newton_poly"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"Run the code on the example above:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [],
"source": [
"# Example: Use of divided differences and the Newton interpolation\n",
"# formula.\n",
"xdata = [0, 2 / 3, 1]\n",
"ydata = [1, 1 / 2, 0]\n",
"F = divdiff(xdata, ydata) # The table of divided differences\n",
"print(\"The table of divided differences:\\n\", F)\n",
"\n",
"x = np.linspace(0, 1, 101) # The x-values in which the polynomial is evaluated\n",
"p = newton_interpolation(F, xdata, x)\n",
"plt.plot(x, p) # Plot the polynomial\n",
"plt.plot(xdata, ydata, \"o\") # Plot the interpolation points\n",
"plt.title(\"The interpolation polynomial p(x)\")\n",
"plt.grid(True)\n",
"plt.xlabel(\"x\");"
]
}
],
"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.9.1"
}
},
"nbformat": 4,
"nbformat_minor": 4
}