Master Course in Numerical Mathematics 65

Since zx = 2x − 2 − y and zy = 2y + 1 − x, a point (x0; y0), such that

5 Optimization
zx (x0; y0) = zy (x0; y0) = 0
In optimization problems
has to solve we are
the linear seeking the
equation minimum or a relative minimum of a real-valued function
system
f . (For searching a maximum " of f , we will#simply
" search
# "a minimum # of −f .)
2 −1 x 2
If f depends on one variable only, you know ·that f0 has=a relative minimum
. at x0 , if (4.1)
−1 2 y0 −1
f 0 (x0 ) = 0 and f 00 (x0 ) > 0.
The unique solution (x0; y0) = (1, 0) can be determined by Cramer’s rule or by
Gaussian
If f depends on elimination.
two variables, Ityou
turns
know outthatthatf has a relative minimum at (x0 , y0 ), if
" #! " #!
∂f F (x ; y ) F
xx 0 0∂f xy 0 0 (x ; y ) 2 −1 ∂ 2f
det
(x0 , y0 ) = 0 and (x0 , y0 ) = 0 and = det
det(Hess(f )) > 0 and =(x30 ,>y00,
) > 0,
∂x Fyx(x0; y0∂y
) Fyy (x0; y0) −1 2 ∂x2
wherewhich
Hess(fgives
) denotes the Hessian matrix
that (x of the second order derivatives of f .
evidence 0; y0 ) = (1, 0) is either a local maximum or mini-
mum.
However, Fxx (x0; y0)the
Since application,
in many = 2derivatives
> 0, it is may
a minimum.
be too expensive to compute. The next section
will therefore introduce an optimization method, which does not use derivatives.

Remark Note that the minimum (x0; y0) = (1, 0) in example 4.1 can be found
5.1 The Nelder-Mead Simplex Method
analytically, since the equation system (4.1) is linear. In general, when z =
F (x; y)is iscourtesy
This section a function withDr.
of Prof. exponential, logarithmic,
A. Rodenhausen, or sine and cosine terms, the
HAW Hamburg.
equation system (4.1) is non-linear and can only be solved by use of numerical
methods. Even powers xn and y n lead to non-linear systems in (4.1) for n > 2.

The basic idea of the Nelder-Mead algorithm

The principal idea of the Nelder-Mead simplex method is to construct a se-
quence of simplexes
T0, T1, . . . , Tm ⊂ Rn .

y v3 z
v3
x
y
x v4
Figure 4.1: A simplex in R2
is a triangle (left). In R3 it is
a tetrahedron (right).
v1
v2 v1
v2

14
66 HAW Hamburg — Faculty of Life Sciences — Prof. Dr. Anna Rodenhausen

A simplex in Rn is given by n + 1 vertices v1, v2, . . . , vn+1. The simplex is the

convex hull of the points v1, v2, . . . , vn+1. i. e., it is the smallest convex set in
Rn that contains v1 , v2, . . . , vn+1. A simplex is a triangle for n = 2, and it is a
tetrahedron for n = 3 (see figure 4.1).
For any Nelder-Mead simplex Tk we find Tk+1 by replacing the vertex W by
another vertex with a smaller function value. The new vertex is found by one
of the four steps described in the following for a function z = F (x1; x2) with
n = 2 variables. The steps are in principal the same for n > 2.

Four possible steps to replace the worst vertex

Let B, G, and W (standing for best, good, and worst) denote the vertices of the
initial triangle in R2 . By definition, the function values at B, G, and W satisfy

F (B) ≤ F (G) ≤ F (W ).

Step R: Reflection of W at M B
Let M denote the midpoint of the edge connecting
R
B and G, and let R be the reflection point of W T2
at M. Than M
B+G W
T1
M= and R = 2M − W.
2
If F (R) < F (W ) is satisfied, replace W by R. G

Step E: Extention of R to E
B
E When F (R) < F (W ) was satis-
R fied in step R, it might be the case
T2
that F (E) < F (R), where
M
W
T1 E = 2R − M.

G
If this is the case, replace W not by
R but even by E.

15
Master Course in Numerical Mathematics 67

Step C: Contraction of W
When the steps R and E did not succeed, determine
W −M
C1 , C2 = M ± .
2
Try whether F (C1) < F (W ) or F (C2) < F (W ). When both C1 and C2 have
the desired property, take the point with the smaller function value as C. Re-
place W by C.
B B

C2 C2
T2 T2
C1 M C1 M
W W
T1
T1

G G

B Step S: Shrinking of W to S
When steps R, E, and C have not been successful,
S T2 shrink the whole triangle T1 in the direction of the best
M vertex B. Explicitly, replace G by M and W by
W T1
B+W
S= .
2
G

The size of the simplexes T0, T1, . . . , Tm, . . . is finally decreasing for increasing
m. When the maximum distance between two vertices of Tm is less than a
prescribed precision, the algorithm terminates and returns the best vertex of Tm
as a local minimum of z = F (x1; x2; . . . ; xn).

16
68 HAW Hamburg — Faculty of Life Sciences — Prof. Dr. Anna Rodenhausen

The logical structure behind the steps R, E, C, and S

if ( F(R) < F(W) ) if ( F(E) < F(R) )

replace W by E. step E

else
replace W by R. step R

if ( F(C1) < F(W) ) if ( F(C2) < F(C1) )

replace W by C2
step C
else
replace W by C1
else replace G by M and
replace W by S step S

Example 4.2: — Running the Nelder-Mead simplex algorithm —

Consider again the function f (x; y) = x2 + y 2 − 2x + y − xy from example
0.1 with a local minimum point at (x0; y0 ) = (1; 0). We initialize the simplex
algorithm with the triangle T0 spanned by the vertices (−1; −2), (−3; −3), and
(−4; 0). The vertices of the resulting sequence of triangles T0, T1, T2, . . . , T11
and their function values are given in the table. The second column of the table
shows a letter R, E or C indicating which replacement was recently made (the
shrinking step S does not occur in this example):
k step B f (B) G f (G) W f (W )
0 − (−1.00; −2.00) 3.00 (−3.00; −3.00) 12.00 (−4.00; 0.00) 24.00
1 R (−1.00; −2.00) 3.00 (−3.00; −3.00) 12.00 (0.00; −5.00) 20.00
2 C (−1.00; −2.00) 3.00 (−1.00; −3.75) 9.56 (−3.00; −3.00) 12.00
3 R (−1.00; −2.00) 3.00 (1.00; −2.75) 6.56 (−1.00; −3.75) 9.56
4 E (2.00; 0.38) −0.23 (−1.00; −2.00) 3.00 (1.00; −2.75) 6.56
5 R (2.00; 0.38) −0.23 (0.00; 1.12) 2.39 (−1.00; −2.00) 3.00
6 C (2.00; 0.38) −0.23 (0.00; −0.62) −0.23 (0.00; 1.12) 2.39
7 C (0.50; 0.50) −0.25 (2.00; 0.38) −0.23 (0.00; −0.62) −0.23
8 C (0.62; −0.09) −0.89 (0.50; 0.50) −0.25 (2.00; 0.38) −0.23
9 C (1.28; 0.29) −0.92 (0.62; −0.09) −0.89 (0.50; 0.50) −0.25
10 R (1.28; 0.29) −0.92 (0.62; −0.09) −0.89 (1.41; −0.30) −0.62
11 C (1.18; −0.30) −0.94 (1.28; 0.29) −0.92 (0.62; −0.09) −0.89
.. .. .. .. ..
. . . . .
25 C (1.00; 0.00) −1.00 (1.01; 0.01) −1.00 (1.01; −0.01) −1.00

17
Master Course in Numerical Mathematics 69

After 25 steps of optimization, two decimals after the decimal point of the best
vertex B = (1.00; 0.00) are exact. Figure 4.2 gives an illustration of the trian-
gles T0 , T1, T2, . . . , T11.

T11
T8
T6

-4 -3 -2 -1 T9 2 3
T7
T10

T5

T4
T0
-2

T3
T2 -3

T1 -4

-5

Figure 4.2: The triangles T0 , T1 , T2 , . . . , T11 , which are obtained with the initial triangle T0 spanned by
(−1; −2), (−3; −3) and (−4; 0). The minimum is represented by the little white disc inside of T11 . It
can be observed that the triangles capture the minimum point systematically.

Calling the MATLAB function fminsearch

The MATLAB function fminsearch is called by the command line

[min, fval] = fminsearch(’func two var’, x).

18
 Numerical Mathematics (25w) Lecture Notes Prof. Dr.-Ing. T. Schiemann 19

6 Numerical Differentiation

Differentiation is usually not a big problem, because the derivatives of a function defined by ana-
lytical expressions can always be determined. However, we may want to compute derivative values
numerically, if the function to be differentiated is not defined by a formula but by a table of values,
e.g. originating from a measurement.

On the other hand, the formulas for numerical differentiation, which we will see in this chapter, may
be used to derive other formulas, which require a general formula for the derivative of a function.
We will see this in the chapter about partial differential equations at the end of the semester.

The following two sections introduce the most basic formulas for numerical computation of first and
second order derivatives using the approximation by difference quotients.

6.1 First order derivative

The first order derivative of a function f : R → R is defined by

f (x + h) − f (x)
f 0 (x) = lim .
h→0 h
If h is small enough, we can set
f (x + h) − f (x)
f 0 (x) ≈ .
h

This is the approximation of f 0 (x) by the right difference, because on the x-axis, the position x + h
is right of the position x.

In the same way we can construct the approximation of f 0 (x) by the left difference:
f (x) − f (x − h)
f 0 (x) ≈ .
h

When we add the formulas for the left and the right difference, we obtain the approximation of f 0 (x)
by the central difference:
f (x + h) − f (x − h)
f 0 (x) ≈ .
2h

ln x
Example: f (x) = , x0 = 0.2, h = 0.01
2 + cos x
1 ln 0.2 ln 0.19


Left difference: f 0 (x0 ) ≈ 0.01
· 2+cos 0.2
− 2+cos 0.19
≈ 100 · (−0.540068 + 0.556918) ≈ 1.685
1 ln 0.21 ln 0.19


Central difference: f 0 (x0 ) ≈ 0.02
· 2+cos 0.21
− 2+cos 0.19
≈ 100 · (−0.524054 + 0.556918) ≈ 1.643

The exact value is f 0 (x0 ) ≈ 1.642.

20 Numerical Mathematics (25w) Lecture Notes Prof. Dr.-Ing. T. Schiemann

6.2 Error-order of a numerical method

The formulas for numerical differentiation deliver only an approximation of the true value of the
derivative. We are therefore observing an error E. In the case of numerical differentiation, the error
depends on the width h, which has been chosen in the formulas above.

It is important to study, how the error-function E(h) is developing, when we change the value of h.

For typical numerical methods, the error-function is close to a power-function:

E(h) ≈ C · hr ,

where C is some constant factor and r is a constant exponent, which is called the order of the error.

For r = 1, we speak of a linear error-order: E(h) = C · h. In this case, the error E is proportional
to the width h. Reducing h by a certain factor would reduce the error by the same factor.

For r = 2, we speak of a quadratic error-order: E(h) = C · h2 . Reducing h by a certain factor would

reduce the error by the square of this factor, e.g. setting h to 50% of its original value would reduce the
error to 25% of its original value. Methods with quadratic error-order are therefore usually preferred
compared to methods with linear error-order. (But better methods could be computationally more
expensive...)

It is important to note, that all error-considerations are made under the assumption that
the error is “small enough”, which means that we are close to the true value. In practical
situations, it is not always clear, if this is the case, because the true value is unknown.

The error-order r can be determined, if we know the error E for at least two different values of h:
log(E1 /E2 )
r=
log(h1 /h2 )
This formula can best be applied, if h1 and h2 are related by a factor of 10, e.g. h1 = 0.01 and
h2 = 0.001, because the logarithms of h1 /h2 is then an integer number. If the errors are small
enough, their quotient E1 /E2 is usually close to a power of 10, so that also the logarithm of E1 /E2
can easily be evaluated.

Example: (The values are taken from the table on the next page)
h1 = 0.01, h2 = 0.001, E1 = 3.4E-4, E2 = 3.4E-6 ⇒ h1 /h2 = 10, E1 /E2 = 100 ⇒ r = 2

The strict mathematical way of determining the error order would be a proof that the error-function
E(h) is a power-function with certain values for C and r. This proof would be independent of certain
h- or E-values and would instead use the definition of the numerical method only.

6.3 Error-order for approximation of first order derivatives

Instead of a precise error analysis we will only look at some examples as “mathematical experiments”.

The table below shows the left, right and central differences for the h-values h1 = 0.01 and h2 = 0.001
for some functions f (x) at the position x0 = 3. As we know the exact derivatives of these functions,
we can compute the errors E1 and E2 .

From this table we conclude that the approximation by the left and the right differences show nearly
identical errors, and both of them have a linear error order, because the reduction of the width h by
 Numerical Mathematics (25w) Lecture Notes Prof. Dr.-Ing. T. Schiemann 21

a factor of 10 does also reduce the error by a factor of 10. This is the typical linear/proportional
behaviour.

The approximation by the central difference has a quadratic error order, because the reduction of
the width h by a factor of 10 reduces the error by a factor of 100(= 102 ).

Absolute Errors at x0 = 3
Exact Left difference Right difference Central difference
Function f 0 (x) h1 = 0.01 h2 = 0.001 h1 = 0.01 h2 = 0.001 h1 = 0.01 h2 = 0.001
x2 6.00 1.00E-02 1.00E-03 1.00E-02 1.00E-03 1.28E-13 6.61E-13
x3 27.00 8.99E-02 9.00E-03 9.00E-02 9.00E-03 1.00E-04 1.00E-06
x4 108.00 5.39E-01 5.40E-02 5.41E-01 5.40E-02 1.20E-03 1.20E-05
1/x -0.11 3.72E-04 3.70E-05 3.69E-04 3.70E-05 1.23E-06 1.23E-08
ex 20.09 1.00E-01 1.00E-02 1.01E-01 1.00E-02 3.35E-04 3.35E-06
Absolute errors for approximation of the first order derivatives for some func-
tions at the position x0 = 3.

For the function f (x) = x2 , the central difference does even deliver exact values.

For polynomials, it can easily be shown that the central difference has a quadratic error order: let
E(h) be the error for the function f (x) = x3 :

(x + h)3 − (x − h)3
E(h) = 3x2
− |{z}
| 2h
{z } exact
central diff.
x3 + 3x2 h + 3xh2 + h3 − (x3 − 3x2 h + 3xh2 − h3 )
= − 3x2
2h

6x2 h + 2h3
= − 3x2
2h

= 3x2 + h2 − 3x2

= h2

Extending this example towards a Taylor-series for some general function f leads to a general proof
of the error orders for the difference-formulas.
22 Numerical Mathematics (25w) Lecture Notes Prof. Dr.-Ing. T. Schiemann

6.4 Second order derivative

For approximation of the second order derivative f 00 at a position x0 we first are first considering the
function
f (x0 − h) − 2f (x0 ) + f (x0 + h) f (x0 + h) − f (x0 − h)
p2 (x) := 2
(x − x0 )2 + (x − x0 ) + f (x0 )
2h 2h
for some value h ∈ R.

p2 has the property to interpolate f at the positions x0 , x0 − h, and x0 + h: p2 (x0 ) = f (x0 ),

p2 (x0 − h) = f (x0 − h), and p2 (x0 + h) = f (x0 + h).

We therefore use the derivatives of p2 as approximations of the derivatives of f and we obtain

the central difference for the first order derivative and the following formula for the second order
derivative:
f (x + h) − 2f (x) + f (x − h)
f 00 (x) ≈
h2

From the following table we conclude that this approximation of f 00 has a quadratic error order.

Exact Absolute Errors

Function f 00 (x) h=0.01 h=0.001
x4 108.00 2.00E-04 2.00E-06
1/x 0.007 8.23E-07 8.29E-09
ex 20.09 1.67E-04 1.67E-06
Absolute errors for approximation of the second order derivatives for some
functions at the position x = 3.