Suppose you have a list of \(n+1\) given point \((x_i, y_i)\) with \(i \in \{0, \dots, n\}\) and \(\forall i,j \in \{0, \dots, n\}: i \neq j \Rightarrow x_i \neq x_j\).
Now you want to find a polynomial \(\displaystyle p(x) = \sum_{i=0}^n a_i \cdot x^i\) that goes through all of the given points. This means:
\(\forall i \in \{0, \dots, n\}: p(x_i) = y_i\)
Existence and uniqueness
Theorem: When you have a list of \(n+1\) points with mutually different \(x_i\) there is exactly one polynom of degree \(\leq n\) that fits those points.
Proof:
You can formulate a linear system of equations:
A matrix of the form of \(A\) is called Vandermonde matrix. When the data points \(x_i\) are mutually different, it is known that the Vandermonde matrix is invertible (source).
This means, the solution \((a_0 \dots a_n)^T\) of this linear equation gives the polynom \(p(x) = \sum_{i=0}^n a_i x^i\)
So the solution exists and is unique \(\blacksquare\)
Straight forward interpolating polynomials
For this algorithm, I'll find the polynomial in its monomial from \(p(x) = \sum_{i=0}^n a_i x^i\). I'll use the matrix \(A\) from section "Uniqueness".
You might want to take a look at my article about Gaussian elimination.
#!/usr/bin/env python
# -*- coding: utf-8 -*-
def pprintGaus(A):
""" Pretty print a n×n matrix with a result vector n×1. """
n = len(A)
for i in range(0, n):
line = ""
for j in range(0, n+1):
line += str(A[i][j]) + "\t"
if j == n-1:
line += "| "
print(line)
print("")
def pprintPolynomial(A):
""" Pretty print a polynomial. """
line = ""
for i in range(len(x)-1, -1, -1):
if x[i] != 0:
if i == 0:
line += "+" + str(x[i])
else:
if x[i] == 1:
line += "+ x^" + str(i) + "\t"
elif x[i] == -1:
line += "- x^" + str(i) + "\t"
else:
line += "+" + str(x[i]) + "·x^" + str(i) + "\t"
print(line)
def gauss(A):
""" Solve a linear sysem of equations given by a n×n matrix
with a result vector n×1. """
n = len(A)
for i in range(0,n):
# Search for maximum in this column
maxEl = abs(A[i][i])
maxRow = i
for k in range(i+1,n):
if abs(A[k][i]) > maxEl:
maxEl = A[k][i]
maxRow = k
# Swap maximum row with current row (column by column)
for k in range(i,n+1):
tmp = A[maxRow][k]
A[maxRow][k] = A[i][k]
A[i][k] = tmp
# Make all rows below this one 0 in current column
for k in range(i+1,n):
c = -A[k][i]/A[i][i]
for j in range(i,n+1):
if i==j:
A[k][j] = 0
else:
A[k][j] += c * A[i][j]
# Solve equation Ax=b for an upper triangular matrix A
x=[0 for i in range(n)]
for i in range(n-1,-1,-1):
x[i] = A[i][n]/A[i][i]
for k in range(i-1,-1,-1):
A[k][n] -= A[k][i] * x[i]
return x
def setGauss(points):
""" Create a system of equations for gaussian elimination from
a set of points. """
n = len(points) - 1
A = [[0 for i in range(n+2)] for j in range(n+1)]
for i in range(n+1):
x = points[i]["x"]
for j in range(n+1):
A[i][j] = x**j
A[i][n+1] = points[i]["y"]
return A
if __name__ == "__main__":
from fractions import Fraction
# Read input data
points = []
points.append({"x": Fraction(-1), "y": Fraction(1)})
points.append({"x": Fraction(1), "y": Fraction(1)})
points.append({"x": Fraction(2), "y": Fraction(2)})
A = setGauss(points)
# Print input
pprintGaus(A)
# Calculate solution
x = gauss(A)
# Print result
pprintPolynomial(x)
It is also interesting to get the value of \(p(x)\) at any given point \(x \in \mathbb{R}\):
def evaluatePolynomial(p, x):
y = 0
xi = 1
for i, a in enumerate(p):
y += a * xi
xi *= x
return y
Time complexity to get the polynomial: \(\frac{1}{3} n^3 + \mathcal{O}(n^2)\) (where \(n\) is the number of points) Space complexity for the polynomial: \(\mathcal{O}(n^2)\) (where \(n\) is the number of points)
Time complexity to evaluate the value of \(p(x)\) for any \(x \in \mathbb{R}\): \(\mathcal{O}(n)\).
Polynomials
\(\displaystyle \mathbb{R}_n[X] := \left \{p:\mathbb{R} \rightarrow \mathbb{R} | p(x) = \sum_{i=0}^n a_i \cdot x^i \text{ with } a_0, \dots, a_n \in \mathbb{R} \right \}\)
So \(\mathbb{R}_n[X]\) are all polynomials with real coefficients and degree not higher than latex]n\(. $\mathbb{R}_n[X]\) forms a vector space for \(n \in \mathbb{N}_0\). The degree of that vector space is \(n+1\).
What does that mean?
This means you can use different bases for polynomials. The base we usually use for \(\mathbb{R}_n[X]\) is \(\{x^0, x^1, x^2, x^3, \dots\}\), but you can switch the base.
Lagrange interpolation
Define \(n+1\) polynomials like this:
\(\displaystyle L_{i}(x) := \prod_{j=0 \atop j \neq i}^n \frac{x-x_j}{x_i - x_j}\)
This means:
\(\displaystyle L_{i}(x_i) := \prod_{j=0 \atop j \neq i}^n \frac{x_i-x_j}{x_i - x_j} = 1\)
and
\(\displaystyle L_{i}(x_p) := \prod_{j=0 \atop j \neq i}^n \frac{x_p-x_j}{x_i - x_j} = 0 \;\; p \in \{0, \dots, n\} \setminus \{i\}\)
So the polynomial \(y_i \cdot L_{i}(x)\) fits the point \((x_i, y_i)\) and is zero for all other points. This implies that \(\displaystyle p(x) = \sum_{i=0}^n y_i \cdot L_i(x)\) is an interpolation of our data points. The degree of \(p(x)\) is not higher than \(n\). You can see that easily when you take a look at the definition of \(p(x)\) and \(L_i(x)\).
The polynomials \(L_i(x)\) form another base for \(\mathbb{R}_n[X]\).
Lagranges way to interpolate polynomials can be implemented like this:
def lagrangeInterpolation(points):
p = []
for i in range(len(points)):
Li = {"y": points[i]["y"], "polynomial":[]}
for j in range(len(points)):
if j == i:
continue
Li["polynomial"].append({
"sub": points[j]["x"],
"divisor": points[i]["x"] - points[j]["x"]
})
p.append(Li)
return p
def evaluateLagrangePolynomial(p, x):
y = 0
for Li in p:
prod = 1
for term in Li["polynomial"]:
prod *= (x - term["sub"])/term["divisor"]
y += Li["y"]*prod
return y
Time complexity to get the polynomial: \(n^2 + \mathcal{O}(n)\) (where \(n\) is the number of points) Space complexity for the polynomial: \(\mathcal{O}(n^2)\) (where \(n\) is the number of points)
Time complexity to evaluate the value of \(p(x)\) for any \(x \in \mathbb{R}\): \(\mathcal{O}(n^2)\).
Newton interpolation
Define
\(\displaystyle N_i(x) := \prod_{j=0}^{i-1} (x-x_j)\)
So you know that
\(N_i(x) = 0 \Leftrightarrow \exists p \in \{1, \dots, i\}: x = x_{i-p}\)
So your polynomial is
\(\displaystyle p(x) = \sum_{i=0}^n c_i \cdot N_i(x)\)
for the correct \(c_i\). You can do this by solving the following system of equations. Please note that you don't have to store that matrix to get those \(c_i\), divided differences are better.
You can get this lower triangular matrix like this:
def getGaussSystemForNewton(points):
n = len(points) - 1
A = [[0 for i in range(n+2)] for j in range(n+1)]
for j in range(0,n+2):
for i in range(j,n+1):
if j == 0:
A[i][j] = 1
else:
A[i][j] = A[i][j-1]*(points[i]["x"]-points[j-1]["x"])
if j == n+1:
for i in range(0, n):
A[i][j] = points[i]["y"]
return A
From my previous posts about solving equations of upper triangular matrices and lower unitriangular matrices you know that the space complexity of this is in \(\Theta(n^2)\).
According to Wikipedia, you can use Horner's method to evaluate this Polynom in \(\mathcal{O}(n)\). But I really don't want to implement this today.
Interactive example
- Click to add points.
- Ctrl+Click to remove point.
- You can enter a function that is valid JavaScript in the text field.