5.1. Polynomial Collocation (Interpolation/Extrapolation)#

References:

Introduction#

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. (Other important choices are rational functions and “trigonometric polynomials”: sums of multiples of sines and cosines.) Such polynomials can then be used to aproximate derivatives and integrals.

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\). 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)\).

In fact there are infinitely many such polynomials: given one, add to it any polynomial with zeros at all of the \(n+1\) notes. So to make the problem well-posed, we seek the collocating polynomial of lowest degree.

Theorem 5.1 (Existence and uniqueness of a collocating polynomial)

Given \(n+1\) distinct values \(x_i\), \(0 \leq i \leq n\), and corresponding \(y\)-values \(y_i\), there is a unique polynomial \(P\) of degree at most \(n\) with \(P(x_i) = y_i\) for \(0 \leq i \leq 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\).)

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. 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.

In general, this method starts by assuming that the function wanted is a sum of unknown multiples of a collection of known functions. 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\).
(Note: any of the \(c_i\) could be zero, including \(c_n\), in which case the degree is less than \(n\).)
The unknown factors (\(c_0 \cdots c_n\)) are the undetermined coefficients.

Next one states the problem as a system of equations for these undetermined coefficients, and solves them.
Here, we have \(n+1\) conditions to be met:

\[P(x_i) = \sum_{j=0}^n c_j x_i^j = y_i, \quad 0 \leq i \leq 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, and so is equivalent to the case of all \(y_i = 0\) having only the solution with all \(c_i = 0\).

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\). 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\). The theorem is proven.

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. (Aside: here as in many places, Numpy mimics the names and functionality of corresponding Matlab tools.)

Briefly, the algorithm is this (indexing from 0 !)

  • Create the \(n+1 \times n+1\) matrix \(V\) with elements

\[v_{i,j} = x_i^j,\; 0 \leq i \leq n, \, 0 \leq j \leq n\]

and the \(n+1\)-element column vector \(y\) with elements \(y_i\) as above.

  • Solve \(V c = y\) for the vector of coefficients \(c_j\) as above.

I use the name \(V\) because this is called the Vandermonde Matrix.

from numpy import array, linspace, zeros, zeros_like, exp
from  matplotlib.pyplot import figure, plot, title, grid, legend

Example 5.1

As usual, I concoct a first example with known correct answer, by using a polynomial as \(f\):

\[ f(x) = 4 + 7x -2x^2 - 5x^3 + 2x^4 \]

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.)

def f(x):
    return 4 + 7*x - 2*x**2 -5*x**3 + 3*x**4

xnodes = array([1., 2., 7., 5., 4.])  # They do not need to be in order
nnodes = len(xnodes)
n = nnodes-1
print(f"The x nodes 'x_i' are {xnodes}")
ynodes = zeros_like(xnodes)
for i in range(nnodes):
    ynodes[i] = f(xnodes[i])
print(f"The y values at the nodes are {ynodes}")

The x nodes 'x_i' are [1. 2. 7. 5. 4.]
The y values at the nodes are [   7.   18. 5443. 1239.  448.]

The Vandermonde matrix:

V = zeros([nnodes, nnodes])
for i in range(nnodes):
    for j in range(nnodes):
        V[i,j] = xnodes[i]**j

Solve, using our functions seen in earlier sections and gathered in Notebook for generating the module numericalMethods

from numericalMethods import rowReduce, backwardSubstitution

---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
Cell In[7], line 1
----> 1 from numerical_methods import rowReduce, backwardSubstitution

ImportError: cannot import name 'rowReduce' from 'numerical_methods' (/Users/brenton/Library/CloudStorage/OneDrive-CollegeofCharleston/numerical-methods-and-analysis/python-edition/main/numerical_methods.py)

(U, z) = rowReduce(V, ynodes)
c = backwardSubstitution(U, z)
print(f"The coefficients of P are {c}")

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[8], line 1
----> 1 (U, z) = rowReduce(V, ynodes)
      2 c = backwardSubstitution(U, z)
      3 print(f"The coefficients of P are {c}")

NameError: name 'rowReduce' is not defined

We can also check the resulting values of the polynomial:

P = c[0] + c[1]*xnodes + c[2]*xnodes**2 + c[3]*xnodes**3 + c[4]*xnodes**4

for (x, y, P_i) in zip(xnodes, ynodes, P):
    print(f"P({x}) should be {y}; it is {P_i}")

Functions for computing the coefficients and evaluating the polynomials#

We will use this procedure several times, so it time to put it into a functions — and add a pretty printer for polynomials.

def fitpolynomial(x, y):
    """Compute the coefficients c_i of the polynomial of lowest degree that collocates the points (x[i], y[i]).
    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,
    in which case the array has some trailing zeroes.
    The polynomial is thus p(x) = c[0] + c[1]x + ... c[d] x^d where n =length(x)-1, the nominal degree.
    """
    nnodes = len(x)
    n = nnodes - 1
    V = zeros([nnodes, nnodes])
    for i in range(nnodes):
        for j in range(nnodes):
             V[i,j] = x[i]**j
    (U, z) = rowReduce(V, y)
    c = backwardSubstitution(U, z)
    return c

def evaluatepolynomial(x, coeffs):
    # Evaluate the polynomial with coefficients in coeffs  at the points in x.
    npoints = len(x)
    ncoeffs = len(coeffs)
    n = ncoeffs - 1
    powers = linspace(0, n, n+1)
    y = zeros_like(x)
    for i in range(npoints):
        y[i] = sum(coeffs * x[i]**powers)
    return y

def showpolynomial(c):
    print("P(x) = ", end="")
    n = len(c)-1
    print(f"{c[0]:.4}", end="")
    if n > 0:
        coeff = c[1]
        if coeff > 0:
            print(f" + {coeff:.4}x", end="")
        elif coeff < 0:
            print(f" - {-coeff:.4}x", end="")
    if n > 1:
        for j in range(2, len(c)):
            coeff = c[j]
            if coeff > 0:
                print(f" + {coeff:.4}x^{j}", end="")
            elif coeff < 0:
                print(f" - {-coeff:.4}x^{j}", end="")
    print()

While debugging, redo the first example:

c_new = fitpolynomial(xnodes, ynodes)
print(c_new)

P_i_new = evaluatepolynomial(xnodes, c_new)

print(P_i_new)

print(f"The values of P(x_i) are {P_i_new}")

showpolynomial(c_new)

Example 5.2 (\(f(x)\) not a polynomial of degree \(\leq 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: \(x_0 = 1\), \(x_1 = 2\), \(x_2 = 3\) and \(x_3 = 4\)

# Reduce the degree of $P$ to at most 3:
n = 3
xnodes = array([-2.0, 0., 1.0, 2.])
ynodes = zeros_like(xnodes)
for i in range(len(xnodes)):
    ynodes[i] = f(xnodes[i])
print(f"n is now {n}, the nodes are now {xnodes}, with f(x_i) values {ynodes}")
c = fitpolynomial(xnodes, ynodes)
print(f"The coefficients of P are now {c}")
showpolynomial(c)
print(f"The values of P at the nodes are now {evaluatepolynomial(xnodes, c)}")

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.

x = linspace(min(xnodes), max(xnodes))  # defaulting to 50 points, for graphing

figure(figsize=[12,6])
plot(x, f(x), label="y=f(x)")
plot(xnodes, ynodes, "*", label="nodes")
P_n_x = evaluatepolynomial(x, c)
plot(x, P_n_x, label="y = P_n(x)")
legend()
grid(True);

P_error = f(x) - P_n_x
figure(figsize=[12,6])
title("Error in P_n(x)")
plot(x, P_error, label="y=f(x)")
plot(xnodes, zeros_like(xnodes), "*")
grid(True);

Example 5.3 (\(f(x)\) not a polynomial at all)

\(f(x) = e^x\) with five nodes, equally spaced from \(-1\) to \(1\)

def g(x): return exp(x)
a_g = -1.0
b_g = 1.0

n = 3
xnodes_g = linspace(a_g, b_g, n+1)
ynodes_g = zeros_like(xnodes_g)
for i in range(len(xnodes_g)):
    ynodes_g[i] = g(xnodes_g[i])
print(f"{n=}")
print(f"node x values {xnodes_g}")
print(f"node y values {ynodes_g}")
c_g = fitpolynomial(xnodes_g, ynodes_g)
print(f"The coefficients of P are {c_g}")
showpolynomial(c_g)
P_values = evaluatepolynomial(c_g, xnodes_g)
print(f"The values of P(x_i) are {P_values}")

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.

x_g = linspace(a_g - 0.25, b_g + 0.25)  # Go a bit beyond the nodes in each direction

figure(figsize=[14,10])
title("With $g(x) = e^x$")
plot(x_g, g(x_g), label="y = $g(x)$")
plot(xnodes_g, ynodes_g, "*", label="nodes")
P_g = evaluatepolynomial(x_g, c_g)
plot(x_g, P_g, label=f"y = $P_{n}(x)$")
legend()
grid(True);

P_error = g(x_g) - P_g
figure(figsize=[14,10])
title(f"Error in $P_{n}(x)$ for $g(x) = e^x$")
plot(x_g, P_error)
plot(xnodes_g, zeros_like(xnodes_g), "*")
grid(True);

P_error = f(x) - P_m_x
figure(figsize=[14,10])
title("Error in $P_1(x)$ for $f(x)$")
plot(x, P_error, label="y=f(x)")
plot(xnodes, zeros_like(xnodes), "*")
grid(True);

This work is licensed under Creative Commons Attribution-ShareAlike 4.0 International