MTH618 Wk 5 Lect 1

OPTIMISATION METHODS
The materials in this topic is particularly useful in modeling large-scale real-world
problems.
Problems, such as optimizing production plans for different industries (microchips,
pharmaceuticals, cars, aluminum, steel, chemicals), optimizing usage of transportation
systems (usage of runways in airports, tracks of subways), efficiency in running of power
plants, optimal shipping (delivery services, shipping of containers, shipping goods from
factories to warehouses and from warehouses to stores), designing optimal financial
portfolios and others are all applications of the optimisation methods. But the size
of the problem sometimes requires the use of optimization software. More recently,
environmental concerns have put new aspects into the picture, where an important
concern, added to these problems, is the minimization of environmental impact. The
main task becomes to model these problems correctly.

Interpolation
This part falls under Numerical Analysis which which deals with methods that
construct new unknown function values from known function values. The knowledge
gained here is applied to spline interpolation and is also useful for understanding
numericalintegration and differentiation.
Given ordered pairs of points

(x0 , y0 ), (x1 , y1 ), (x2 , y2 ), ··· (xn , yn ) ,

we need to find a function f such that

f (x0 ) = y0 , f (x1 ) = y1 , f (x2 ) = y2 , ··· f (xn ) = yn .

In other words, a function which touches all the given points. This process is called
interpolation.
The standard idea in interpolation now is to find a polynomial pn (x) of degree n (or
less) such that
f = pn .
We call pn an interpolation polynomial and x0 , . . . , xn nodes. If f is another type
of mathematical function, we call pn an polynomial approximation. If we use pn to
get values of f within x0 ≤ x ≤ xn , it is known as interpolation and if outside this
interval x0 ≤ x ≤ xn , this is called extrapolation. Note that for the given data, pn is
unique.

1
Lagrange Interpolation
Linear Lagrange Interpolation
• This is interpolation by the straight line through (x0 , y0 ), (x1 , y1 ).
• The linear Lagrange polynomial is
p1 = L 0 y 0 + L 1 y 1
with L0 and L1 to be linear polynomials (straight lines) which satisfy
L0 (x0 ) = 1 L1 (x0 ) = 0
L0 (x1 ) = 0 L1 (x1 ) = 1.

• Solving (Exercise!) gives

x−x1 x−x0
L0 (x) = x0 −x1
, L1 (x) = x1 −x0
.

This gives the linear Lagrange polynomial

x − x1 x − x0
p1 (x) = L0 (x) · y0 + L1 (x) · y1 = · y0 + · y1
x0 − x1 x1 − x0

Example. Compute the 4DP -value of ln 9.2 from ln 9.0 = 2.1972, ln 9.5 = 2.2513 by
linear Lagrange interpolation and determine the error, using ln 9.2 = 2.2192 (4D).
Solution. Note that
x0 = 9.0, x1 = 9.5, y0 = ln 9.0, y1 = ln 9.5 .
So
x − 9.5 x − 9.5
L0 (x) = = = −2.0(x − 9.5)
9.0 − 9.5 −0.5
x − 9.0
L1 (x) = = 2.0(x − 9.0).
9.5 − 9.0
Then L0 (9.2) = 0.6 and L1 (9.2) = 0.4. Hence
ln 9.2 ≈ p1 (9.2) = L0 (9.2) · y0 + L1 (9.2) · y1
= (0.6)(2.1972) + (0.4)(2.2513) = 2.2188.
The error is
 = a − ã = 2.2192 − 2.2188 = 0.0004
which shows that linear interpolation here does not reach 4D-accuracy; it would suffice
for 3D-accuracy.

2
Quadratic Lagrange Interpolation
• Quadratic interpolation involves extending the idea of linear interpolation above.
• Interpolating points (x0 , y0 ), (x1 , y1 ), (x2 , y2 ) using a 2nd-degree polynomial (quadratic
equation) is
p2 (x) = L0 (x) · y0 + L1 (x) · y1 + L2 (x) · y2
with L0 (x0 ) = 1, L1 (x1 ) = 1, L2 (x2 ) = 1 and otherwise 0.
• Denote
lk (x) = (x − x0 )(x − x1 ) · · · (x − xn )
with (x − xk ) term missing, then
l0 (x) (x − x1 )(x − x2 )
L0 (x) = =
l0 (x0 ) (x0 − x1 )(x0 − x2 )
l1 (x) (x − x0 )(x − x2 )
L1 (x) = =
l1 (x1 ) (x1 − x0 )(x1 − x2 )
l2 (x) (x − x0 )(x − x1 )
L2 (x) = = .
l2 (x0 ) (x2 − x0 )(x2 − x1 )

Example. Compute the 4D-value of ln 9.2 from

ln 9.0 = 2.1972, ln 9.5 = 2.2513, ln 11.0 = 2.3979 .

by quadratic Lagrange interpolation.

Solution. Calculating
(x − 9.5)(x − 11.0)
L0 (x) = = x2 − 20.5x + 104.5
(9.0 − 9.5)(9.0 − 11.0)
(x − 9.0)(x − 11.0) 1
L1 (x) = =− (x2 − 20x + 99)
(9.5 − 9.0)(9.5 − 11.0) 0.75
(x − 9.0)(x − 9.5) 1
L2 (x) = = (x2 − 18.5x + 85.5).
(11.0 − 9.0)(11.0 − 9.5) 3
Then
L0 (9.2) = 0.5400, L1 (9.2) = 0.4800, L2 (9.2) = −0.0200
and
ln 9.2 ≈ p2 (9.2) = (0.5400)(2.1972) + (0.4800)(2.2513) + (−0.0200)(2.3979)
= 2.2192.

3
General Lagrange Interpolation Polynomial
n n
X X lk (x)
f (x) ≈ pn (x) = Lk (x)yk = · yk
k=0 k=0
lk (xk )

where Lk (xk ) = 1 and Lk is 0 at the other nodes and

l0 (x) = (x − x1 )(x − x2 ) · · · (x − xn )
lk (x) = (x − x0 ) · · · (x − xk−1 )(x − xk+1 ) · · · (x − xn )
ln (x) = (x − x0 )(x − x1 ) · · · (x − xn−1 ).
The error estimate is then

f (n+1) (t)
n (x) = f (x) − pn (x) = (x − x0 )(x − x1 ) · · · (x − xn )
(n + 1)!

For this course, we will only utilise Linear and Quadratic Lagrange Interpolation.

4
Newton’s Divided Difference Interpolation
For the given data (x0 , y0 ), · · · , (xn , yn ), we find pn (x) but in a different form than
Lagrange’s method. Practically more important are Newton’s forms of pn (x), which we
can use for solving ODEs.
Let pn−1 (x) be the (n − 1)st Newton polynomial (which will be found later), thus
pn−1 satisfies

pn−1 (x0 ) = y0 , pn−1 (x1 ) = y1 , ··· pn−1 (xn ) = yn .

Let use write
pn (x) = pn−1 (x) + gn (x),
so
gn (x) = pn (x) − pn−1 (x).
Here gn (x) has to be determined so that

pn (x0 ) = y0 , pn (x1 ) = y1 , ··· pn (xn ) = yn .

Note that
gn (x) = an (x − x0 )(x − x1 ) · · · (x − xn ).
So we need to determine the constant an . Setting x = xn and substituting,
yn − pn−1 (xn )
an = .
(x − x0 )(x − x1 ) · · · (x − xn )
Writing ak instead of an , we obtain the kth divided difference, recursively defined
below.
y1 − y0
a1 = f [x0 , x1 ] =
x1 − x0
f [x1 , x2 ] − f [x0 , x1 ]
a2 = f [x0 , x1 , x2 ] =
x2 − x 0
..
.
f [x1 , . . . , xk ] − f [x0 , . . . , xk−1 ]
ak = f [x0 , . . . , xk ] =
xk − x0
• If n = 1, we have a1 as defined above and gives the Newton interpolation polyno-
mial of 1st degree
p1 (x) = y0 + (x − x0 ) · f [x0 , x1 ]

• If n = 2, then with p1 , we find a2 as defined above, then second Newton polynomial

p2 (x) = f0 + (x − x0 ) · f [x0 , x1 ] + (x − x0 )(x − x1 ) · f [x0 , x1 , x2 ]

• For n = k, we obtain

pk (x) = pk−1 (x) + (x − x0 )(x − x1 ) · · · (x − xk−1 ) · f [x0 , x1 , . . . , xk ]

• With p0 (x) = y0 , by repeated application with k = 1, . . . , n, this finally gives

Newton’s divided difference interpolation formula

f (x) ≈ pn (x)

5
Example. Compute f (9.2) where f (x) = ln x, from the values shown in the following
table.
xj yj = f (xj )
x0 = 8.0 y0 = 2.079442
x1 = 9.0 y1 = 2.197225
x2 = 9.5 y2 = 2.251292
x3 = 11.0 y3 = 2.397895
Solution. The computations are:
y1 − y0 2.197225 − 2.079442
f [x0 , x1 ] = = = 0.117783
x1 − x0 9.0 − 8.0
y2 − y1
f [x1 , x2 ] = = 0.108134
x2 − x1
y3 − y2
f [x2 , x3 ] = = 0.097735
x3 − x2
Then find
f [x1 , x2 ] − f [x1 , x0 ]
f [x0 , x1 , x2 ] = = −0.006433
x2 − x0
f [x2 , x3 ] − f [x1 , x2 ]
f [x1 , x2 , x3 ] = = −0.005200
x3 − x1
Then
f [x1 , x2 , x3 ] − f [x0 , x1 , x2 ]
f [x0 , x1 , x2 , x3 ] = = 0.000411
x3 − x0
Tabulating the values, we see the calculation pattern forming:
xj yj = f (xj ) f [xi , xi+1 ] f [xi , xi+1 , xi+2 ] f [xi , xi+1 , xi+2 , xi+3 ]
x0 = 8.0 y0 = 2.079442
0.117783
x1 = 9.0 y1 = 2.197225 −0.006433
0.108134 0.000411
x2 = 9.5 y2 = 2.251292 −0.005200
0.097735
x3 = 11.0 y3 = 2.397895
Note that
y0 = 2.079442
a1 = f [x0 , x1 ] = 0.117783
a2 = f [x0 , x1 , x2 ] = −0.006433
a3 = f [x0 , x1 , x2 , x3 ] = 0.000411
then
f (x) ≈ p3 (x)
= 2.079442 + 0.117783(x − 8.0) − 0.006433(x − 8.0)(x − 9.0)
+0.000411(x − 8.0)(x − 9.0)(x − 9.5).
At x = 9.2,
f (9.2) ≈ p3 (9.2) = 2.219208.
The exact value to 6DP is f (9.2) = ln 9.2 = 2.219203. Note that we can nicely see how
the accuracy increases from term to term:
p1 (9.2) = 2.220782, p2 (9.2) = 2.219238, p3 (9.2) = 2.219208.

6
 MTH618 Wk 5 Lect 2

Spline Interpolation
Given data (function values, points in the xy-plane) (x0 , y0 ), (x1 , y1 ), ..., (xn , yn ) can
be interpolated by a polynomial Pn (x) of degree n or less so that the curve of Pn (x)
passes through these n + 1 points (xi , yi ).
Note that if n is large, then there may be trouble with the curve; Pn (x) may tend to
oscillate for x between the nodes x0 , . . . , xn . Hence, we must be prepared for numeric
instability.
A famous example by C. Runge (1856-1927) for which the maximum error even
approaches ∞ as n → ∞ (with nodes kept equidistant and their number increased) is
given in the first diagram. The second illustration shows the increase of the oscillation
with n for some other function that is piecewise linear.

These undesirable oscillations are avoided by the method of splines initiated by I. J.

Schoenberg in 1946. This method is widely used in practice. It also laid the foundation
for much of modern Computer Aided Design (CAD).
Instead of using a single high-degree polynomial Pn over the entire interval a ≤
x ≤ b in which the nodes lie, we use n number of low degree polynomials e.g., cubic,
polynomials
q0 (x), q1 (x), ··· , qn−1 (x),

1
one over each subinterval between adjacent nodes. So
x0 to x1 =⇒ use q0 (x)
x1 to x2 =⇒ use q1 (x)
.. .. ..
. . .
xn−1 to xn =⇒ use qn−1 (x).
From this, we construct an interpolation function



 qo (x), x 0 ≤ x ≤ x1
q1 (x), x1 ≤ x ≤ x2

g(x) = . .

 .. ..


q (x), x
n−1 ≤x≤x n−1 n

called a spline, by fitting these polynomials together into a single continuous curve
passing through the data points, that is,
g(x0 ) = y0 , g(x1 ) = y1 , ··· , g(xn ) = yn .
Thus spline interpolation is piecewise polynomial interpolation.
The simplest qj ’s would be linear polynomials. However, the curve of a piecewise
linear continuous function has corners and would be of little interest in general—think
of designing the body of a car or a ship.

Cubic Splines
We shall consider cubic splines because these are the most important ones in appli-
cations. By definition, a cubic spline g(x) interpolating data (x0 , y0 ), . . . , (xn , yn ) is a
continuous function on the interval x0 ≤ x ≤ xn that has continuous first and second
derivatives and satisfies the interpolation conditions
g(x0 ) = y0 , g(x1 ) = y1 , ··· , g(xn ) = yn .

Furthermore, between adjacent nodes, g(x) is given by a polynomial qj (x) of degree 3

or less. In addition to the above conditions, we also require that
g 0 (x0 ) = k0 , g 0 (xn ) = kn

(given tangent directions of g(x) at the endpoints of the interval x0 ≤ x ≤ xn .

Theorem. (Existence and Uniqueness of Cubic Splines) Let (x0 , y0 ), (x1 , y1 ), · · · , (xn , yn )
with given (arbitrarily spaced) xj and given yj = f (xj ), j = 0, 1, · · · , n. Let k0 and kn
be any given numbers. Then there is one and only one cubic spline g(x) corresponding
to the nodes and satisfying
g(x0 ) = y0 , g(x1 ) = y1 , ··· , g(xn ) = yn
and
g 0 (x0 ) = k0 , g 0 (xn ) = kn .
Remark. Note that if we are trying to interpolate a function f (x) by constructing g(x),
then the conditions
g 0 (x0 ) = k0 , g 0 (xn ) = kn

2
become
g 0 (x0 ) = f 0 (x0 ), g 0 (xn ) = f 0 (x)n
which are known as clamped conditions. Other conditions of interest are the free or
natural conditions
g 00 (x0 ) = 0, g 00 (xn ) = 0
(geometrically: zero curvature at the ends).

Determination of the Splines

Let k0 and kn be given. Obtain k1 , · · · , kn by solving the linear system of equations
given by

cj−1 · kj−1 + 2(cj−1 + cj ) · kj+1 + cj · kj+1 = 3 c2j−1 ∆fj + c2j ∆fj+1

 

where

1 1
cj = =
hj xj+1 − xj
1 1
cj−1 = =
hj−1 xj − xj−1
hj = xj+1 − xj
and j = 1, . . . , n − 1.
Recall that the spline g(x) to be found consists of n cubic polynomials q0 , . . . , qn−1 .
We write these polynomial in the form

qj (x) = aj0 + aj1 · (x − xj ) + aj2 · (x − xj )2 + aj3 · (x − xj )3

where j = 0, . . . , n − 1. Using Taylor’s formula, we obtain

aj0 = qj (xj ) = yj
aj1 = qj0 (xj ) = kj
1 3 1
aj2 = qj00 (xj ) = 2 · (yj+1 − yj ) − · (kj+1 + 2kj )
2 hj hj
1 2 1
aj3 = qj000 (xj ) = 3 · (yj − yj+1 ) + 2 · (kj+1 − kj )
6 hj hj

Derivations:
• For aj0 , note that from the polynomial formula,
qj (xj ) = aj0 + aj1 · (xj − xj ) + aj2 · (xj − xj )2 + aj3 · (xj − xj )3
= aj0 .

• For aj1 , we use the 1st derivative of qj (x),

qj0 (x) = aj1 + 2aj2 (x − xj ) + 3aj3 (x − xj )2
qj0 (xj ) = aj1 + 2aj2 (xj − xj ) + 3aj3 (xj − xj )2
qj0 (xj ) = aj1 .

3
 • For aj2 , we use the 2nd derivative of qj (x),

qj00 (x) = 2aj2 + 6aj3 (x − xj )

qj00 (xj ) = 2aj2 + 6aj3 (xj − xj )
qj00 (xj ) = 2aj2
and
qj00 (xj ) = −6c2j · f (xj ) + 6c2j · f (xj+1 ) − 4cj kj − 2cj kj+1
where yj = f (xj ) from the proof to obtain
1 1
aj2 = qj00 (x) = −6c2j · f (xj ) + 6c2j · f (xj+1 ) − 4cj kj − 2cj kj+1

2 2
= −3c2j · f (xj ) + 3c2j · f (xj+1 ) − 2cj kj − cj kj+1
3 1
= 2 (−f (xj ) + f (xj+1 )) − (2kj + kj+1 ) .
hj hj

• For aj3 , we use the 2rd derivative of qj (x) again,

qj00 (x) = 2aj2 + 6aj3 (x − xj )

qj00 (xj+1 ) = 2aj2 + 6aj3 (xj+1 − xj )
= 2aj2 + 6aj3 hj
and
qj00 (xj+1 ) = 6c2j · f (xj ) − 6c2j · f (xj+1 ) + 2cj kj + 4cj kj+1
where yj = f (xj ) from the proof to obtain
qj (xj+1 ) = qj (xj+1 )
2aj2 + 6aj3 hj = 6c2j · f (xj ) − 6c2j · f (xj+1 ) + 2cj kj + 4cj kj+1 .

Since 2aj2 = 6c2j (−f (xj ) + f (xj+1 )) − cj (4kj + 2kj+1 ), then subtracting 2aj2 out,
we obtain

[2aj2 + 6aj3 hj ] − 2aj2 = 6c2j · f (xj ) − 6cj · f (xj+1 ) + 2cj kj + 4cj kj+1 − 2aj2 .
 

6aj3 hj = 12c2j · f (xj ) − 12c2j · f (xj+1 ) + 6cj kj + 6cj kj+1

12 12 6 6
6aj3 hj = 2
· f (xj ) − 2 · f (xj+1 ) + kj + kj+1
hj hj hj hj
2 2 6 6
aj3 = 3 · f (xj ) − 3 · f (xj+1 ) + 2 kj + 2 kj+1
hj hj hj hj
2 6
= 3 · [f (xj ) − ·f (xj+1 )] + 2 [kj + kj+1 ] .
hj hj

1
Remark. Note that for equidistant nodes of distance hj = h, we can write cj = c = h
and this simplifies linear system to solve to

3
kj−1 + 4kj + kj+1 = (yj+1 − yj−1 )
h
where j = 1, . . . , n − 1.

4
Example. Spline Interpolation. Equidistant Nodes. Interpolate f (x) = x4 on
the interval −1 ≤ x ≤ 1 by the cubic spline g(x) corresponding to the nodes x0 = −1,
x1 = 0, x2 = -1 and satisfying the clamped conditions g0(−1) = f 0(−1) and g0(1) = f 0(1).

Solution. In the standard notation, the given data are

f (−1) = 1 =⇒ x0 = −1, y0 = 1
f (0) = 0 =⇒ x1 = 0, y1 = 0
f (1) = 1 =⇒ x2 = 1, y2 = 1.

Note that h = 1 and n = 2, so our spline consists of n = 2 polynomials

qj (x) = aj0 + aj1 · (x − xj ) + aj2 · (x − xj )2 + aj3 · (x − xj )3
q0 (x) = a00 + a01 (x − (−1)) + a02 (x − (−1))2 + a03 (x − (−1))3
= a00 + a01 (x + 1) + a02 (x + 1)2 + a03 (x + 1)3
q1 (x) = a10 + a11 x + a12 x2 + a13 x3 .

Since xj are equidistant, we find kj from

3
kj−1 + 4kj + kj+1 = (yj+1 − yj−1 ).
h
Since n = 2, the system of linear equations is a single equation (with j = 1 and h = 1)
k0 + 4k1 + k2 = 3(y2 − y0 ).
From the clamped conditions

k0 = f 0 (x0 ) = 4x30 = 4(−1)3 = −4

k2 = f 0 (x2 ) = 4x32 = 4(1)3 = 4

then
(−4) + 4k1 + (4) = 3(1 − 1)
k1 = 0.

Then we can now obtain the coefficients of q0 , namely,

a00 = y0 = 1
a01 = k0 = −4
3 1
a02 = 2 (y1 − y0 ) − (k1 + 2k0 ) = 3(0 − 1) − (0 − 8) = 5
1 1
2 1
a03 = 3 (y0 − y1 ) + (k1 − k0 ) = 2(1 − 0) + (0 − (−4)) = −2.
1 12
Similarly, for the coefficients of q1 , we obtain

a10 = y1 = 0
a11 = k1 = 0
a12 = 3(y2 − y1 ) − (k2 + 2k1 ) = 3(1 − 0) − (4 + 0) = −1
a13 = 2(y1 − y2 ) + (k2 − k1 ) = 2(0 − 1) + (4 − 0) = 2.

5
So listing out the polynomials

q0 (x) = 1 − 4(x + 1) + 5(x + 1)2 − 2(x + 1)3

= −x2 − 2x3
q1 (x) = (0) + (0)x + (−1)x2 + (2)x3
= −x2 + 2x3 .
Then (
−x2 − 2x3 , −1 ≤ x ≤ 0
g(x) =
−x2 + 2x3 , 0 ≤ x ≤ 1.

Example. Natural Spline. Arbitrarily Spaced Nodes. Find a spline approximation and
a polynomial approximation for the curve of the cross section of the circular shaped
Shrine of the Book in Jerusalem shown below.

Solution. Thirteen points, about equally distributed along the contour give the
data:
xj −5.8 −5.0 −4.0 −2.5 −1.5 −0.8 0 0.8 1.5 2.5 4.0 5.0 5.8
yj 0 1.5 1.8 2.2 2.7 3.5 3.9 3.5 2.7 2.2 1.8 1.5 0
• The polynomial approximation for the curve is
P12 (x) = 3.9 − 0.65083x2 + 0.033858x4 + 0.011041x6 − 0.0014010x8
+0.000055595x10 − 0.00000071867x12
which is useless as shown in the figure above.

6
 • The coefficients of the cubic spline approximation are given below.
j x-interval aj0 aj1 aj2 aj3
0 0.0 − 0.8 3.9 0.00 −0.61 −0.015
1 0.8 − 1.5 3.5 −1.01 −0.65 0.66
2 1.5 − 2.5 2.7 −0.95 0.73 −0.27
3 2.5 − 4.0 2.2 −0.32 −0.091 0.084
4 4.0 − 5.0 1.8 −0.027 0.29 −0.56
5 5.0 − 5.8 1.5 −1.13 −1.39 0.58
Note that spline is symmetric, so the left of the y-axis is not listed since the polynomials
can be reflected to the negative side. The spline follows practically the contour of the
roof, with a small error near the nodes −0.8 and 0.8.