Chapter 1

Solving Equations
• Stewart platform
– As a six-degree-of-freedom, robot, the Stewart
platform can be placed at any point and
inclination in three dimensional space that is
within its reach
– A two-dimensional is discussed in the Reality
Check
• Given (x, y, θ), determine (p1, p2, p3)
• No closed-form solution is known

By Ethan Arnold, CC BY-SA 3.0,

https://commons.wikimedia.org/w/index.php?curid=8864535
 1.1 The Bisection Method
• Bracketing a root
– f(x) be a continuous function for x [a,b]
– f(a)f(b)<0 => exists a number r [a,b] such that
f(r)=0
– Bisection method

(tolerance)
• Example 1.1
Use Bisection Method to find the root of the function
f ( x ) = x + x − 1 on the interval [0,1]
3

>> f=@(x) x^3+x-1; >> f(a)*f(x0)

>> a=0;b=1;
>> f(a)*f(b) ans =
More iterations to improve
ans = 0.3750 the accuracy
-1 >> x1=(x0+b)/2
=> Programming

>> x0=(a+b)/2 x1 =

x0 = 0.7500

0.5000
• How accurate and how fast?
– [a0,b0] is the staring interval
– After n bisection steps, bn-an=(b0-a0)/2n , and the
solution is xc=(bn+an)/2
– Error< (bn-an)/2= (b0-a0)/2n+1
– Total function evaluations=n+2
– Definition1.3
• A solution is correct within p decimal places if the error
is less than 0.5x10-p
π = 3.141592653589793

πˆ1 ≈ 3.14  π − πˆ1 ≈ 0.00159

0.159 × 10−2 < 0.5 × 10−2  2 decimal places

πˆ 2 ≈ 3.141  π − πˆ 2 ≈ 0.00059
0.59 × 10−3 > 0.5 × 10−3
0.059 ×10−2 < 0.5 × 10−2  2 decimal places

πˆ3 ≈ 3.142  π − πˆ3 ≈ 0.00039

0.39 × 10−3 < 0.5 × 10−3  3 decimal places
• Example 1.2
Use Bisection Method to find a root of
f ( x ) = cos x − x in the interval [0,1] to within six correct
places

b − a 1− 0
En < n+1 = n+1 < 0.5 × 10−6 >> log2(1/0.5e-6)
2 2
1 ans =
2n+1 >
0.5 × 10−6 20.931568569324174
 1 
log 2 2n+1 > log 2  −6 
 0.5 × 10 
n > 19.93
 – See & run “ex1p2.m”

Computed Solutions(xc)
1 Error bound & errors
100

0.8
10-2

0.6
10-4
0.4

10-6
0.2

0 10-8
0 5 10 15 20 0 5 10 15 20
• Bisection summary
– A stable method that always converge
– We can know how many steps is required for a
desired precision
– By contrast, more high-power algorithms are
often less predictable
• Not necessarily converge
• Need to define stopping criteria

 1.2 Fixed-Point Iteration
• Fixed points of a function
– Def. 1.4: The real number r is a fixed point of the
function g if g(r)=r
• g(x)=cos(x)=> r=0.7390851332
• g(x)=x3=> r=-1, 0, 1
– Once equation is written as g(x)=x, fixed-point
iteration proceeds as inital guess  x0
xi +1 = g ( xi )
– Not necessarily converge
 – Recalculate EXAMPLE 1.1: x3+x-1=0
• x0=0.5, r=0.682327803828020

Three possible arrangements:

(a) >> g=@(x) [1-x(1)^3, (1-x(2))^(1/3), (1+2*x(3)^3)/(1+3*x(3)^2)];
x = 1 − x3 >> x=[0.5 0.5 0.5];
>> x=g(x)
(b)
x=
x = 3 1− x
(c) 0.8750 0.7937 0.7143

3 x3 + x = 1+2 x 3 >> x=g(x)

x ( 3 x 2 + 1) = 1+2 x 3
x=
1+2 x 3
x= 0.3301 0.5909 0.6832
1 + 3x 2
 – Geometry of fixed-point iteration

1+2 x 3
(a) g ( x ) = 1 − x 3 (b) g ( x ) = 3 1 − x (c) g ( x ) =
1 + 3x 2
– Linear convergence of fixed-point iteration
• Geometric view
– |g’(r)| makes the crucial difference between divergence and
convergence

3 5 1 3
g1 ( x ) = − x + g2 ( x ) = − x +
2 2 2 2
r =1 r =1
3 1
g1 ( r ) = > 1
′ g1 ( r ) = < 1
′
2 2
 • In terms of equations, it helps to write g1(x) and g2(x) in
terms of x-r
– The error is ei=|xi-r| for i-th iteration

3 5 1 3
g1 ( x ) = − x + g2 ( x ) = − x +
2 2 2 2
3 1
= − ( x − 1) + 1 = − ( x − 1) + 1
2 2
3 1
g1 ( x ) − 1 = − ( x − 1) g 2 ( x ) − 1 = − ( x − 1)
2 2
3 1
xi +1 − 1 = − ( xi − 1) xi +1 − 1 = − ( xi − 1)
2 2
3 1
ei +1 = ei ei +1 = ei
2 2
– Def. 1.5 :
ei +1
lim = S < 1  Linear convergence with rate S
i →∞ e
i

– Def. 1.7: A iterative method is called locally

convergence to r if the method converges to r for
initial guesses sufficiently close to r

– Theorem 1.6: If S=|g’(r)|<1, fixed-point iteration

converges linearly with rate S to the fixed point r
for initial guesses sufficiently close to r
 • Proof of Theorem 1.6

Mean Value Theorem xi +1 = g ( xi )

f (b ) − f ( a ) xi +1 ≈ xi ≈ r
f ′(c) =
b−a
g ( xi ) − g ( r ) xi +1 − r
c ∈ [ a, b ] g′( c ) = =
xi − r xi − r
c ∈ [ xi , r ]  c ≈ r
xi +1 − r = g ′ ( c ) [ xi − r ]
ei +1 = g ′ ( c ) ei
ei +1 ≈ g ′ ( r ) ei
 – Revisit EXAMPLE 1.1: x3+x-1=0
• x0=0.5, r=0.682327803828020

(c)
(a) (b)
3 3 x 3 + x = 1+2 x 3
x = 1− x x = 3 1− x
g ′ ( r ) = 1.3966 > 1 x ( 3 x 2 + 1) = 1+2 x 3
g ′ ( r ) = 0.716 < 1
1+2 x 3
x=
1 + 3x 2
g′( r ) = 0
• Example 1.3
Explain why the FPI g(x)=cos(x) converges

• Example 1.4
Use FPI to find a root of cos(x)=sin(x)
• Example 1.5
Find the fixed point of g(x)=2.8x-x2

g ′ ( 0 ) = 2.8
g ′ (1.8 ) = −0.8
• Example 1.6
Calculate by using FPI
x= 2
2
x =2
(b)
(a)
2
2 =x
=x x
x 2
2 + x = 2x
xi +1 = x
xi x+2 x
=x
=> won’t converge 2
=> rapidly converge
• YBC 7289 is a Babylonian clay tablet
– Some time in the range from 1800–1600 BC
» The Babylonians’ method is not known, but some
speculate it is based on the FPI algorithm in base 60
» In the first century A.D., Heron of Alexandria used the FPI
algorithm to calculate 720
 – Stopping criteria

1.

2.

3.

 1.3 Limits of Accuracy
• Forward and backward error
– Forward error: |r-xa|
– Backward error: |f(xa)|
– Well-conditioned problems
• Both errors approach zero, e.g.,
>> syms x
>> f1=(x-1)*(x-2)*(x-3)*(x-4)*(x-5);
>> expand(f1)
ans = x^5 - 15*x^4 + 85*x^3 - 225*x^2 + 274*x - 120
>> f1=@(x) x^5 - 15*x^4 + 85*x^3 - 225*x^2 + 274*x - 120;
>> fzero(f1,5)
ans = 5
>> f1(5)
ans =0
 – Ill-conditioned problems
• Forward error can hardly approach zero
• Example 1.7
3
4 8  2
f ( x ) = x3 − 2 x 2 + x− =x− 
3 27  3
>> f=@(x) x^3-2*x^2+4/3*x-8/27;
>> fzero(f,1)

ans =

0.666664247279230

>> f(0.666664247279230)

ans =

0
 • Wilkinson polynomial

>> syms x
>> f2=(x-1)*(x-2)*(x-3)*(x-4)*(x-5)*(x-6)*(x-7)*(x-8)*(x-9)*(x-10)*...
(x-11)*(x-12)*(x-13)*(x-14)*(x-15)*(x-16)*(x-17)*(x-18)*(x-19)*(x-20);
>> expand(f2)

ans =

x^20 - 210*x^19 + 20615*x^18 - 1256850*x^17 + 53327946*x^16 - 1672280820*x^15 +

40171771630*x^14 - 756111184500*x^13 + 11310276995381*x^12 - 135585182899530*x^11
+ 1307535010540395*x^10 - 10142299865511450*x^9 + 63030812099294896*x^8 -
311333643161390640*x^7 + 1206647803780373360*x^6 - 3599979517947607200*x^5 +
8037811822645051776*x^4 - 12870931245150988800*x^3 + 13803759753640704000*x^2 -
8752948036761600000*x + 2432902008176640000
>> f2=@(x) x^20 - 210*x^19 + 20615*x^18 - 1256850*x^17 + 53327946*x^16 -
1672280820*x^15 + 40171771630*x^14 - 756111184500*x^13 + 11310276995381*x^12 -
135585182899530*x^11 + 1307535010540395*x^10 - 10142299865511450*x^9 +
63030812099294896*x^8 - 311333643161390640*x^7 + 1206647803780373360*x^6 -
3599979517947607200*x^5 + 8037811822645051776*x^4 - 12870931245150988800*x^3 +
13803759753640704000*x^2 - 8752948036761600000*x + 2432902008176640000;

>> fzero(f2,16)

ans =

16.014680305804578

>> f2(16.014680305804578)

ans =

-1.089172930560000e+11

 1.4 Newton’s Method
• Given a function f(x), finding f(xr)=0
– Initial guess: x0
– Express f(x) with Taylor series @ x0
f ( x ) = f ( x0 ) + f ′ ( x0 )( x − x0 ) +

– Find the root of the truncated linear equation

f ( x0 ) + f ′ ( x0 )( x1 − x0 ) = 0
f ( x0 )
x1 = x0 −
f ′ ( x0 )
– Iteration formula
f ( xi )
xi +1 = xi − , i = 0,1,2
f ′ ( xi )

-1 x3 = 142.35
f ( x)
f(m)

-2

-3
x0 = 68
-4

-5
50 100 150 200
x
m [kg]
• Example 1.11
Find the Newton’s Method formula for the
3
equation x + x − 1 = 0

f ( x ) = x3 + x − 1
f ′ ( x ) = 3x 2 + 1
f ( xi ) xi3 + xi − 1
xi +1 = xi − = xi −
f ′ ( xi ) 3 xi2 + 1
2 xi3 + 1
= 2 >> roots([1 0 1 -1])
3 xi + 1
ans =

-0.341163901914009 + 1.161541399997252i
-0.341163901914009 - 1.161541399997252i
0.682327803828020 + 0.000000000000000i
 – See & run ex1p11.m
0: xc=-7.000000e-01 error=1.382328e+00
1: xc=1.271255e-01 error=5.552023e-01
2: xc=9.576781e-01 error=2.753503e-01
3: xc=7.348278e-01 error=5.249999e-02 =>Quadratic convergence
4: xc=6.845918e-01 error=2.263967e-03
5: xc=6.823322e-01 error=4.370376e-06
6: xc=6.823278e-01 error=1.631228e-11
Computed Solutions(x )
c
10

2 x2
0

-2 x1
x0
-4
-1 -0.5 0 0.5 1 1.5 2
 – Revisit Example 1.2
• Solve f ( x ) = cos ( x ) − x = 0

Computed Solutions(xc) Errors

1 100
Bisection Bisection
Newton Newton

0.8
10-5

0.6
10-10
0.4

10-15
0.2

0 10-20
0 5 10 15 20 0 5 10 15 20
• Quadratic convergence of Newton’s method
– Def. 1.10
An iteration method is quadratic convergence if
ei +1
lim 2 = M < ∞
i →∞ e
i

– Theorem 1.11
Newton’s Method is locally and quadratically
convergent to r &

f ′′ ( r )
M=
2 f ′( r )
 – Proof of Theorem 1.11
• To prove convergence
– Newton’s method is a particular form of Fixed-Point Iteration

FPI: xi +1 = g ( xi )
f ( xi )
Newton's: xi +1 = xi −
f ′ ( xi )
f ( x)
 g ( x) = x −
f ′( x)
2
 f ′ ( x )  − f ( x ) f ′′ ( x ) f ( x ) f ′′ ( x )
g′( x) = 1 − 2
= 2
 f ′ ( x )   f ′ ( x ) 
f ( r ) f ′′ ( r )
g′( r ) = 2
= 0 <1
 f ′ ( r ) 

=> Newton’s Method is locally convergent according to Theorem 1.6

• To prove quadratic convergence
– Expand f(x) at xi
f ′′ ( ci )
f ( x ) = f ( xi ) + f ′ ( xi )( x − xi ) + ( x − xi ) ci ∈ [ x, xi ]
2
,
2
( r − xi )
2

f ( r ) = f ( xi ) + f ′ ( xi )( r − xi ) + f ′′ ( ci )
2
( r − xi )
2

0 = f ( xi ) + f ′ ( xi )( r − xi ) + f ′′ ( ci )
2
− f ( xi ) ( r − xi ) f ′′ ( ci )
2

= ( r − xi ) +
f ′ ( xi ) 2 f ′ ( xi )
f ( xi ) ( r − xi ) f ′′ ( ci )
2

xi − −r =
f ′ ( xi ) 2 f ′ ( xi )
f ′′ ( ci )
( xi − r )
2
xi +1 − r =
2 f ′ ( xi )
f ′′ ( ci ) 2 f ′′ ( r ) 2
ei +1 = ei ≈ ei
2 f ′ ( xi ) 2 f ′(r )
– Algorithm for finding

f ( x) = x − a = 0
2

f ′( x ) = 2x
a
xi +
f ( xi ) 2
xi − a xi
xi +1 = xi − = xi − =
f ′ ( xi ) 2 xi 2
For a = 2  (1.14 )
• Linear convergence of Newton’s method
– Newton's Method does not always converge
quadratically
– FPI does not always converge linearly
• Example 1.12
Use Newton’s Method to find a root of

• Example 1.13
Use Newton’s Method to find a root of
– Theorem 1.13
If f contains a root r of multiplicity m>1, then Modified
Newton’s Method
mf ( x )
xi +1 = xi −
f ′( x )
converges locally and quadratically to r

f ( x ) = x3
• Example 1.15
Apply Newton’s Method to

with starting guess x0=1/2

=> Like FPI, Newton’s method may not

converge to a root
• Summary of Newton’s Method

f ( x ) = x3 − 3sin ( x )

>> fx=@(x) x.^3-3*sin(x); % function >> x=x-fx(x)./fxp(x) % second iteration

fxp=@(x) 3*x.^2-3*cos(x);% derivative
x=
>> x=[-1 -0.82 0.5 1]; % initial guess -1.6500 36.6081 0.0082 1.6500
>> x=x-fx(x)./fxp(x) % first iteration
>> x=x-fx(x)./fxp(x) % third iteration
x=
-2.1054 54.9121 -0.1975 2.1054 x=
-1.4714 24.4006 -0.0000 1.4714

 1.5 Root-Finding Without Derivatives
• Secant method and variants
– Secant method: Use first order polynomial g(x) to
interpolate f(x)
1 ⊕
x1
0 f ( x) f ( x1 ) − f ( x0 )
x2 g ( x ) = f ( x0 ) + ( x − x0 )
-1 x1 − x0
f (m)

-2 f ( x1 ) − f ( x0 )
g ( x2 ) = f ( x0 ) + ( x2 − x0 )
g ( x) x1 − x0
-3
x1 − x0
x2 = x1 − f ( x0 )
-4 f ( x1 ) − f ( x0 )
⊕ x0
-5
50 100 m 150 200
 – Iteration formula & convergence rate

Initial guesses: x0 & x1

xi − xi −1
xi +1 = xi − f ( xi ) , i = 1,2,
f ( xi ) − f ( xi −1 )

α −1
f ′′ ( r )
ei +1 ≈ eiα
2 f ′( r )
1+ 5
where α = ≈ 1.62  Superlinear convergence
2
• Example 1.16
Apply the Secant Method with starting guesses
x0=0, x1=1 to find the root of

>> fx=@(x) x^3+x-1;

x0=0;x1=1;
x2=x1-(x1-x0)*fx(x1)/(fx(x1)-fx(x0))

x2 =

0.5000

>> x3=x2-(x2-x1)*fx(x2)/(fx(x2)-fx(x1))

x3 =

0.6364
– Muller’s method and IQI
• Use three points to draw the parabola

Muller: IQI:
y = P2 ( x ) x = P2 ( y )

y = f ( x)
• Brent’s method
– Bisection + Secant + IQI
>> f=@(x) x^3+x-1;
fzero(f,[0 1],optimset('Display','iter'))

Func-count x f(x) Procedure

2 1 1 initial
3 0.5 -0.375 bisection
4 0.636364 -0.105935 interpolation
5 0.68491 0.00620153 interpolation
6 0.682225 -0.000246683 interpolation
7 0.682328 -5.43508e-07 interpolation
8 0.682328 1.50102e-13 interpolation
9 0.682328 0 interpolation

Zero found in the interval [0, 1]

ans =
0.682327803828019
>> f=@(x) x^3+x-1;
fzero(f,1,optimset('Display','iter'))

Search for an interval around 1 containing a sign change:

Func-count a f(a) b f(b) Procedure
1 1 1 1 1 initial interval
3 0.971716 0.88924 1.02828 1.11556 search
5 0.96 0.844736 1.04 1.16486 search
7 0.943431 0.783145 1.05657 1.23606 search
9 0.92 0.698688 1.08 1.33971 search
11 0.886863 0.584404 1.11314 1.4924 search
13 0.84 0.432704 1.16 1.7209 search
15 0.773726 0.236918 1.22627 2.07028 search
16 0.68 -0.005568 1.22627 2.07028 search
Search for a zero in the interval [0.68, 1.22627]:
Func-count x f(x) Procedure
16 0.68 -0.005568 initial
17 0.681465 -0.00206575 interpolation
18 0.682329 2.05001e-06 interpolation
19 0.682328 -1.51068e-09 interpolation
20 0.682328 -1.22125e-15 interpolation
21 0.682328 0 interpolation

Zero found in the interval [0.68, 1.22627]

ans =

0.6823
• Revisit Example 1.2
– Solve
• Use IQI to solve Example 1.2
– Solve
X 0 = [ 0 0.5 1]
X 1 = [ 0.5 0.75 1]
x0 = 0.5 → x1 = 0.75 → x2 = 0.739 →

6
x = −0.0777 y 2 − 0.6036 y + 0.7390
4

0 f ( x)
-2

-4
x = −0.1412 y 2 − 0.6088 y + 0.75
-6
-6 -5 -4 -3 -2 -1 0 1 2
Computed Solutions(xc)

Bisection
Newton
FPI
Secant
IQI
Errors

Bisection
Newton
FPI
Secant
IQI
• A circuit example
( )
I D = I S eVD VT − 1 ≈ I S eVD VT , when VD >> VT
I S = 10−15 A
VT ≈ 27 mV
100 Ω

ID
+
4V
VD
−
