MATH256 Numerical Methods
Numerical Methods
Hello, dear friend, you can consult us at any time if you have any questions, add WeChat: THEend8_
Numerical Methods MATH256
Contents
* 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 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.