Suppose you have a list of n+1n+1 given point (xi,yi)(xi,yi) with i∈0,…,ni∈0,…,n and ∀i,j∈0,…,n:i≠j⇒xi≠xj∀i,j∈0,…,n:i≠j⇒xi≠xj.
Now you want to find a polynomial p(x)=n∑i=0ai⋅xip(x)=n∑i=0ai⋅xi that goes through all of the given points. This means:
∀i∈0,…,n:p(xi)=yi∀i∈0,…,n:p(xi)=yi
Existence and uniqueness ¶
Theorem: When you have a list of n+1n+1 points with mutually different xixi there is exactly one polynom of degree ≤n≤n that fits those points.
Proof:
You can formulate a linear system of equations:
(x00⋯xn0 ⋮⋱⋮ x0n⋯xnn)⏟:=A∈R(n+1)×(n+1)(a0 ⋮ an)=(y0 ⋮ yn)(x00⋯xn0 ⋮⋱⋮ x0n⋯xnn):=A∈R(n+1)×(n+1)(a0 ⋮ an)=(y0 ⋮ yn)
A matrix of the form of AA is called Vandermonde matrix. When the data points xixi are mutually different, it is known that the Vandermonde matrix is invertible (source).
This means, the solution (a0…an)T(a0…an)T of this linear equation gives the polynom p(x)=∑ni=0aixip(x)=∑ni=0aixi
So the solution exists and is unique ◼■
Straight forward interpolating polynomials ¶
For this algorithm, I'll find the polynomial in its monomial from p(x)=∑ni=0aixip(x)=∑ni=0aixi. I'll use the matrix AA 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)p(x) at any given point x∈Rx∈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: 13n3+O(n2)13n3+O(n2) (where nn is the number of points) Space complexity for the polynomial: O(n2)O(n2) (where nn is the number of points)
Time complexity to evaluate the value of p(x)p(x) for any x∈Rx∈R: O(n)O(n).
Polynomials ¶
$\displaystyle \mathbb{R}n[X] := \left {p:\mathbb{R} \rightarrow \mathbb{R} | p(x) = \sum \right }$}^n a_i \cdot x^i \text{ with } a_0, \dots, a_n \in \mathbb{R
So Rn[X]Rn[X] are all polynomials with real coefficients and degree not higher than latex]n..\mathbb{R}_n[X]formsavectorspaceforformsavectorspaceforn \in \mathbb{N}_0.Thedegreeofthatvectorspaceis.Thedegreeofthatvectorspaceisn+1$.
What does that mean?
This means you can use different bases for polynomials. The base we usually use for Rn[X]Rn[X] is x0,x1,x2,x3,…x0,x1,x2,x3,…, but you can switch the base.
Lagrange interpolation ¶
Define n+1n+1 polynomials like this:
Li(x):=n∏j=0j≠ix−xjxi−xjLi(x):=n∏j=0j≠ix−xjxi−xj
This means:
Li(xi):=n∏j=0j≠ixi−xjxi−xj=1Li(xi):=n∏j=0j≠ixi−xjxi−xj=1
and
Li(xp):=n∏j=0j≠ixp−xjxi−xj=0p∈0,…,n∖iLi(xp):=n∏j=0j≠ixp−xjxi−xj=0p∈0,…,n∖i
So the polynomial yi⋅Li(x)yi⋅Li(x) fits the point (xi,yi)(xi,yi) and is zero for all other points. This implies that p(x)=n∑i=0yi⋅Li(x)p(x)=n∑i=0yi⋅Li(x) is an interpolation of our data points. The degree of p(x)p(x) is not higher than nn. You can see that easily when you take a look at the definition of p(x)p(x) and Li(x)Li(x).
The polynomials Li(x)Li(x) form another base for Rn[X]Rn[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: n2+O(n)n2+O(n) (where nn is the number of points) Space complexity for the polynomial: O(n2)O(n2) (where nn is the number of points)
Time complexity to evaluate the value of p(x)p(x) for any x∈Rx∈R: O(n2)O(n2).
Newton interpolation ¶
Define
Ni(x):=i−1∏j=0(x−xj)Ni(x):=i−1∏j=0(x−xj)
So you know that
Ni(x)=0⇔∃p∈1,…,i:x=xi−pNi(x)=0⇔∃p∈1,…,i:x=xi−p
So your polynomial is
p(x)=n∑i=0ci⋅Ni(x)p(x)=n∑i=0ci⋅Ni(x)
for the correct cici. 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 cici, divided differences are better.
(10 1(x1−x0) 1(x2−x0)(x2−x0)(x2−x1) ⋮⋮⋱ 1(xn−x0)⋯∏n−1i=0(xn−xi) )⋅(c0 ⋮ cn)=(f0 ⋮ fn)(10 1(x1−x0) 1(x2−x0)(x2−x0)(x2−x1) ⋮⋮⋱ 1(xn−x0)⋯∏n−1i=0(xn−xi) )⋅(c0 ⋮ cn)=(f0 ⋮ fn)
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 Θ(n2).
According to Wikipedia, you can use Horner's method to evaluate this Polynom in 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.