微信客服：xiaoxionga100
微信客服：ITCS521
MATH256 2022-math256代写
时间：2023-04-13
Numerical Methods (5 Parts)MATH256 2022-23
Contents
4 Polynomial Interpolation 4-1
4.1 Theorem of Uniqueness . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4-2
4.2 The Lagrange interpolation formula . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4-2
4.2.1 Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4-3
4.3 Newton’s polynomial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4-4
4.4 Theorem on Divided Difference . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4-6
4.5 Divided difference tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4-8
4.5.1 Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4-8
4.5.2 Further examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4-9
4.6 Evenly spaced nodes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4-9
4.7 Bounding of Interpolation Error . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4-9
4.8 Chebyshev polynomials . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4-11
4.8.1 Examples of first few Chebyshev polynomials . . . . . . . . . . . . . . . . . . . . 4-11
4.8.2 General expression for Tn(t) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4-12
4.8.3 Maxima and minima of Chebyshev polynomials . . . . . . . . . . . . . . . . . . . 4-13
4.8.4 Roots of Chebyshev polynomials . . . . . . . . . . . . . . . . . . . . . . . . . . . 4-14
4.8.5 Example of the normalised Chebyshev polynomial T̂n(t) in Maple . . . . . . . . . . 4-15
4.8.6 Theory of Minimax Property of T̂n(t) . . . . . . . . . . . . . . . . . . . . . . . . 4-15
4.8.7 Interpolation using Chebyshev nodes . . . . . . . . . . . . . . . . . . . . . . . . . 4-16
4.8.8 Maple Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4-17
4.9 Splines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4-17
* Download and save these files before reading on [Optional V1 files contain more graphics]
§4.2 lagrange.mw §4.2 V1_P2_uniquePractice.mw
§4.3 newton_ip.mw §4.3 V1_Lagrange_Newton_Practice.mw
§4.5.2 d_tab.mw §4.7 V1_Newton_NonsmoothF_Error.mw
§4.8.1 chebyshev_demo.mw §4.8.7 V1_Newton_NonsmoothF_Cheby.mw
§4.9 splines.mw §4.9 V1_splines_xy.mw
4 Polynomial Interpolation
It is often useful to fit a polynomial through a set of data points. Given data of (n+ 1) points
{(x0, y0), (x1, y1), . . . , (xn, yn)}
(with xj ≠ xq for j ̸= q) we aim to find a polynomial P of degree n, or Pn(x), such that
P (xj) = yj , j = 0, 1, . . . , n.
The data may come from an experiment, or yj may represent values of some given function y that we
want to approximate using a simpler function, in which case yj = y(xj). The process of finding P is called
solving the interpolation problem, and the points xj are called nodes. We start with a very simple (but
crucial) theorem.
4-1
Numerical Methods (5 Parts)MATH256 2022
4.1 Theorem of Uniqueness
The solution to the interpolation problem is unique. [Note: existence is a separate matter]
Proof
Suppose that P (x) and Q(x) are two solutions to the interpolation problem. Then the difference
R(x) = P (x)−Q(x)
is a polynomial of degree n or smaller, such that R(xj) = 0 for j = 0, 1, . . . , n. Since there are n+ 1 such
values, this is only possible if R(x) = 0 for all x, meaning P (x) and Q(x) are identical; there is at most
only one solution.
Remarks
• The above is for uniqueness and doesn’t prove that a solution exists. Existence is shown below by
construction (the best way of a proof).
• The solution clearly exists if n = 1, because then P (x) is just the straight line through (x0, y0) and
(x1, y1) — a situation we met previously in False-position and Secant methods.
4.2 The Lagrange interpolation formula
Consider the function
Lj(x) =
n∏
q=0
q ̸=j
x− xq
xj − xq =
(x− x0) · · · (x− xj−1)(x− xj+1) · · · (x− xn)
(xj − x0) · · · (xj − xj−1)(xj − xj+1) · · · (xj − xn) .
· · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · ·
Clearly, Lj(xj) = 1, and if p ̸= j there is a factor of (x− xp). Therefore,
Lj(xp) =
{
1 if p = j,
0 if p ̸= j.
The Lagrange interpolation formula is
P (x) =
n∑
j=0
yjLj(x).
This solves the interpolation problem because
P (xp) =
n∑
j=0
yjLj(xp) = yp
for p = 0, 1, . . . , n.
· · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · ·
4-2
Numerical Methods (5 Parts)MATH256 2022-23
Remarks:
• Lj(x) is called a Lagrange polynomial, or a Lagrange basis polynomial.
• Since Lj(x) is a polynomial of degree n, it follows that P (x) is a polynomial whose degree is at most
n (some terms may cancel).
• The Lagrange interpolation formula shows that a solution to the interpolation problem is certain to
exist; we can always fit a polynomial of degree n through n+ 1 data points.
4.2.1 Example
Use the Lagrange formula to find the quadratic polynomial P2 = P (x) that passes through the points
(0, 5), (1, 1) and (2,−1).
Solution
· · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · ·
We have xj = j, for j = 0, 1 and 2, so
L0(x) =
x− x1
x0 − x1
x− x2
x0 − x2 =
x− 1
−1
x− 2
−2 =
(x− 1)(x− 2)
2
L1(x) =
x− x0
x1 − x0
x− x2
x1 − x2 = −x(x− 2)
and
L2(x) =
x− x0
x2 − x0
x− x1
x2 − x1 =
x
2 (x− 1).
Then
P (x) =
2∑
j=0
yjLj(x) =
5(x− 1)(x− 2)
2 − x(x− 2)−
x(x− 1)
2
= x2 − 5x+ 5.
It is easy to check this is correct:
p(0) = 5,
p(1) = 1,
p(2) = 4− 10 + 5 = −1.
· · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · ·
Of course we can solve small interpolation problems from first principles; e.g. write P (x) = ax2 + bx+ c
and use the method of undetermined coefficients:
P (0) = c = 5, P (1) = a+ b+ 5 = 1 and P (2) = 4a+ 2b+ 5 = −1
and solve for a and b. However, the Lagrange formula gives us a systematic way to get the soluton or to
solve the problem for any number of points; see lagrange.mw.
4-3
MATH256 2021–22 Numerical Methods (5 Parts)
4.3 Newton’s polynomial
A downside of the Lagrange formula is that the entire calculation must be repeated if we need to add an
extra interpolation point. We can overcome this problem using Newton’s polynomial
Pn−1(x) = a0 + a1(x− x0) + a2(x− x0)(x− x1) + · · ·+ an−1
n−2∏
j=0
(x− xj),
and importantly
Pn(x) = Pn−1(x) + an
n−1∏
j=0
(x− xj) i.e.
Pn(x) = a0 + a1(x− x0) + a2(x− x0)(x− x1) + · · ·+ an
n−1∏
j=0
(x− xj).
The idea is that all terms on the right-hand side have a factor (x − x0) except the first, so we can set
x = x0 to determine a0. Then we set x = x1 to find a1, etc. Each coefficient aj is obtained without
changing the earlier coefficients. Before proceeding with this, we make one important observation:
xn appears in only one term (the last), so the leading order coefficient in the Newton polynomial
is an.
Now we determine the coefficients aj . The key is to get a convenient and easily memorable formula.
What we find below is that we can derive recurrent formulas but they are not in the most convenient format
(next subsection will reformulate them to a nice form).
· · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · ·
Since P (x0) = a0, we choose
a0 = y0.
Similarly,
P (x1) = a0 + a1(x1 − x0)
= y0 + a1(x1 − x0),
and, since we want P (x1) = y1, we get an unsurprising formula:
a1 =
y1 − y0
x1 − x0 .
The coefficient a1 is a divided difference: a ratio of a change in y to a change in x. For the next term, we
have
P (x2) = a0 + a1(x2 − x0) + a2(x2 − x0)(x2 − x1)
i.e.
y2 = y0 +
y1 − y0
x1 − x0 (x2 − x0) + a2(x2 − x0)(x2 − x1).
Rearranging, we find that
a2 =
1
x2 − x1
(
y2 − y0
x2 − x0 −
y1 − y0
x1 − x0
)
,
which is a second divided difference (a divided difference of divided differences). There is a pattern to
the rearrangements, which we can exploit. Consider the calculation for a3.
4-4
MATH256 2021–22 Numerical Methods (5 Parts)
• Set x = x3. All terms after a3 disappear from the RHS.
P (x3) = a0 + a1(x3 − x0) + a2(x3 − x0)(x3 − x1)
+ a3(x3 − x0)(x3 − x1)(x3 − x2).
• Use P (x3) = y3 and take a0 = y0 to the LHS:
y3 − y0 = a1(x3 − x0) + a2(x3 − x0)(x3 − x1)
+ a3(x3 − x0)(x3 − x1)(x3 − x2).
• All terms on the RHS now have a factor (x3 − x0); divide through:
y3 − y0
x3 − x0 = a1 + a2(x3 − x1) + a3(x3 − x1)(x3 − x2).
Now we have a divided difference on the LHS.
• Take a1 to the LHS and use a1 = (y1 − y0)/(x1 − x0):
y3 − y0
x3 − x0 −
y1 − y0
x1 − x0 = a2(x3 − x1) + a3(x3 − x1)(x3 − x2).
• All terms on the RHS now have a factor (x3 − x1). Divide through:
1
x3 − x1
(
y3 − y0
x3 − x0 −
y1 − y0
x1 − x0
)
= a2 + a3(x3 − x2).
Now we have a second divided difference on the LHS.
• Take a2 to the LHS:
1
x3 − x1
(
y3 − y0
x3 − x0 −
y1 − y0
x1 − x0
)
− 1
x2 − x1
(
y2 − y0
x2 − x0 −
y1 − y0
x1 − x0
)
= a3(x3 − x2).
Evidently, dividing by (x3 − x2) will produce a third divided difference.
· · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · ·
To generalise this, we introduce some new notation.
The leading coefficient in the Newton polynomial through the points (x0, y0), (x1, y1), . . .,
(xn, yn) (i.e. an) is denoted by y[x0, x1, . . . , xn].
Note that earlier coefficients are defined through interpolation problems with fewer points; for example a1 is
the leading coefficient in the Newton polynomial through (x0, y0) and (x1, y1). In general this is not the
same as the linear coefficient in the Newton polynomial through points up to (xn, yn).
4-5
MATH256 2021–22 Numerical Methods (5 Parts)
We know the first few of these (up to y[x0, x1, x2, x3]) are divided differences:
y[x0] = a0 = y0 (of order 0),
y[x0, x1] = a1 =
y1 − y0
x1 − x0 (of order 1),
y[x0, x1, x2] = a2=
1
x2 − x1
(
y2 − y0
x2 − x0 −
y1 − y0
x1 − x0
)
(of order 2),
... ... ... ...
From the formulae, it is clear that
y[x0, x1] = y[x1, x0] and y[x0, x1, x2] = y[x2, x1, x0].
In fact, the ordering of the points can never matter because the interpolation problem has a unique solution.
If we rearrange the points we must still end up with the same coefficient an (the intermediate steps will
change). Hence,
y[x0, x1, . . . , xn] = y[xn, xn−1, . . . , x0], etc.
Now we can determine a recurrence relation for an = y[x0, . . . , xn]. We expect this to be a divided difference
of (n− 1)th order divided differences. Consider Newton’s formula with x = xn; that is
P (xn) = a0 + (xn − x0)a1 + · · ·+ an−1
n−2∏
j=0
(xn − xj) + an
n−1∏
j=0
(xn − xj).
Assume that
aj = y[x0, x1, . . . , xj ] for j = 1, . . . , n− 1
(we already know this holds for j = 0, 1, 2 and 3). Then the process of rearrangements (move a0 to the
LHS, divide by (xn − x1), etc.) produces a divided difference on the left-hand side; this won’t involve xn−1.
That is,
y[x0, x1, . . . , xn−2, xn] = an−1 + an(xn − xn−1).
Since an−1 = y[x0, x1, . . . , xn−1], it now follows that
an = y[x0, x1, . . . , xn] =
y[x0, x1, . . . , xn−2, xn]− y[x0, x1, . . . , xn−1]
xn − xn−1 .
Thus we have proved inductively that an is an nth divided difference.
However, in the numerator of the above RHS, the second term is fine but the first one is not convenient (as
remarked earlier) – this will be fixed in the next subsection.
4.4 Theorem on Divided Difference
It turns out that it doesn’t matter which points are ‘missing’ from the differences on the right-hand side (as
long as they are not the same). Therefore we can do even better with the following theorem.
The nth divided difference for the points (x0, y0), . . ., (xn, yn) is given by
an = y[x0, x1, . . . , xn] =
y[x1, x2, . . . , xn−1, xn]− y[x0, x1, . . . , xn−1]
xn − x0 .
4-6
MATH256 2021–22 Numerical Methods (5 Parts)
Proof
· · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · ·
Suppose that P and Q are polynomials of degree n− 1, such that
P (xj) = yj for j = 0, 1, . . . , n− 1
and
Q(xj) = yj for j = 1, . . . , n.
From the last section, the leading coefficients in P and Q are y[x0, . . . , xn−1] and y[x1, . . . , xn] respectively.
Consider the function
deg n︷ ︸︸ ︷
R(x) =
deg n−1︷ ︸︸ ︷
P (x) +
deg n−1︷ ︸︸ ︷[
Q(x)− P (x)]h(x). (∗)
Evaluating at x = xj , we obtain
R(xj) =
y0 + [Q(x0)− y0]h(x0) if j = 0,
P (xn) + [yn − P (xn)]h(xn) if j = n,
yj otherwise (P,Q both pass thru j = 1, . . . , n− 1).
To correct the two end points, we need h(x0) = 0 and h(xn) = 1 which we achieve by choosing
h(x) = x− x0
xn − x0 .
Therefore R(x) solves the interpolation problem for (x0, y0), . . ., (xn, yn), so its leading coefficient is
an = y[x0, . . . , xn].
Since the leading coefficients in P and Q are
y[x0, . . . , xn−1] and y[x1, . . . , xn],
respectively, matching leading coefficients in (∗), we find that
y[x0, . . . , xn] =
y[x1, . . . , xn]− y[x0, . . . , xn−1]
xn − x0 ,
as required, because the above RHS is the leading coefficient for xn in [Q− P ](x− x0)/(xn − x0).
· · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · ·
Further by uniqueness of interpolation polynomials, we must have that
an =
y[x0, x1, . . . , xn−2, xn]− y[x0, x1, . . . , xn−2, xn−1]
xn − xn−1 =
y[x1, x2, . . . , xn−1, xn]− y[x0, x1, . . . , xn−1]
xn − x0 .
Hence there is nothing special about the indices 0 and n on the above RHS; we could ‘single out’ any
pair of points and end up with a similar formula. We can count which members are missing on the list
x0, x1, . . . , xn−2, xn−1, xn to see a pattern of the fraction.
4-7
MATH256 2021–22 Numerical Methods (5 Parts)
4.5 Divided difference tables
A divided difference always satisfies the recurrence
y[x0, . . . , xn] =
⟨d.d. missing first value⟩ − ⟨d.d. missing second value⟩
⟨second missing x value⟩ − ⟨first missing x value⟩ .
Clearly y[x0, x1, . . . , xn−2, xn] missed the first value xn−1 while y[x0, x1, . . . , xn−2, xn−1] missed the second
value xn. So the denominator is xn − xn−1. For the case of indices 0 and n, y[x1, x2, . . . , xn−1, xn] missed
x0 and y[x0, x1, . . . , xn−1] missed xn so the denominator is xn − x0.
To find the coefficients in the Newton polynomial for a set of data points (x0, y0), (x1, y1), . . . (xn, yn),
we construct a table as follows:
p
0 1 2 3 4
x y
0 x0 y[x0]
1 x1 y[x1] y[x0, x1]
2 x2 y[x2] y[x1, x2] y[x0, x1, x2]
j 3 x3 y[x3] y[x2, x3] y[x1, x2, x3] y[x0, x1, x2, x3]
4 x4 y[x4] y[x3, x4] y[x2, x3, x4] y[x1, x2, x3, x4] y[x0, x1, x2, x3, x4]
... ... ... ... ... ... ... . . .
(Placing the x values to the left of the y values is just for convenience.) Let the entry in row j and column
p be labelled Aj,p = y[xj−p, xj−p+1, . . . , xj ]. Then
Aj,0 = y[xj ] = y(xj) = yj , Aj,1 = y[xj−1, xj ], Aj,2 = y[xj−2, xj−1, xj ].
Elsewhere, Aj,p is constructed using Aj,p−1 and Aj−1,p−1 (i.e. the entry directly to the left and the entry
above this). We just need to be careful to get the correct x values in the denominator. We make two
observations:
• xj does not appear at all before row j. Therefore the ‘missing’ node in Aj−1,p−1 for Aj,p is xj .
• In row j, we start with xj and add additional nodes in decreasing order of index (i.e. xj−1, xj−2, . . .).
Therefore the ‘missing’ node in Aj,p−1 for Aj,p is xj−p.
Hence, for p ̸= 0 e.g. (j,p)=(2,1), (j,p)=(4,2), (j,p)=(4,4) :
Aj,p =
Aj,p−1 −Aj−1,p−1
xj − xj−p .
Clearly the numerator is ‘easy’ (left vs left’s previous row), and the denominator seems less trivial but it is
simply the difference between the tail point of and its head e.g. xj−xj−p for Aj,p = y[xj−p, xj−p+1, . . . , xj ].
4.5.1 Example
Create a divided table for the data
(1.1, 3.7), (0.5, 1.2), and (1.8,−1.4),
and hence generate the Newton polynomial P (x) that passes through these points.
4-8
MATH256 2021–22 Numerical Methods (5 Parts)
Solution
· · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · ·
We create the table as follows:
p
0 1 2
x y
0 1.1 3.7
j 1 0.5 1.2 4.167
2 1.8 −1.4 −2.0 −8.810
The three calculated entries are obtained as follows:
A1,1 =
A1,0 −A0,0
x1 − x0 =
1.2− 3.7
0.5− 1.1 ≈ 4.167
A2,1 =
A2,0 −A1,0
x2 − x1 =
−1.4− 1.2
1.8− 0.5 = −2
A2,2 =
A2,1 −A1,1
x2 − x0 =
−2− 4.167
1.8− 1.1 ≈ −8.810
The Newton polynomial is
P (x) = 3.7 + 4.167(x− 1.1)− 8.810(x− 1.1)(x− 0.5).
· · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · ·
4.5.2 Further examples
See d_tab.mw.
4.6 Evenly spaced nodes
Interpolating polynomials simplify if the nodes are equally spaced, i.e.
x1 − x0 = x2 − x1 = · · · = ∆x.
We shall skip this material as no new maths is involved.
4.7 Bounding of Interpolation Error
Suppose we interpolate the function f(x) for a ≤ x ≤ b using a polynomial P (x) such that P (xj) = f(xj).
We would like to say something about f(x)− P (x) at other points, i.e. how large is the error? We can
derive a useful expression for this provided f(x) is (at least) n+ 1 times differentiable, by considering the
function
F (t) = f(t)− P (t)− [f(x)− P (x)] n∏
j=0
t− xj
x− xj , x ̸= xj .
· · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · ·
4-9
MATH256 2021–22 Numerical Methods (5 Parts)
Step 1: We prove that there exists a point ξ ∈ [a, b] such that F (n+1)(ξ) = 0.
There are at least n+ 2 points in the interval [a, b] at which F (t) = 0: t = xj for j = 0, . . . , n, and t = x.
Between each two roots there must be a turning point, so there are at least n+1 points at which F ′(t) = 0.
There must be at one least turning point of F ′(t) between each of these points, so there are at least n
points at which F ′′(t) = 0. Repeatedly applying this argument shows there exists at least one point ξ, say,
such that F (n+1)(ξ) = 0. Since the location of ξ depends on x, we write ξ = ξ(x).
Step 2: We obtain an expression for F (n+1)(t) and then set t = ξ(x).
We have
F (t) = f(t)− P (t)− f(x)− P (x)n∏
j=0
(x− xj)
[
tn+1 − (x0 + x1 + · · ·+ xn)tn + · · ·
]
so differentiating n+ 1 times yields
F (n+1)(t) = f (n+1)(t)− f(x)− P (x)n∏
j=0
(x− xj)
(n+ 1)!.
Finally, we set t = ξ(x) to eliminate F , and rearrange to obtain
f(x)− P (x) = f
(n+1) (ξ(x))
(n+ 1)!
n∏
j=0
(x− xj).
· · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · ·
Remarks
• Although we derived this under the assumption that x is distinct from xj , both sides vanish if x = xj ,
so it is correct for all x ∈ [a, b].
• We don’t know ξ, but the ratio f (n+1) (ξ(x)) /(n + 1)! usually turns out to be relatively small for
large n (at least if f is a ‘reasonable’ function).
• For equally spaced nodes, the product term
Kn(x) =
n∏
j=0
(x− xj)
can take very large values, particularly near the end-points of the interpolation interval (see figure
4.1). This in turn can cause the approximation to break down if a large number of points are used
(see the demonstration in d_tab.mw).
4-10
MATH256 2021–22 Numerical Methods (5 Parts)
2.
x
2
4
6
8
10
−4 −2 0 2 4
n = 10
lo
g
1
0
|K
n
(x
)|
10
x
2
4
6
8
10
−4 −2 0 2 4
n = 20
lo
g
1
0
|K
n
(x
)|
Figure 4.1: Logarithmic plots of the product |Kn(x)| for x ∈ [−5, 5], using equally spaced nodes.
4.8 Chebyshev polynomials
To improve the accuracy of polynomial interpolation, we can choose the points xj in a more subtle way.
First we will consider interpolation on the interval [−1, 1]; later we can generalise the results to deal with
arbitrary intervals. The aim is then to choose the points xj such the product
Kn(x) =
n∏
j=0
(x− xj)
remains small for all x ∈ [−1, 1]. To achieve this, we consider the function Tn, which has the property that
Tn(cos θ) = cos(nθ), n = 0, 1, . . .
The idea here is that | cos(nθ)| ≤ 1, so we have a bound on the modulus, and cos(nθ) can be expanded to
show that Tn is a polynomial. These are called Chebyshev polynomials of the first kind.∗
4.8.1 Examples of first few Chebyshev polynomials
· · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · ·
The first two Chebyshev polynomials are given by
T0(cos θ) = 1 ⇒ T0(t) = 1
and
T1(cos θ) = cos θ ⇒ T1(t) = t.
At higher orders, we can use addition formulae. Thus
T2(cos θ) = cos(2θ)
= cos2 θ − sin2 θ
= 2 cos2 θ − 1,
∗Many different spellings of the name ‘Chebyshev’ appear in the literature. In general, one should look under ‘T’ (Tchebichef,
etc.) as well as ‘C’, particularly in indexes of older books.
4-11
MATH256 2021–22 Numerical Methods (5 Parts)
so that T2(t) = 2t2 − 1 and
T3(cos θ) = cos(3θ)
= cos(θ + 2θ)
= cos θ cos(2θ)− sin θ sin(2θ)
= cos θ
[
2 cos2 θ − 1]− 2 sin2 θ cos θ
= cos θ
[
2 cos2 θ − 1]− 2(1− cos2 θ) cos θ
= 4 cos3 θ − 3 cos θ,
meaning that T3(t) = 4t3 − 3t.
· · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · ·
Evidently, the leading coefficient can’t be a fraction; in fact it turns out to be 2n−1 for each n (Problem 6.5).
We can make a polynomial with a smaller maximum modulus by dividing through:
T̂n(t) =
Tn(t)
2n−1 .
This is called a normalised Chebyshev polynomial. Experimenting with the plots in chebyshev_demo.mw
suggests there is no monic polynomial† Qn of degree n such that
max
−1≤t≤1
|Qn(t)| < max−1≤t≤1
∣∣T̂n(t)∣∣,
for n = 2 or n = 3. This is called the minimax property (minimising the maximum modulus). We now
need to find a general expression for Tn(t) and show that the minimax property holds for all n.
4.8.2 General expression for Tn(t)
Below we consider two ways of deriving a formula for computing Tn(t).
(1) By recursive formula — simple derivation for a recursive formula:
Since T0(cos θ) = cos0 θ, T0(t) = 1 and T1(cos θ) = cos1 θ, T1(t) = t,
we have T2(cos 2θ) = 2 cos2 θ − 1, T2(t) = 2t2 − 1.
For Tn+1(cos θ) = cos((n+ 1)θ), by keeping cos terms after expansion, we get
cos((n+ 1)θ) = cos((n− 1)θ + 2θ)
= cos(n− 1)θ cos(2θ)− sin(n− 1)θ sin(2θ)
= Tn−1T2 − 2 cos θ sin θ sin(n− 1)θ
= Tn−1T2 − 2T1[cos θ cos(n− 1)θ − cosnθ]
Tn+1(t) = Tn−1T2 − 2T1[T1Tn−1 − Tn]
= 2T1Tn + (T2 − 2T 21 )Tn−1
= 2tTn + (2t2 − 1− t2)Tn−1
= 2tTn − Tn−1.
(2) By binomial coefficients — a direct and non-recursive formula
· · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · ·
†Monic polynomials have leading coefficient 1: antn + an−1tn−1 + · · ·+ a0 is a monic polynomial of degree n if and only if
an = 1.
4-12
MATH256 2021–22 Numerical Methods (5 Parts)
Let θ ∈ R, and consider De Moivre’s formula
einθ = (cos θ + i sin θ)n
=
n∑
j=0
(
n
j
)
(i sin θ)j cosn−j θ,
where the binomial coefficient is given by (
n
j
)
= n!
j!(n− j)! .
Taking the real part, we obtain
cos(nθ) =
n∑
j=0
j even
(
n
j
)
(i sin θ)j cosn−j θ.
Next, we write j = 2p. The upper limit is then ⌊n/2⌋, where ⌊·⌋ means ‘round down’, so
cos(nθ) =
⌊n/2⌋∑
p=0
(
n
2p
)
(i sin θ)2p cosn−2p θ
=
⌊n/2⌋∑
p=0
(
n
2p
)
(cos2 θ − 1)p cosn−2p θ.
Since Tn(cos θ) = cos(nθ), we now have
Tn(t) =
⌊n/2⌋∑
p=0
(
n
2p
)(
t2 − 1
)p
tn−2p.
We can check this against our earlier results:
T0(t) =
0∑
p=0
(
0
2p
)(
t2 − 1
)p
t−2p = 1,
T1(t) =
0∑
p=0
(
1
2p
)(
t2 − 1
)p
t1−2p = t,
T2(t) =
1∑
p=0
(
2
2p
)(
t2 − 1
)p
t2−2p = 2t2 − 1,
T3(t) =
1∑
p=0
(
3
2p
)(
t2 − 1
)p
t3−2p = t3 + 3(t2 − 1)t = 4t3 − 3t.
· · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · ·
4.8.3 Maxima and minima of Chebyshev polynomials
· · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · ·
4-13
MATH256 2021–22 Numerical Methods (5 Parts)
For t ∈ [−1, 1] (the range for cos θ), we know that |Tn(t)| ≤ 1, because
|Tn(cos θ)| = | cos(nθ)| ≤ 1.
The extreme values ±1 are achieved by setting
θ = jπ
n
, j ∈ Z,
which gives t = χj , where
χj = cos
(
jπ
n
)
.
Evidently, χ0 = χ2n = 1 after which the points χj are repeated. Also,
χ2n−j = cos
((2n− j)π
n
)
= cos(2π) cos
(
jπ
n
)
= χj ,
so there are n+ 1 distinct values for χj , obtained by taking j = 0, . . . , n. Two of these are the end points
χ0 = 1 and χn = −1; there are n− 1 stationary points in between. Since Tn(t) is a polynomial of degree
n, there can be no other stationary points. Note that
Tn(χj) = Tn
(
cos(jπ/n)
)
= cos
(
n× jπ
n
)
= (−1)j ,
so χ1, χ3, . . . are minima, whereas χ0, χ2, . . . are maxima.
· · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · ·
4.8.4 Roots of Chebyshev polynomials
To find the points where Tn(t) = 0, we use the formula Tn(cos θ) = cos(nθ) again, this time choosing θ so
that cos(nθ) = 0.
· · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · ·
Thus,
nθ =
(
j + 12
)
π ⇒ θ =
(
j + 12
) π
n
, j ∈ Z,
which gives us
tj = cos
((
j + 12
) π
n
)
.
Again, these repeat starting at j = 2n, but also
t2n−1−j = cos
((
2n− j − 1 + 12
) π
n
)
= cos
((
−j − 12
) π
n
)
= tj ,
meaning only the points t0, t1, . . . , tn−1 are distinct. Since Tn(t) is a polynomial of degree n, it follows
that these are all of the roots.
· · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · ·
4-14
MATH256 2021–22 Numerical Methods (5 Parts)
4.8.5 Example of the normalised Chebyshev polynomial T̂n(t) in Maple
Set a value for n and then use the command
plot( ChebyshevT( n , t ) / 2^(n-1) , t = -1.1 .. 1.1 )
to make a plot of a normalised Chebyshev polynomial in Maple. Observe the pattern of roots and stationary
points.
4.8.6 Theory of Minimax Property of T̂n(t)
Amongst all the monic polynomials of degree n (with real coefficients), the normalised Chebyshev polynomial
T̂n(t) =
Tn(t)
2n−1
has the smallest maximum modulus for t ∈ [−1, 1].
Proof
· · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · ·
Clearly there is only one monic polynomial of degree zero, so we may assume without loss that n ≥ 1.
Suppose there exists a monic polynomial Q(t) of degree n such that
max
−1≤t≤1
|Q(t)| < max
−1≤t≤1
∣∣T̂n(t)∣∣ (∗)
and consider the difference
R(t) = T̂n(t)−Q(t).
Since T̂n(t) and Q(t) are monic, the terms in tn cancel, meaning R(t) is a polynomial whose degree is at
most n− 1. At the maxima of the Chebyshev polynomial (section 4.8.3), we have
T̂n(χ2j) = max−1≤t≤1 T̂n(t) =
1
2n−1
and at the minima
T̂n(χ2j−1) = min−1≤t≤1 T̂n(t) = −
1
2n−1 .
Hence, R(χj) > 0 when j is even, and R(χj) < 0 when j is odd. Between each pair of these values, there
must be a point at which R(t) = 0. However, there are n + 1 points χj , so there must be n points at
which R(t) = 0. This is impossible unless R(t) = 0 for all t, in which case Q(t) = T̂n(t), but this is also
impossible in view of (∗). Therefore we have a contradiction, so our assumption (the existence of Q(t)
satisfying (∗)) is false.
· · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · ·
4-15
MATH256 2021–22 Numerical Methods (5 Parts)
4.8.7 Interpolation using Chebyshev nodes
Recall that the error in interpolating f(x) using a polynomial P (x) of degree n is given by
f(x)− P (x) = f
(n+1) (ξ(x))
(n+ 1)! Kn(x),
where
Kn(x) =
n∏
j=0
(x− xj)
(see section 4.7). We will minimise the maximum value of |Kn| by creating a Chebyshev polynomial on the
right-hand side. Since there are n+ 1 nodes xj , we will use the roots of Tn+1(t). These are given by
tj = cos
((
j + 12
) π
n+ 1
)
, j = 0, . . . , n.
Since cos θ is a decreasing function for 0 ≤ θ ≤ π, these are arranged in descending order. Therefore we
map the interval [−1, 1] (for t) to [a, b] (for x) using the linear transformation
x = a− b2 t+
a+ b
2 .
In particular, t = −1 and t = 1 correspond to x = b and x = a, respectively. The Chebyshev nodes for
[a, b] are given by
xj =
a− b
2 tj +
a+ b
2 , j = 0, 1, . . . n.
Once we have calculated the nodes, we can construct the Newton polynomial using a table of divided
differences. Let us now consider the error in this approximation.
· · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · ·
Using Chebyshev nodes and writing Kn in terms of t, we obtain
Kn(x) =
n∏
j=0
(
a− b
2 t+
a+ b
2 −
a− b
2 tj −
a+ b
2
)
=
n∏
j=0
(
a− b
2 (t− tj)
)
=
(
a− b
2
)n+1 n∏
j=0
(t− tj)
=
(
a− b
2
)n+1
T̂n+1(t),
because the points tj are the roots of T̂n+1(t). Finally, since |T̂n+1(t)| ≤ 2−n, we arrive at the final result
|f(x)− P (x)| ≤
∣∣f (n+1) (ξ(x))∣∣
(n+ 1)! × 2
(
b− a
4
)n+1
.
· · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · ·
4-16
MATH256 2021–22 Numerical Methods (5 Parts)
Remarks
• By switching from equally spaced nodes to Chebyshev nodes, we have replaced the factor ∏nj=0(x−xj)
(about which we know very little) with 2((a− b)/4)n+1 which can easily be calculated.
• The error formula now shows explicitly that increasing b − a makes accurate interpolation more
difficult, but increasing n yields a better result (provided f (n+1)(ξ(x)) behaves itself).
• Chebyshev nodes will generally produce much better results than equally spaced nodes.
4.8.8 Maple Examples
See d_tab.mw.
4.9 Splines
If we need to fit a curve through experimental data, we may not have the option of choosing the distribution
of the points, so we can no longer use Chebyshev nodes. Instead, we can use a piecewise fit, with a function
called a spline between each pair of points. This avoids the need for polynomials of very high degree.
Denote the spline from xj to xj+1 by Sj(x), so that
Sj(xj) = yj and Sj(xj+1) = yj+1, j = 0, 1, . . . , n− 1;
see figure 4.2. For a smooth curve, we need continuous gradient and curvature‡, so the first and second
derivatives at the end of each spline must match those at the start of the next, i.e.
S′j(xj+1) = S′j+1(xj+1) and S′′j (xj+1) = S′′j+1(xj+1), j = 0, 1, . . . , n− 2.
This gives us four conditions for each spline (two data values and two matching derivatives), except the
last, which has only two. The simplest function that can satisfy four such conditions is a cubic polynomial.
Before trying to determine these, let us consider the numbers of equations and unknowns.
• We have n cubic polynomials so there are 4n coefficients to find.
• Fitting a spline through a data point generates an equation that does not involve coefficients from
any other spline. We can use this to eliminate one unknown. Each spline passes through two data
points, so 2n unknowns are eliminated and 2n remain.
• Matching first and second derivatives at interior points relates each spline to an adjacent spline. Since
there are n− 1 interior points, this generates a system of 2(n− 1) equations.
• If we make two assumptions, we will have 2(n− 1) equations for 2(n− 1) unknowns. It turns out
that we can eliminate half of the coefficients from this system.
Next we find equations for the splines. We start by writing
Sj(x) = aj(x− xj)3 + bj(x− xj)2 + cj(x− xj) + dj
Clearly Sj(xj) = yj implies dj = yj so there are only 3n unknowns left: {aj , bj , cj}.
· · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · ·
‡The curvature of a function s(x) is given by κ(x) = s′′(x)
(
1 + [s′(x)]2
)−3/2.
4-17
MATH256 2021–22 Numerical Methods (5 Parts)
xx0 x1 x2 x3
x
y0
y1
y2
y3
S0
S1
S2
Matching
derivatives here
No matching
derivatives here
Figure 4.2: A cubic spline interpolation through the points (xj , yj), j = 0, 1, 2, 3.
0th Derivatives. Evaluating at x = xj+1 and writing hj = xj+1 − xj , we obtain
Sj(xj+1) = yj+1 = ajh3j + bjh2j + cjhj + yj ,
and using the definition of a divided difference this becomes
ajh
2
j + bjhj + cj = y[xj , xj+1], j = 0, 1, . . . , n− 1. (i)
Note Sj(xj) = yj was satisfied and Sj(xj+1) = yj+1 is also satisfied for j = 0, 1, . . . , n− 1.
1st Derivatives. Now consider the derivatives of the splines:
S′j(x) = 3aj(x− xj)2 + 2bj(x− xj) + cj
and
S′j+1(x) = 3aj+1(x− xj+1)2 + 2bj+1(x− xj+1) + cj+1.
Matching S′j(x) = S′j+1(x) at x = xj+1, we find that
3ajh2j + 2bjhj + cj = cj+1, j = 0, 1, . . . , n− 1 (ii)
where j = 0, 1, . . . , n− 2 is expected but we add an extra match at xn to allow j = n− 1 i.e. to allow Sn
to exist (since Sn(x) sits beyond xn i.e. it does not exist).
2nd Derivatives. For the second derivative,
S′′j (x) = 6aj(x− xj) + 2bj
and
S′′j+1(x) = 6aj+1(x− xj+1) + 2bj+1.
Equating these at x = xj+1 yields
3ajhj + bj = bj+1, j = 0, 1, . . . , n− 1. (iii)
4-18
MATH256 2021–22 Numerical Methods (5 Parts)
Here as before, j = 0, 1, . . . , n− 2 is expected but we add an extra match at xn to allow j = n− 1 (since
Sn(x) does not exist).
To summarise, from (i),(ii),(iii), there are 3n equations and the total number of unknowns is 3n+ 2.
Here two of the unknowns bn, cn are extra, each having its own equation in (ii) and (iii) respectively.
• Theoretically we could fix any two of these 3n + 2 unknowns to have a square system and hence a
hopefully unique solution by solving an 3n × 3n system. Below we first simplify these 3n equations by
as much elimination as possible. Then we decide which 2 unknowns should be fixed to leave us a most
convenient system.
We can use aj from (iii) to eliminate aj from (i)-(ii) and then cj next. We are then left with
(i)⇒ bj+1 + 2bj3 hj + cj = y[xj , xj+1] ⇒ cj = y[xj , xj+1]−
bj+1 + 2bj
3 hj ,
⇒ cj+1 = y[xj+1, xj+2]− bj+2 + 2bj+13 hj ....................................................... (∗)
(ii)⇒ (bj+1 + bj)hj + cj = cj+1........................................................................... (†)
Note that the index j = 0, 1, . . . , n− 2 holds in (*) for cj+1 while it was j = 0, 1, . . . , n− 1 with cj . Hence,
(bj+1 + bj)hj + y[xj , xj+1]− bj+1 + 2bj3 hj︸ ︷︷ ︸
cj
= y[xj+1, xj+2]− bj+2 + 2bj+13 hj+1︸ ︷︷ ︸
cj+1
is only valid for j = 0, 1, . . . , n− 2. Gathering all bj terms, we obtain
(†)⇒ hj+1bj+2 + 2(hj+1 + hj)bj+1 + hjbj = 3
(
y[xj+1, xj+2]− y[xj , xj+1]
)
.
· · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · ·
The range of indices j = 0, 1, . . . , n−2 implies that we have n−1 equations for n+1 unknowns b0, b1, . . . , bn
or we have two extra unknowns. This is caused by the fact that we do not impose conditions on the first
and last splines involving end points. To get around this, we must make a ‘closure assumption’ at end
points, the simplest of which is (mainly to maintain a nicely ‘structured’ matrix)
b0 = bn = 0.
The equations for bj can then be written in matrix-vector form (note: A is symmetric, diagonally dominant)
2(h0 + h1) h1
h1 2(h1 + h2) h2
h2 2(h2 + h3) h3
. . . . . . . . .
hn−3 2(hn−3 + hn−2) hn−2
hn−2 2(hn−2 + hn−1)
b1
b2
b3
...
bn−2
bn−1
= 3
y[x1, x2]− y[x0, x1]
y[x2, x3]− y[x1, x2]
y[x3, x4]− y[x2, x3]
...
y[xn−2, xn−1]− y[xn−3, xn−2]
y[xn−1, xn]− y[xn−2, xn−1]
.
4-19
MATH256 2021–22 Numerical Methods (5 Parts)
If you think about it, setting b1 = b2 = 0 would ‘ruin’ the nice tri-diagonal structure of matrix A.
After solving this for bj ’s, we can construct the splines (for j = 0, 1, . . . , n− 1) using
aj =
bj+1 − bj
3hj
and cj = y[xj , xj+1]− bj+1 + 2bj3 hj .
Remarks
• The coefficient bn, though used in computing an−1, cn−1, doesn’t appear directly in the expressions
for the splines, but it does have a meaning:
S′′n−1(xn) = 6an−1hn−1 + 2bn−1 = 2bn,
so its value is half the second derivative at the right end of the curve.
• Similarly, S′′0 (x0) = 2b0, so the closure assumption b0 = bn = 0 causes the second derivative to vanish
at both ends. Other closure assumptions are possible (see Problem 7.3).
• Now is the time to review how the error analysis from §4.7 (Bounding of Interpolation Error) applies
to Sj(x). For some ξ ∈ [a, b], the upper bound will involve |f (4)(ξ)| which is expected to be MUCH
smaller than |f (n+1)(ξ)| since n can be arbitrarily large.
• Final thought: Here Sj(x) is a degree-3 polynomial, smoothly pieced together by continuity conditions
at interior nodes. Hence if we want a difference degree (e.g. quadratic of degree 2 or quartic of
degree 4, or quintic of degree 5) for Sj(x), a similar procedure can be devised to construct them.
However, linear polynomials for Sj(x) cannot have continuity at nodes – why?
Examples
See splines.mw.
4-20
Contents
4 Polynomial Interpolation 4-1
4.1 Theorem of Uniqueness . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4-2
4.2 The Lagrange interpolation formula . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4-2
4.2.1 Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4-3
4.3 Newton’s polynomial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4-4
4.4 Theorem on Divided Difference . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4-6
4.5 Divided difference tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4-8
4.5.1 Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4-8
4.5.2 Further examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4-9
4.6 Evenly spaced nodes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4-9
4.7 Bounding of Interpolation Error . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4-9
4.8 Chebyshev polynomials . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4-11
4.8.1 Examples of first few Chebyshev polynomials . . . . . . . . . . . . . . . . . . . . 4-11
4.8.2 General expression for Tn(t) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4-12
4.8.3 Maxima and minima of Chebyshev polynomials . . . . . . . . . . . . . . . . . . . 4-13
4.8.4 Roots of Chebyshev polynomials . . . . . . . . . . . . . . . . . . . . . . . . . . . 4-14
4.8.5 Example of the normalised Chebyshev polynomial T̂n(t) in Maple . . . . . . . . . . 4-15
4.8.6 Theory of Minimax Property of T̂n(t) . . . . . . . . . . . . . . . . . . . . . . . . 4-15
4.8.7 Interpolation using Chebyshev nodes . . . . . . . . . . . . . . . . . . . . . . . . . 4-16
4.8.8 Maple Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4-17
4.9 Splines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4-17
* Download and save these files before reading on [Optional V1 files contain more graphics]
§4.2 lagrange.mw §4.2 V1_P2_uniquePractice.mw
§4.3 newton_ip.mw §4.3 V1_Lagrange_Newton_Practice.mw
§4.5.2 d_tab.mw §4.7 V1_Newton_NonsmoothF_Error.mw
§4.8.1 chebyshev_demo.mw §4.8.7 V1_Newton_NonsmoothF_Cheby.mw
§4.9 splines.mw §4.9 V1_splines_xy.mw
4 Polynomial Interpolation
It is often useful to fit a polynomial through a set of data points. Given data of (n+ 1) points
{(x0, y0), (x1, y1), . . . , (xn, yn)}
(with xj ≠ xq for j ̸= q) we aim to find a polynomial P of degree n, or Pn(x), such that
P (xj) = yj , j = 0, 1, . . . , n.
The data may come from an experiment, or yj may represent values of some given function y that we
want to approximate using a simpler function, in which case yj = y(xj). The process of finding P is called
solving the interpolation problem, and the points xj are called nodes. We start with a very simple (but
crucial) theorem.
4-1
Numerical Methods (5 Parts)MATH256 2022
4.1 Theorem of Uniqueness
The solution to the interpolation problem is unique. [Note: existence is a separate matter]
Proof
Suppose that P (x) and Q(x) are two solutions to the interpolation problem. Then the difference
R(x) = P (x)−Q(x)
is a polynomial of degree n or smaller, such that R(xj) = 0 for j = 0, 1, . . . , n. Since there are n+ 1 such
values, this is only possible if R(x) = 0 for all x, meaning P (x) and Q(x) are identical; there is at most
only one solution.
Remarks
• The above is for uniqueness and doesn’t prove that a solution exists. Existence is shown below by
construction (the best way of a proof).
• The solution clearly exists if n = 1, because then P (x) is just the straight line through (x0, y0) and
(x1, y1) — a situation we met previously in False-position and Secant methods.
4.2 The Lagrange interpolation formula
Consider the function
Lj(x) =
n∏
q=0
q ̸=j
x− xq
xj − xq =
(x− x0) · · · (x− xj−1)(x− xj+1) · · · (x− xn)
(xj − x0) · · · (xj − xj−1)(xj − xj+1) · · · (xj − xn) .
· · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · ·
Clearly, Lj(xj) = 1, and if p ̸= j there is a factor of (x− xp). Therefore,
Lj(xp) =
{
1 if p = j,
0 if p ̸= j.
The Lagrange interpolation formula is
P (x) =
n∑
j=0
yjLj(x).
This solves the interpolation problem because
P (xp) =
n∑
j=0
yjLj(xp) = yp
for p = 0, 1, . . . , n.
· · · · · · · · · · · · · · · · · · · · · · · · · · · ·
4-2
Numerical Methods (5 Parts)MATH256 2022-23
Remarks:
• Lj(x) is called a Lagrange polynomial, or a Lagrange basis polynomial.
• Since Lj(x) is a polynomial of degree n, it follows that P (x) is a polynomial whose degree is at most
n (some terms may cancel).
• The Lagrange interpolation formula shows that a solution to the interpolation problem is certain to
exist; we can always fit a polynomial of degree n through n+ 1 data points.
4.2.1 Example
Use the Lagrange formula to find the quadratic polynomial P2 = P (x) that passes through the points
(0, 5), (1, 1) and (2,−1).
Solution
· · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · ·
We have xj = j, for j = 0, 1 and 2, so
L0(x) =
x− x1
x0 − x1
x− x2
x0 − x2 =
x− 1
−1
x− 2
−2 =
(x− 1)(x− 2)
2
L1(x) =
x− x0
x1 − x0
x− x2
x1 − x2 = −x(x− 2)
and
L2(x) =
x− x0
x2 − x0
x− x1
x2 − x1 =
x
2 (x− 1).
Then
P (x) =
2∑
j=0
yjLj(x) =
5(x− 1)(x− 2)
2 − x(x− 2)−
x(x− 1)
2
= x2 − 5x+ 5.
It is easy to check this is correct:
p(0) = 5,
p(1) = 1,
p(2) = 4− 10 + 5 = −1.
· · · · · · · · · · · · · · · · · · · · · · · · · · · ·
Of course we can solve small interpolation problems from first principles; e.g. write P (x) = ax2 + bx+ c
and use the method of undetermined coefficients:
P (0) = c = 5, P (1) = a+ b+ 5 = 1 and P (2) = 4a+ 2b+ 5 = −1
and solve for a and b. However, the Lagrange formula gives us a systematic way to get the soluton or to
solve the problem for any number of points; see lagrange.mw.
4-3
MATH256 2021–22 Numerical Methods (5 Parts)
4.3 Newton’s polynomial
A downside of the Lagrange formula is that the entire calculation must be repeated if we need to add an
extra interpolation point. We can overcome this problem using Newton’s polynomial
Pn−1(x) = a0 + a1(x− x0) + a2(x− x0)(x− x1) + · · ·+ an−1
n−2∏
j=0
(x− xj),
and importantly
Pn(x) = Pn−1(x) + an
n−1∏
j=0
(x− xj) i.e.
Pn(x) = a0 + a1(x− x0) + a2(x− x0)(x− x1) + · · ·+ an
n−1∏
j=0
(x− xj).
The idea is that all terms on the right-hand side have a factor (x − x0) except the first, so we can set
x = x0 to determine a0. Then we set x = x1 to find a1, etc. Each coefficient aj is obtained without
changing the earlier coefficients. Before proceeding with this, we make one important observation:
xn appears in only one term (the last), so the leading order coefficient in the Newton polynomial
is an.
Now we determine the coefficients aj . The key is to get a convenient and easily memorable formula.
What we find below is that we can derive recurrent formulas but they are not in the most convenient format
(next subsection will reformulate them to a nice form).
· · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · ·
Since P (x0) = a0, we choose
a0 = y0.
Similarly,
P (x1) = a0 + a1(x1 − x0)
= y0 + a1(x1 − x0),
and, since we want P (x1) = y1, we get an unsurprising formula:
a1 =
y1 − y0
x1 − x0 .
The coefficient a1 is a divided difference: a ratio of a change in y to a change in x. For the next term, we
have
P (x2) = a0 + a1(x2 − x0) + a2(x2 − x0)(x2 − x1)
i.e.
y2 = y0 +
y1 − y0
x1 − x0 (x2 − x0) + a2(x2 − x0)(x2 − x1).
Rearranging, we find that
a2 =
1
x2 − x1
(
y2 − y0
x2 − x0 −
y1 − y0
x1 − x0
)
,
which is a second divided difference (a divided difference of divided differences). There is a pattern to
the rearrangements, which we can exploit. Consider the calculation for a3.
4-4
MATH256 2021–22 Numerical Methods (5 Parts)
• Set x = x3. All terms after a3 disappear from the RHS.
P (x3) = a0 + a1(x3 − x0) + a2(x3 − x0)(x3 − x1)
+ a3(x3 − x0)(x3 − x1)(x3 − x2).
• Use P (x3) = y3 and take a0 = y0 to the LHS:
y3 − y0 = a1(x3 − x0) + a2(x3 − x0)(x3 − x1)
+ a3(x3 − x0)(x3 − x1)(x3 − x2).
• All terms on the RHS now have a factor (x3 − x0); divide through:
y3 − y0
x3 − x0 = a1 + a2(x3 − x1) + a3(x3 − x1)(x3 − x2).
Now we have a divided difference on the LHS.
• Take a1 to the LHS and use a1 = (y1 − y0)/(x1 − x0):
y3 − y0
x3 − x0 −
y1 − y0
x1 − x0 = a2(x3 − x1) + a3(x3 − x1)(x3 − x2).
• All terms on the RHS now have a factor (x3 − x1). Divide through:
1
x3 − x1
(
y3 − y0
x3 − x0 −
y1 − y0
x1 − x0
)
= a2 + a3(x3 − x2).
Now we have a second divided difference on the LHS.
• Take a2 to the LHS:
1
x3 − x1
(
y3 − y0
x3 − x0 −
y1 − y0
x1 − x0
)
− 1
x2 − x1
(
y2 − y0
x2 − x0 −
y1 − y0
x1 − x0
)
= a3(x3 − x2).
Evidently, dividing by (x3 − x2) will produce a third divided difference.
· · · · · · · · · · · · · · · · · · · · · · · · · · · ·
To generalise this, we introduce some new notation.
The leading coefficient in the Newton polynomial through the points (x0, y0), (x1, y1), . . .,
(xn, yn) (i.e. an) is denoted by y[x0, x1, . . . , xn].
Note that earlier coefficients are defined through interpolation problems with fewer points; for example a1 is
the leading coefficient in the Newton polynomial through (x0, y0) and (x1, y1). In general this is not the
same as the linear coefficient in the Newton polynomial through points up to (xn, yn).
4-5
MATH256 2021–22 Numerical Methods (5 Parts)
We know the first few of these (up to y[x0, x1, x2, x3]) are divided differences:
y[x0] = a0 = y0 (of order 0),
y[x0, x1] = a1 =
y1 − y0
x1 − x0 (of order 1),
y[x0, x1, x2] = a2=
1
x2 − x1
(
y2 − y0
x2 − x0 −
y1 − y0
x1 − x0
)
(of order 2),
... ... ... ...
From the formulae, it is clear that
y[x0, x1] = y[x1, x0] and y[x0, x1, x2] = y[x2, x1, x0].
In fact, the ordering of the points can never matter because the interpolation problem has a unique solution.
If we rearrange the points we must still end up with the same coefficient an (the intermediate steps will
change). Hence,
y[x0, x1, . . . , xn] = y[xn, xn−1, . . . , x0], etc.
Now we can determine a recurrence relation for an = y[x0, . . . , xn]. We expect this to be a divided difference
of (n− 1)th order divided differences. Consider Newton’s formula with x = xn; that is
P (xn) = a0 + (xn − x0)a1 + · · ·+ an−1
n−2∏
j=0
(xn − xj) + an
n−1∏
j=0
(xn − xj).
Assume that
aj = y[x0, x1, . . . , xj ] for j = 1, . . . , n− 1
(we already know this holds for j = 0, 1, 2 and 3). Then the process of rearrangements (move a0 to the
LHS, divide by (xn − x1), etc.) produces a divided difference on the left-hand side; this won’t involve xn−1.
That is,
y[x0, x1, . . . , xn−2, xn] = an−1 + an(xn − xn−1).
Since an−1 = y[x0, x1, . . . , xn−1], it now follows that
an = y[x0, x1, . . . , xn] =
y[x0, x1, . . . , xn−2, xn]− y[x0, x1, . . . , xn−1]
xn − xn−1 .
Thus we have proved inductively that an is an nth divided difference.
However, in the numerator of the above RHS, the second term is fine but the first one is not convenient (as
remarked earlier) – this will be fixed in the next subsection.
4.4 Theorem on Divided Difference
It turns out that it doesn’t matter which points are ‘missing’ from the differences on the right-hand side (as
long as they are not the same). Therefore we can do even better with the following theorem.
The nth divided difference for the points (x0, y0), . . ., (xn, yn) is given by
an = y[x0, x1, . . . , xn] =
y[x1, x2, . . . , xn−1, xn]− y[x0, x1, . . . , xn−1]
xn − x0 .
4-6
MATH256 2021–22 Numerical Methods (5 Parts)
Proof
· · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · ·
Suppose that P and Q are polynomials of degree n− 1, such that
P (xj) = yj for j = 0, 1, . . . , n− 1
and
Q(xj) = yj for j = 1, . . . , n.
From the last section, the leading coefficients in P and Q are y[x0, . . . , xn−1] and y[x1, . . . , xn] respectively.
Consider the function
deg n︷ ︸︸ ︷
R(x) =
deg n−1︷ ︸︸ ︷
P (x) +
deg n−1︷ ︸︸ ︷[
Q(x)− P (x)]h(x). (∗)
Evaluating at x = xj , we obtain
R(xj) =
y0 + [Q(x0)− y0]h(x0) if j = 0,
P (xn) + [yn − P (xn)]h(xn) if j = n,
yj otherwise (P,Q both pass thru j = 1, . . . , n− 1).
To correct the two end points, we need h(x0) = 0 and h(xn) = 1 which we achieve by choosing
h(x) = x− x0
xn − x0 .
Therefore R(x) solves the interpolation problem for (x0, y0), . . ., (xn, yn), so its leading coefficient is
an = y[x0, . . . , xn].
Since the leading coefficients in P and Q are
y[x0, . . . , xn−1] and y[x1, . . . , xn],
respectively, matching leading coefficients in (∗), we find that
y[x0, . . . , xn] =
y[x1, . . . , xn]− y[x0, . . . , xn−1]
xn − x0 ,
as required, because the above RHS is the leading coefficient for xn in [Q− P ](x− x0)/(xn − x0).
· · · · · · · · · · · · · · · · · · · · · · · · · · · ·
Further by uniqueness of interpolation polynomials, we must have that
an =
y[x0, x1, . . . , xn−2, xn]− y[x0, x1, . . . , xn−2, xn−1]
xn − xn−1 =
y[x1, x2, . . . , xn−1, xn]− y[x0, x1, . . . , xn−1]
xn − x0 .
Hence there is nothing special about the indices 0 and n on the above RHS; we could ‘single out’ any
pair of points and end up with a similar formula. We can count which members are missing on the list
x0, x1, . . . , xn−2, xn−1, xn to see a pattern of the fraction.
4-7
MATH256 2021–22 Numerical Methods (5 Parts)
4.5 Divided difference tables
A divided difference always satisfies the recurrence
y[x0, . . . , xn] =
⟨d.d. missing first value⟩ − ⟨d.d. missing second value⟩
⟨second missing x value⟩ − ⟨first missing x value⟩ .
Clearly y[x0, x1, . . . , xn−2, xn] missed the first value xn−1 while y[x0, x1, . . . , xn−2, xn−1] missed the second
value xn. So the denominator is xn − xn−1. For the case of indices 0 and n, y[x1, x2, . . . , xn−1, xn] missed
x0 and y[x0, x1, . . . , xn−1] missed xn so the denominator is xn − x0.
To find the coefficients in the Newton polynomial for a set of data points (x0, y0), (x1, y1), . . . (xn, yn),
we construct a table as follows:
p
0 1 2 3 4
x y
0 x0 y[x0]
1 x1 y[x1] y[x0, x1]
2 x2 y[x2] y[x1, x2] y[x0, x1, x2]
j 3 x3 y[x3] y[x2, x3] y[x1, x2, x3] y[x0, x1, x2, x3]
4 x4 y[x4] y[x3, x4] y[x2, x3, x4] y[x1, x2, x3, x4] y[x0, x1, x2, x3, x4]
... ... ... ... ... ... ... . . .
(Placing the x values to the left of the y values is just for convenience.) Let the entry in row j and column
p be labelled Aj,p = y[xj−p, xj−p+1, . . . , xj ]. Then
Aj,0 = y[xj ] = y(xj) = yj , Aj,1 = y[xj−1, xj ], Aj,2 = y[xj−2, xj−1, xj ].
Elsewhere, Aj,p is constructed using Aj,p−1 and Aj−1,p−1 (i.e. the entry directly to the left and the entry
above this). We just need to be careful to get the correct x values in the denominator. We make two
observations:
• xj does not appear at all before row j. Therefore the ‘missing’ node in Aj−1,p−1 for Aj,p is xj .
• In row j, we start with xj and add additional nodes in decreasing order of index (i.e. xj−1, xj−2, . . .).
Therefore the ‘missing’ node in Aj,p−1 for Aj,p is xj−p.
Hence, for p ̸= 0 e.g. (j,p)=(2,1), (j,p)=(4,2), (j,p)=(4,4) :
Aj,p =
Aj,p−1 −Aj−1,p−1
xj − xj−p .
Clearly the numerator is ‘easy’ (left vs left’s previous row), and the denominator seems less trivial but it is
simply the difference between the tail point of and its head e.g. xj−xj−p for Aj,p = y[xj−p, xj−p+1, . . . , xj ].
4.5.1 Example
Create a divided table for the data
(1.1, 3.7), (0.5, 1.2), and (1.8,−1.4),
and hence generate the Newton polynomial P (x) that passes through these points.
4-8
MATH256 2021–22 Numerical Methods (5 Parts)
Solution
· · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · ·
We create the table as follows:
p
0 1 2
x y
0 1.1 3.7
j 1 0.5 1.2 4.167
2 1.8 −1.4 −2.0 −8.810
The three calculated entries are obtained as follows:
A1,1 =
A1,0 −A0,0
x1 − x0 =
1.2− 3.7
0.5− 1.1 ≈ 4.167
A2,1 =
A2,0 −A1,0
x2 − x1 =
−1.4− 1.2
1.8− 0.5 = −2
A2,2 =
A2,1 −A1,1
x2 − x0 =
−2− 4.167
1.8− 1.1 ≈ −8.810
The Newton polynomial is
P (x) = 3.7 + 4.167(x− 1.1)− 8.810(x− 1.1)(x− 0.5).
· · · · · · · · · · · · · · · · · · · · · · · · · · · ·
4.5.2 Further examples
See d_tab.mw.
4.6 Evenly spaced nodes
Interpolating polynomials simplify if the nodes are equally spaced, i.e.
x1 − x0 = x2 − x1 = · · · = ∆x.
We shall skip this material as no new maths is involved.
4.7 Bounding of Interpolation Error
Suppose we interpolate the function f(x) for a ≤ x ≤ b using a polynomial P (x) such that P (xj) = f(xj).
We would like to say something about f(x)− P (x) at other points, i.e. how large is the error? We can
derive a useful expression for this provided f(x) is (at least) n+ 1 times differentiable, by considering the
function
F (t) = f(t)− P (t)− [f(x)− P (x)] n∏
j=0
t− xj
x− xj , x ̸= xj .
· · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · ·
4-9
MATH256 2021–22 Numerical Methods (5 Parts)
Step 1: We prove that there exists a point ξ ∈ [a, b] such that F (n+1)(ξ) = 0.
There are at least n+ 2 points in the interval [a, b] at which F (t) = 0: t = xj for j = 0, . . . , n, and t = x.
Between each two roots there must be a turning point, so there are at least n+1 points at which F ′(t) = 0.
There must be at one least turning point of F ′(t) between each of these points, so there are at least n
points at which F ′′(t) = 0. Repeatedly applying this argument shows there exists at least one point ξ, say,
such that F (n+1)(ξ) = 0. Since the location of ξ depends on x, we write ξ = ξ(x).
Step 2: We obtain an expression for F (n+1)(t) and then set t = ξ(x).
We have
F (t) = f(t)− P (t)− f(x)− P (x)n∏
j=0
(x− xj)
[
tn+1 − (x0 + x1 + · · ·+ xn)tn + · · ·
]
so differentiating n+ 1 times yields
F (n+1)(t) = f (n+1)(t)− f(x)− P (x)n∏
j=0
(x− xj)
(n+ 1)!.
Finally, we set t = ξ(x) to eliminate F , and rearrange to obtain
f(x)− P (x) = f
(n+1) (ξ(x))
(n+ 1)!
n∏
j=0
(x− xj).
· · · · · · · · · · · · · · · · · · · · · · · · · · · ·
Remarks
• Although we derived this under the assumption that x is distinct from xj , both sides vanish if x = xj ,
so it is correct for all x ∈ [a, b].
• We don’t know ξ, but the ratio f (n+1) (ξ(x)) /(n + 1)! usually turns out to be relatively small for
large n (at least if f is a ‘reasonable’ function).
• For equally spaced nodes, the product term
Kn(x) =
n∏
j=0
(x− xj)
can take very large values, particularly near the end-points of the interpolation interval (see figure
4.1). This in turn can cause the approximation to break down if a large number of points are used
(see the demonstration in d_tab.mw).
4-10
MATH256 2021–22 Numerical Methods (5 Parts)
2.
x
2
4
6
8
10
−4 −2 0 2 4
n = 10
lo
g
1
0
|K
n
(x
)|
10
x
2
4
6
8
10
−4 −2 0 2 4
n = 20
lo
g
1
0
|K
n
(x
)|
Figure 4.1: Logarithmic plots of the product |Kn(x)| for x ∈ [−5, 5], using equally spaced nodes.
4.8 Chebyshev polynomials
To improve the accuracy of polynomial interpolation, we can choose the points xj in a more subtle way.
First we will consider interpolation on the interval [−1, 1]; later we can generalise the results to deal with
arbitrary intervals. The aim is then to choose the points xj such the product
Kn(x) =
n∏
j=0
(x− xj)
remains small for all x ∈ [−1, 1]. To achieve this, we consider the function Tn, which has the property that
Tn(cos θ) = cos(nθ), n = 0, 1, . . .
The idea here is that | cos(nθ)| ≤ 1, so we have a bound on the modulus, and cos(nθ) can be expanded to
show that Tn is a polynomial. These are called Chebyshev polynomials of the first kind.∗
4.8.1 Examples of first few Chebyshev polynomials
· · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · ·
The first two Chebyshev polynomials are given by
T0(cos θ) = 1 ⇒ T0(t) = 1
and
T1(cos θ) = cos θ ⇒ T1(t) = t.
At higher orders, we can use addition formulae. Thus
T2(cos θ) = cos(2θ)
= cos2 θ − sin2 θ
= 2 cos2 θ − 1,
∗Many different spellings of the name ‘Chebyshev’ appear in the literature. In general, one should look under ‘T’ (Tchebichef,
etc.) as well as ‘C’, particularly in indexes of older books.
4-11
MATH256 2021–22 Numerical Methods (5 Parts)
so that T2(t) = 2t2 − 1 and
T3(cos θ) = cos(3θ)
= cos(θ + 2θ)
= cos θ cos(2θ)− sin θ sin(2θ)
= cos θ
[
2 cos2 θ − 1]− 2 sin2 θ cos θ
= cos θ
[
2 cos2 θ − 1]− 2(1− cos2 θ) cos θ
= 4 cos3 θ − 3 cos θ,
meaning that T3(t) = 4t3 − 3t.
· · · · · · · · · · · · · · · · · · · · · · · · · · · ·
Evidently, the leading coefficient can’t be a fraction; in fact it turns out to be 2n−1 for each n (Problem 6.5).
We can make a polynomial with a smaller maximum modulus by dividing through:
T̂n(t) =
Tn(t)
2n−1 .
This is called a normalised Chebyshev polynomial. Experimenting with the plots in chebyshev_demo.mw
suggests there is no monic polynomial† Qn of degree n such that
max
−1≤t≤1
|Qn(t)| < max−1≤t≤1
∣∣T̂n(t)∣∣,
for n = 2 or n = 3. This is called the minimax property (minimising the maximum modulus). We now
need to find a general expression for Tn(t) and show that the minimax property holds for all n.
4.8.2 General expression for Tn(t)
Below we consider two ways of deriving a formula for computing Tn(t).
(1) By recursive formula — simple derivation for a recursive formula:
Since T0(cos θ) = cos0 θ, T0(t) = 1 and T1(cos θ) = cos1 θ, T1(t) = t,
we have T2(cos 2θ) = 2 cos2 θ − 1, T2(t) = 2t2 − 1.
For Tn+1(cos θ) = cos((n+ 1)θ), by keeping cos terms after expansion, we get
cos((n+ 1)θ) = cos((n− 1)θ + 2θ)
= cos(n− 1)θ cos(2θ)− sin(n− 1)θ sin(2θ)
= Tn−1T2 − 2 cos θ sin θ sin(n− 1)θ
= Tn−1T2 − 2T1[cos θ cos(n− 1)θ − cosnθ]
Tn+1(t) = Tn−1T2 − 2T1[T1Tn−1 − Tn]
= 2T1Tn + (T2 − 2T 21 )Tn−1
= 2tTn + (2t2 − 1− t2)Tn−1
= 2tTn − Tn−1.
(2) By binomial coefficients — a direct and non-recursive formula
· · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · ·
†Monic polynomials have leading coefficient 1: antn + an−1tn−1 + · · ·+ a0 is a monic polynomial of degree n if and only if
an = 1.
4-12
MATH256 2021–22 Numerical Methods (5 Parts)
Let θ ∈ R, and consider De Moivre’s formula
einθ = (cos θ + i sin θ)n
=
n∑
j=0
(
n
j
)
(i sin θ)j cosn−j θ,
where the binomial coefficient is given by (
n
j
)
= n!
j!(n− j)! .
Taking the real part, we obtain
cos(nθ) =
n∑
j=0
j even
(
n
j
)
(i sin θ)j cosn−j θ.
Next, we write j = 2p. The upper limit is then ⌊n/2⌋, where ⌊·⌋ means ‘round down’, so
cos(nθ) =
⌊n/2⌋∑
p=0
(
n
2p
)
(i sin θ)2p cosn−2p θ
=
⌊n/2⌋∑
p=0
(
n
2p
)
(cos2 θ − 1)p cosn−2p θ.
Since Tn(cos θ) = cos(nθ), we now have
Tn(t) =
⌊n/2⌋∑
p=0
(
n
2p
)(
t2 − 1
)p
tn−2p.
We can check this against our earlier results:
T0(t) =
0∑
p=0
(
0
2p
)(
t2 − 1
)p
t−2p = 1,
T1(t) =
0∑
p=0
(
1
2p
)(
t2 − 1
)p
t1−2p = t,
T2(t) =
1∑
p=0
(
2
2p
)(
t2 − 1
)p
t2−2p = 2t2 − 1,
T3(t) =
1∑
p=0
(
3
2p
)(
t2 − 1
)p
t3−2p = t3 + 3(t2 − 1)t = 4t3 − 3t.
· · · · · · · · · · · · · · · · · · · · · · · · · · · ·
4.8.3 Maxima and minima of Chebyshev polynomials
· · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · ·
4-13
MATH256 2021–22 Numerical Methods (5 Parts)
For t ∈ [−1, 1] (the range for cos θ), we know that |Tn(t)| ≤ 1, because
|Tn(cos θ)| = | cos(nθ)| ≤ 1.
The extreme values ±1 are achieved by setting
θ = jπ
n
, j ∈ Z,
which gives t = χj , where
χj = cos
(
jπ
n
)
.
Evidently, χ0 = χ2n = 1 after which the points χj are repeated. Also,
χ2n−j = cos
((2n− j)π
n
)
= cos(2π) cos
(
jπ
n
)
= χj ,
so there are n+ 1 distinct values for χj , obtained by taking j = 0, . . . , n. Two of these are the end points
χ0 = 1 and χn = −1; there are n− 1 stationary points in between. Since Tn(t) is a polynomial of degree
n, there can be no other stationary points. Note that
Tn(χj) = Tn
(
cos(jπ/n)
)
= cos
(
n× jπ
n
)
= (−1)j ,
so χ1, χ3, . . . are minima, whereas χ0, χ2, . . . are maxima.
· · · · · · · · · · · · · · · · · · · · · · · · · · · ·
4.8.4 Roots of Chebyshev polynomials
To find the points where Tn(t) = 0, we use the formula Tn(cos θ) = cos(nθ) again, this time choosing θ so
that cos(nθ) = 0.
· · · · · · · · · · · · · · · · · · · · · · · · · · · · · · ·
Thus,
nθ =
(
j + 12
)
π ⇒ θ =
(
j + 12
) π
n
, j ∈ Z,
which gives us
tj = cos
((
j + 12
) π
n
)
.
Again, these repeat starting at j = 2n, but also
t2n−1−j = cos
((
2n− j − 1 + 12
) π
n
)
= cos
((
−j − 12
) π
n
)
= tj ,
meaning only the points t0, t1, . . . , tn−1 are distinct. Since Tn(t) is a polynomial of degree n, it follows
that these are all of the roots.
· · · · · · · · · · · · · · · · · · · · · · · · · · · ·
4-14
MATH256 2021–22 Numerical Methods (5 Parts)
4.8.5 Example of the normalised Chebyshev polynomial T̂n(t) in Maple
Set a value for n and then use the command
plot( ChebyshevT( n , t ) / 2^(n-1) , t = -1.1 .. 1.1 )
to make a plot of a normalised Chebyshev polynomial in Maple. Observe the pattern of roots and stationary
points.
4.8.6 Theory of Minimax Property of T̂n(t)
Amongst all the monic polynomials of degree n (with real coefficients), the normalised Chebyshev polynomial
T̂n(t) =
Tn(t)
2n−1
has the smallest maximum modulus for t ∈ [−1, 1].
Proof
· · · · · · · · · · · · · · · · · · · · · · · · · · · · · · ·
Clearly there is only one monic polynomial of degree zero, so we may assume without loss that n ≥ 1.
Suppose there exists a monic polynomial Q(t) of degree n such that
max
−1≤t≤1
|Q(t)| < max
−1≤t≤1
∣∣T̂n(t)∣∣ (∗)
and consider the difference
R(t) = T̂n(t)−Q(t).
Since T̂n(t) and Q(t) are monic, the terms in tn cancel, meaning R(t) is a polynomial whose degree is at
most n− 1. At the maxima of the Chebyshev polynomial (section 4.8.3), we have
T̂n(χ2j) = max−1≤t≤1 T̂n(t) =
1
2n−1
and at the minima
T̂n(χ2j−1) = min−1≤t≤1 T̂n(t) = −
1
2n−1 .
Hence, R(χj) > 0 when j is even, and R(χj) < 0 when j is odd. Between each pair of these values, there
must be a point at which R(t) = 0. However, there are n + 1 points χj , so there must be n points at
which R(t) = 0. This is impossible unless R(t) = 0 for all t, in which case Q(t) = T̂n(t), but this is also
impossible in view of (∗). Therefore we have a contradiction, so our assumption (the existence of Q(t)
satisfying (∗)) is false.
· · · · · · · · · · · · · · · · · · · · · · · · · · · ·
4-15
MATH256 2021–22 Numerical Methods (5 Parts)
4.8.7 Interpolation using Chebyshev nodes
Recall that the error in interpolating f(x) using a polynomial P (x) of degree n is given by
f(x)− P (x) = f
(n+1) (ξ(x))
(n+ 1)! Kn(x),
where
Kn(x) =
n∏
j=0
(x− xj)
(see section 4.7). We will minimise the maximum value of |Kn| by creating a Chebyshev polynomial on the
right-hand side. Since there are n+ 1 nodes xj , we will use the roots of Tn+1(t). These are given by
tj = cos
((
j + 12
) π
n+ 1
)
, j = 0, . . . , n.
Since cos θ is a decreasing function for 0 ≤ θ ≤ π, these are arranged in descending order. Therefore we
map the interval [−1, 1] (for t) to [a, b] (for x) using the linear transformation
x = a− b2 t+
a+ b
2 .
In particular, t = −1 and t = 1 correspond to x = b and x = a, respectively. The Chebyshev nodes for
[a, b] are given by
xj =
a− b
2 tj +
a+ b
2 , j = 0, 1, . . . n.
Once we have calculated the nodes, we can construct the Newton polynomial using a table of divided
differences. Let us now consider the error in this approximation.
· · · · · · · · · · · · · · · · · · · · · · · · · · · · · · ·
Using Chebyshev nodes and writing Kn in terms of t, we obtain
Kn(x) =
n∏
j=0
(
a− b
2 t+
a+ b
2 −
a− b
2 tj −
a+ b
2
)
=
n∏
j=0
(
a− b
2 (t− tj)
)
=
(
a− b
2
)n+1 n∏
j=0
(t− tj)
=
(
a− b
2
)n+1
T̂n+1(t),
because the points tj are the roots of T̂n+1(t). Finally, since |T̂n+1(t)| ≤ 2−n, we arrive at the final result
|f(x)− P (x)| ≤
∣∣f (n+1) (ξ(x))∣∣
(n+ 1)! × 2
(
b− a
4
)n+1
.
· · · · · · · · · · · · · · · · · · · · · · · · · · · ·
4-16
MATH256 2021–22 Numerical Methods (5 Parts)
Remarks
• By switching from equally spaced nodes to Chebyshev nodes, we have replaced the factor ∏nj=0(x−xj)
(about which we know very little) with 2((a− b)/4)n+1 which can easily be calculated.
• The error formula now shows explicitly that increasing b − a makes accurate interpolation more
difficult, but increasing n yields a better result (provided f (n+1)(ξ(x)) behaves itself).
• Chebyshev nodes will generally produce much better results than equally spaced nodes.
4.8.8 Maple Examples
See d_tab.mw.
4.9 Splines
If we need to fit a curve through experimental data, we may not have the option of choosing the distribution
of the points, so we can no longer use Chebyshev nodes. Instead, we can use a piecewise fit, with a function
called a spline between each pair of points. This avoids the need for polynomials of very high degree.
Denote the spline from xj to xj+1 by Sj(x), so that
Sj(xj) = yj and Sj(xj+1) = yj+1, j = 0, 1, . . . , n− 1;
see figure 4.2. For a smooth curve, we need continuous gradient and curvature‡, so the first and second
derivatives at the end of each spline must match those at the start of the next, i.e.
S′j(xj+1) = S′j+1(xj+1) and S′′j (xj+1) = S′′j+1(xj+1), j = 0, 1, . . . , n− 2.
This gives us four conditions for each spline (two data values and two matching derivatives), except the
last, which has only two. The simplest function that can satisfy four such conditions is a cubic polynomial.
Before trying to determine these, let us consider the numbers of equations and unknowns.
• We have n cubic polynomials so there are 4n coefficients to find.
• Fitting a spline through a data point generates an equation that does not involve coefficients from
any other spline. We can use this to eliminate one unknown. Each spline passes through two data
points, so 2n unknowns are eliminated and 2n remain.
• Matching first and second derivatives at interior points relates each spline to an adjacent spline. Since
there are n− 1 interior points, this generates a system of 2(n− 1) equations.
• If we make two assumptions, we will have 2(n− 1) equations for 2(n− 1) unknowns. It turns out
that we can eliminate half of the coefficients from this system.
Next we find equations for the splines. We start by writing
Sj(x) = aj(x− xj)3 + bj(x− xj)2 + cj(x− xj) + dj
Clearly Sj(xj) = yj implies dj = yj so there are only 3n unknowns left: {aj , bj , cj}.
· · · · · · · · · · · · · · · · · · · · · · · · · · · · · · ·
‡The curvature of a function s(x) is given by κ(x) = s′′(x)
(
1 + [s′(x)]2
)−3/2.
4-17
MATH256 2021–22 Numerical Methods (5 Parts)
xx0 x1 x2 x3
x
y0
y1
y2
y3
S0
S1
S2
Matching
derivatives here
No matching
derivatives here
Figure 4.2: A cubic spline interpolation through the points (xj , yj), j = 0, 1, 2, 3.
0th Derivatives. Evaluating at x = xj+1 and writing hj = xj+1 − xj , we obtain
Sj(xj+1) = yj+1 = ajh3j + bjh2j + cjhj + yj ,
and using the definition of a divided difference this becomes
ajh
2
j + bjhj + cj = y[xj , xj+1], j = 0, 1, . . . , n− 1. (i)
Note Sj(xj) = yj was satisfied and Sj(xj+1) = yj+1 is also satisfied for j = 0, 1, . . . , n− 1.
1st Derivatives. Now consider the derivatives of the splines:
S′j(x) = 3aj(x− xj)2 + 2bj(x− xj) + cj
and
S′j+1(x) = 3aj+1(x− xj+1)2 + 2bj+1(x− xj+1) + cj+1.
Matching S′j(x) = S′j+1(x) at x = xj+1, we find that
3ajh2j + 2bjhj + cj = cj+1, j = 0, 1, . . . , n− 1 (ii)
where j = 0, 1, . . . , n− 2 is expected but we add an extra match at xn to allow j = n− 1 i.e. to allow Sn
to exist (since Sn(x) sits beyond xn i.e. it does not exist).
2nd Derivatives. For the second derivative,
S′′j (x) = 6aj(x− xj) + 2bj
and
S′′j+1(x) = 6aj+1(x− xj+1) + 2bj+1.
Equating these at x = xj+1 yields
3ajhj + bj = bj+1, j = 0, 1, . . . , n− 1. (iii)
4-18
MATH256 2021–22 Numerical Methods (5 Parts)
Here as before, j = 0, 1, . . . , n− 2 is expected but we add an extra match at xn to allow j = n− 1 (since
Sn(x) does not exist).
To summarise, from (i),(ii),(iii), there are 3n equations and the total number of unknowns is 3n+ 2.
Here two of the unknowns bn, cn are extra, each having its own equation in (ii) and (iii) respectively.
• Theoretically we could fix any two of these 3n + 2 unknowns to have a square system and hence a
hopefully unique solution by solving an 3n × 3n system. Below we first simplify these 3n equations by
as much elimination as possible. Then we decide which 2 unknowns should be fixed to leave us a most
convenient system.
We can use aj from (iii) to eliminate aj from (i)-(ii) and then cj next. We are then left with
(i)⇒ bj+1 + 2bj3 hj + cj = y[xj , xj+1] ⇒ cj = y[xj , xj+1]−
bj+1 + 2bj
3 hj ,
⇒ cj+1 = y[xj+1, xj+2]− bj+2 + 2bj+13 hj ....................................................... (∗)
(ii)⇒ (bj+1 + bj)hj + cj = cj+1........................................................................... (†)
Note that the index j = 0, 1, . . . , n− 2 holds in (*) for cj+1 while it was j = 0, 1, . . . , n− 1 with cj . Hence,
(bj+1 + bj)hj + y[xj , xj+1]− bj+1 + 2bj3 hj︸ ︷︷ ︸
cj
= y[xj+1, xj+2]− bj+2 + 2bj+13 hj+1︸ ︷︷ ︸
cj+1
is only valid for j = 0, 1, . . . , n− 2. Gathering all bj terms, we obtain
(†)⇒ hj+1bj+2 + 2(hj+1 + hj)bj+1 + hjbj = 3
(
y[xj+1, xj+2]− y[xj , xj+1]
)
.
· · · · · · · · · · · · · · · · · · · · · · · · · · · ·
The range of indices j = 0, 1, . . . , n−2 implies that we have n−1 equations for n+1 unknowns b0, b1, . . . , bn
or we have two extra unknowns. This is caused by the fact that we do not impose conditions on the first
and last splines involving end points. To get around this, we must make a ‘closure assumption’ at end
points, the simplest of which is (mainly to maintain a nicely ‘structured’ matrix)
b0 = bn = 0.
The equations for bj can then be written in matrix-vector form (note: A is symmetric, diagonally dominant)
2(h0 + h1) h1
h1 2(h1 + h2) h2
h2 2(h2 + h3) h3
. . . . . . . . .
hn−3 2(hn−3 + hn−2) hn−2
hn−2 2(hn−2 + hn−1)
b1
b2
b3
...
bn−2
bn−1
= 3
y[x1, x2]− y[x0, x1]
y[x2, x3]− y[x1, x2]
y[x3, x4]− y[x2, x3]
...
y[xn−2, xn−1]− y[xn−3, xn−2]
y[xn−1, xn]− y[xn−2, xn−1]
.
4-19
MATH256 2021–22 Numerical Methods (5 Parts)
If you think about it, setting b1 = b2 = 0 would ‘ruin’ the nice tri-diagonal structure of matrix A.
After solving this for bj ’s, we can construct the splines (for j = 0, 1, . . . , n− 1) using
aj =
bj+1 − bj
3hj
and cj = y[xj , xj+1]− bj+1 + 2bj3 hj .
Remarks
• The coefficient bn, though used in computing an−1, cn−1, doesn’t appear directly in the expressions
for the splines, but it does have a meaning:
S′′n−1(xn) = 6an−1hn−1 + 2bn−1 = 2bn,
so its value is half the second derivative at the right end of the curve.
• Similarly, S′′0 (x0) = 2b0, so the closure assumption b0 = bn = 0 causes the second derivative to vanish
at both ends. Other closure assumptions are possible (see Problem 7.3).
• Now is the time to review how the error analysis from §4.7 (Bounding of Interpolation Error) applies
to Sj(x). For some ξ ∈ [a, b], the upper bound will involve |f (4)(ξ)| which is expected to be MUCH
smaller than |f (n+1)(ξ)| since n can be arbitrarily large.
• Final thought: Here Sj(x) is a degree-3 polynomial, smoothly pieced together by continuity conditions
at interior nodes. Hence if we want a difference degree (e.g. quadratic of degree 2 or quartic of
degree 4, or quintic of degree 5) for Sj(x), a similar procedure can be devised to construct them.
However, linear polynomials for Sj(x) cannot have continuity at nodes – why?
Examples
See splines.mw.
4-20
