CHBE 230 - Lecture 3

CHBE 230 2019W2

Computational Methods
CHBE 230 - Lecture 3

Root finding

Common in engineering problems to solve f(x)=0 or

g(x)=a i.e. f(x)=0 with f(x)=g(x)-a
Different computational methods to approximate roots
A. Graphical method
B. Incremental search
C. Bisection and False position
D. Fixed-point iteration
E. Newton-Raphson Open methods
F. Secant & Modified secant

Optimization
CHBE 230 - Lecture 3

Fixed-point iteration
Iterative scheme based on successive substitutions
Basic principle
f (x) = 0 rewritten as x = g(x) ) xi+1 = g(xi )

Back to our Lecture 1 example: Find a root of the

following 2nd order algebraic equation
2
x 10x + 16 = 0 , solutions are x1 = 2 , x2 = 8
x2i + 16
Simple root finding algorithm: xi+1 = , x0 = G
10
x2i + 16
xi+1Initial
= guess:, x0 = G
10
Stopping criterion: |xi+1 xi | < tol
CHBE 230 - Lecture 3

Illustrate with Python

Let’s try with different initial guesses G:

G = 1 gives root = 2
G = 7 gives root = 2
G = 9 fails to converge
Why does our scheme either give 2 or fail to converge ?

Let’s try an alternative scheme

p
xi+1 = 10xi 16 , valid only if xi > 1.6
Let’s try with different initial guesses G with the alternative
scheme:
G > 2 gives root = 8
G < 2 fails to converge
CHBE 230 - Lecture 3

Graphical explanation: root is located at the intersection

y=x and y=g(x) divergence
Original scheme
y=x
Convergence/
divergence
Alternative
depends on
scheme
the slope of g:
|g’|<1

convergence to 2
x2 x1 G G x1 x2
CHBE 230 - Lecture 3

Why does convergence require |g’|<1?

See bottom of p. 155 in textbook.

Error changes in each iteration: Ei+1 = g′(xi)Ei

Illustrate by graph ……

Ei y=x
y=g(x)

Ei+1
Ei

Root!
xi xi+1
CHBE 230 - Lecture 3

Newton-Raphson method
Better scheme than fixed-point iteration, although
convergence is not guaranteed either
Fundamental to many numerical schemes as
linearization is almost the only tool we have to solve
non-linear problems (a key idea of CHBE 230)
Principle: locally approximate f(x) by a straight line
using its tangent i.e. f’(x)
Equation of tangent line:
0
y f (xi ) = f (xi )(x xi )
Where does it intersect y=0
f (xi )
at xi+1 defined by: xi+1 = xi
f 0 (xi )
CHBE 230 - Lecture 3

Graphical explanation

f(x)

f(xi) slope = f’(xi)

f(xi+1)
x
xi+1 xi

root
CHBE 230 - Lecture 3

Iterations until: x i+1 x i

< tol
xi+1
Two ways to compute f’(xi):
1. if f is not “too complicated”, derive analytically f
and use the analytical expression of f’ => standard
Newton-Raphson

2. approximation f’ using only values of f

A. using xi and a small step size h:
0 f (xi + h) f (xi )
f (xi ) '
h
B. using xi and xi-1 (secant line => the Secant
Method later)
0 f (xi 1 ) f (xi )
f (xi ) '
xi 1 xi
CHBE 230 - Lecture 3

Illustrate with Python

2
Newton-Raphson for 2nd order equation: x 10x + 16 = 0

Let’s try with different initial guesses: fails for G=5!

Why?
CHBE 230 - Lecture 3

Properties of Newton-Raphson algorithm:

Finds a single root
Fast quadratic convergence (when converges)
Convergence not guaranteed, see textbook page 159
Requires (a) initial guess “not too far away” from
solution, or equivalently (b) function “not too curved”.

f(x)

x1 x2
x0 x

root
CHBE 230 - Lecture 3

But Newton-Raphson does work for the following situation,

where the curve is tangential to the x-axis at the root, for
which the closed-interval bracketing methods won’t work:
2 2
f (x) = (x 2) .(x 3) over [1, 4]

f(x)

root root
CHBE 230 - Lecture 3

Secant method: approximate f’ by slope of secant line:

0 f (xi 1 ) f (xi )
f (xi ) '
xi 1 xi

Different from false-position method? Which is better?

CHBE 230 - Lecture 3

In our root-finding functions posted on Canvas, the

Secant Method is implemented using a small h, not
really the secant line:
0 f (xi + h) f (xi )
f (xi ) '
h

Illustrate with Python

With the Bungee jumper problem of determining the

mass of the jumper with g, cd, v at t known
r ✓r ◆
mg gcd
f (m) = tanh t̃ ṽ = 0
cd m
CHBE 230 - Lecture 3

Bonus example: a floating sphere, prob 5.22 of textbook

Volume of part above water

h ⇡h2
V = (3r h)
3
solid Inputs: fluid density,
solid density, sphere
liquid radius r

Questions:
1. formulate the equation for h
2. find h with the different root finding methods
CHBE 230 - Lecture 3

Solution:
1. We write “sum of the forces is equal to 0”
4 3
Fweight = ⇡r ⇢solid g
3 ✓ ◆
2
4 3 ⇡h
Fbuoyancy = ⇡r (3r h) ⇢f luid g
3 3
After some algebra, we get a 3rd order polynomial
equation in h:
3 2 3 ⇢f luid ⇢solid
h 3rh + 4r =0
⇢f luid
2. We use bisection or false position method to find the
root as this equation has a single root.
Please note that there is a solution only if ⇢f luid > ⇢solid
CHBE 230 - Lecture 3

Optimization
Common in engineering problems to find maxima or
minima of f(x). For example:

Minimum of energy dissipation in a mechanical

system (cars, machines, etc)
Maximum of chemical component production in a
chemical reaction

Minimum of f(x) and maximum of -f(x) are equivalent

problems
CHBE 230 - Lecture 3

Relation between min/max f(x) and roots of f’(x)

At any extremum of f(x), we have f’(x)=0
Max if f’(x)=0 and f’’(x)<0; Min if f’(x)=0 and f’’(x)>0
Similar operation, but optimization typically easier than
root finding for f’(x)=0.
Maximum:
f’(x)=0, f’’(x)<0

0
f(x)

Minimum:
f’(x)=0, f’’(x)>0
CHBE 230 - Lecture 3

Minimum bracketing algorithm

First, find bracket [xl, xu] that contains a single local
minimum (otherwise the scheme may fail)
Recall bisection method and adapt: take mid-point and
check f(xr) w.r.t. f(xu) and f(xl). But this is not enough!
f(x)
xl xr=0.5(xl+xu) xu
x
f(xu)
f(xl)
f(xr)<f(xl)
f(xr)
f(xr)<f(xu)
does not tell whether min is in [xl,xr] or [xr,xu]
CHBE 230 - Lecture 3

We need 1 more comparison, i.e. we need 2 interior points

f(x) d
d
xl x2 x1 xu
f(xu) x

f(x1)
f(xl)
f(x2)

Let’s define x1 and x2 as

x1 = xl + d
with d = c(xu xl ) , 0.5 < c < 1
x2 = xu d
CHBE 230 - Lecture 3

Now we compare f(x1) to f(x2)

If f(x2)< f(x1), then minimum is in [xl,x1]. Why?

f(x)

xl x2 x1 xu
f(xu) x

f(x1)
f(xl)
f(x2)

xl x1 xu,new=x1
f(x2)< f(x1) means there is at least one point (x2) in [xl,x1] at which f
falls between the f values at the two bounds.
CHBE 230 - Lecture 3

If f(x1)< f(x2), then minimum is in [x2,xu]

y
xl x2 x1 xu
x
f(xl)
f(x2)
f(x1)
f(xu)
xl,new=x2 x2 xu

Shrink the interval by 100.(1-c)% at each iteration, slightly

slower than bisection (recall 50% for bisection)
Did we have to choose c>0.5? What if c<0.5?
Then x1<x2, and we still eliminate the smaller portion of the
interval. So just a matter of notation; just need 2 inner points
CHBE 230 - Lecture 3

How to set c ?
Labor saving trick: set c such that
old x2 becomes new x1, or
old x1 becomes new x2
c 1-c

xl x2 x1 xu
1-c c
Eliminate
p
c1 1 −c c 5 1
=
= ) c= ' 0.618
1c 1c c 2
CHBE 230 - Lecture 3

Historical note: “divide a segment into 2 segments

such that the ratio of the whole to the greatest equals
the ratio of the greater to the smaller” (Euclid, 300BCE)
c 1-c p
1 c 5 1
= ) c= ' 0.618
c 1 c 2
1

The golden ratio is:

1 5+1
ϕ= =1+c= = 1.61803398874989......
c 2
CHBE 230 - Lecture 3

Illustration for the case: f(x2)< f(x1), minimum in [xl,x1]

c 1-c

xl x2 x1 xu
1-c c
Eliminate
We have x2,old = xl + (1 c)(xu,old xl )
Now x1,new in the interval [xl,x1,old]=[xl,xu,new]
x1,new = xl + c(xu,new xl ) = xl + c.c(xu,old xl )
2
But the golden ratio is constructed such that c = 1 c
x1,new = xl + (1 c)(xu,old xl ) = x2,old
CHBE 230 - Lecture 3

Pseudo code
Inputs: endpoints a,b; tolerance es; max number of
iterations maxit
Outputs: xmin, fmin=f(xmin), estimate error ea, number of
iterations for convergence i p
1+ 5
Initialization: i=0, xl,0=a, xu,0=b, =
2
Iterative loop: i=i+1
1. d=(φ-1)(xu,i-1-xl,i-1), x1,i=xl,i-1+d, x2,i=xu,i-1-d
2. if f(x1,i)<f(x2,i), xl,i=x2,i and xu,i=xu,i-1, xopt=x1,i
else xu,i=x1,i and xl,i=xl,i-1, xopt=x2,i
3. compute ea=100(2-φ)|(xu,i-xl,i)/xopt|
4. If ea<tol or i>=maxit, xmin=xopt, STOP
CHBE 230 - Lecture 3

Golden ratio search saves evaluation of f as x2,new=x1,old

or x1,new=x2,old
The pseudo code (or the code) does not take advantage
of this (How to improve the code? Save functional values)
Only one minimum that can be a local minimum
Error estimate in the case f(x2)< f(x1), xmin=x2 in [xl,x1]
f(x)
Local min can be in either interval

xl x2 x1 xu
f(xu) x

f(x1)
f(xl)
f(x2) max( l , 1)
ea = 100
Δl Δ1 xopt
CHBE 230 - Lecture 3

Error estimation: which of the two Δ’s is bigger?

l = x2 xl = xu ( 1)(xu xl ) xl
= (2 )(xu xl )
1 = x1 x2 = xl + ( 1)(xu xl ) (xu ( 1)(xu xl ))
= (2 3)(xu xl )

Since 2 ' 0.382 and 2 3 ' 0.2361

max( l, 1) = l = (2 )(xu xl )
xu xl
and hence the formula ea = 100(2 )
xopt

Remark: same applies to f(x1)< f(x2)

CHBE 230 - Lecture 3

Illustrate with Python

Example from our textbook

x2
Compute minimum of f(x) = − 2 sin x in [0,4].
10
Python code inside “Root_Finding_functions.py” and
the example inside “Root_Finding_Methods.ipynb”,
both posted on Canvas.