Nonlinear Equations Nonlinear Equations The scientist described - - PowerPoint PPT Presentation

Nonlinear Equations Nonlinear Equations

The scientist described what is: the engineer creates what never was. Theodor von Karman

1

The father of supersonic flight Fall 2010

Problem Description Problem Description

Given a non-linear equation f(x)=0, find a x* q f( ) such that f(x*) = 0. Thus, x* is a root of f(x)=0. Galois theory in math tells us that only l i l f d ≤ 4 b l d ith polynomials of degree ≤ 4 can be solved with close forms using +, −, ×, ÷ and taking roots. General non-linear equations can be solved with General non-linear equations can be solved with iterative methods. Basically, we try to guess the location of a root, Basically, we try to guess the location of a root, and approximate it iteratively. Unfortunately, this process can go wrong, leading to another root or even diverge.

2

Methods Will Be Discussed Methods Will Be Discussed

There are two types of methods, bracketing and There are two types of methods, bracketing and

• pen. The bracketing methods require an

interval that is known to contain a root, while interval that is known to contain a root, while the open method does not. Commonly seen bracketing methods include Commonly seen bracketing methods include the bisection method and the regula falsi method and the open methods are Newton’s method, and the open methods are Newton s method, the secant method, and the fixed-point iteration iteration.

3

Bisection Method: 1/6 Bisection Method: 1/6

Idea: The key is the intermediate value theorem. Idea: The key is the intermediate value theorem. If y = f(x) is a continuous function on [a,b] and r is between f(a) and f(b) then there is a x* such is between f(a) and f(b), then there is a x such that r = f(x*). Th f if f( )×f(b) < 0 d 0 th i Therefore, if f(a)×f(b) < 0 and r = 0, there is a root x* of f(x) = 0 in [a, b].

f(a) y=f(x) f(b) a b x* X

4

f(b)

Bisection Method: 2/6 Bisection Method: 2/6

Note that f(a)×f(b) < 0 guarantees a root in [a,b]. f( ) f( ) g Compute c = (a+b)/2 and f(c). If f(c) = 0, we have found a root. f Otherwise, either f(a)×f(c) < 0 or f(b)×f(c) < 0 but not both. Use [a,c] for the former and [c,b] for th l tt til |f( )| < the latter until |f(c)| < ε .

1 2 f(a) y=f(x) 3 4 f(a) a b x* X 4

5

f(b) x*

Bisection Method: 3/6 Bisection Method: 3/6

Conver Convergence ence: Since the first iteration reduces g the interval length to |b-a|/2, the second to |b-a|/22, the third to |b-a|/23, etc, after the k-th iteration, the interval length is |b-a|/2k. If ε > 0 is a given tolerance value, we have |b-a|/2k

k

< ε or |b-a|/ε < 2k for some k. Taking base 2 logarithm, we have the value of k, the expected number of iterations with accuracy ε:

k b a − ⎡ ⎤ l

Convergence is not fast but guaranteed and steady.

k = ⎡ ⎢ ⎢ ⎤ ⎥ ⎥ log2 ε

6

g g y

Bisection Method: 4/6 Bisection Method: 4/6

Al Algorithm:

• rithm:

Fa = f(a)

g If ABS(b-a) < ε, then stop.

Fa = f(a) Fb = f(b) ! F ! Fa* Fb * Fb must must be < 0 be < 0 DO IF (ABS(b-a) < ε) EXIT

p Note the red framed IF statement.

IF (ABS(b-a) < ε) EXIT c = (a+b)/2 IF (c ==a .OR. c == b) EXIT Fc = f(c)

Without this, an infinite loop can

Fc f(c) IF (Fc == 0) EXIT IF (Fa*Fc < 0) THEN b = c

• ccur if |b-a| is very

small and c can be a b d t

Fb = Fc ELSE a = c

• r b due to

rounding.

Fa = Fc END IF END DO

7

Bisection Method: 5/6 Bisection Method: 5/6

However, there is a catch. ABS(b-a)<ε is an , absolute error and the accuracy of the root is approximately k accurate digits if ε = 10-k. What if the actual root is smaller than ε = 10-k? In this case, the bisection method returns no accurate digit at all! Recall that –log10(|x-x*|/|x*|) is approximately P I P It! the number of accurate digits of x. Prove rove It! It! Thus, if the anticipated root is very close to 0,

• ne may change the absolute error to relative

error.

8

Bisection Method: 6/6 Bisection Method: 6/6

The ABS(b-a) < ε may be changed to S( ) y g ABS(b-a)/MIN(ABS(a),ABS(b)) < ε Thi t t h t ti l bl if th i iti l This test has a potential problem: if the initial bracket contains 0, MIN(ABS(a),ABS(b)) h 0 hi h ld di i i b can approach 0 which would cause a division by zero! Nothing is perfect!

9

Bisection Method Example: 1/3 Bisection Method Example: 1/3

Suppose we wish to find the root closest to 6 of pp the following equation:

i ( ( )) + f x x x x ( ) sin( cos( )) = + + 1

2

Plotting the equation f(x) = 0 around x = 6 is very helpful.

x + 1

very helpful. You may use gnuplot for plotting. Gnuplot is available free on many platforms such as is available free on many platforms such as Windows, MacOS and Unix/Linux.

10

Bisection Method Example: 2/3 Bisection Method Example: 2/3

Function f(x) has a root in [4,6], and f(4) < 0 and f(6) > 0. We may use [4,6] with the Bisection method.

f x x x x ( ) sin( cos( )) = + + 1

2

x + 1

11

Bisection Method Example: 3/3 Bisection Method Example: 3/3

19 iterations, x* ≈ 5.5441017… f(x*) ≈ 0.863277×10-7. , f( )

[5 5 5 625] [5.5,5.5625] f x x ( ) sin( cos( )) +

[5,6] [5.5,6]

[5.5,5.75] [5.5,5.625] f x x ( ) ( ( )) = + 1

2

12

[4,6] [5,6]

New ton’s Method: 1/13 New ton s Method: 1/13

Idea Idea: If a is very close to a root of y = f(x) , the Idea Idea: If a is very close to a root of y f(x) , the tangent line though (a, f(a)) intersects the X-axis at a point b that is hopefully closer to the root. at a point b that is hopefully closer to the root. The line through (a, f(a)) with slope f’(a) is:

f ( ) y f a x a f a − − = ( ) ( )

'

y=f(x) (a,f(a)) y=f(x) a b x*

13

a b x* y - f(a) = f’(a)×(x-a)

New ton’s Method: 2/13 New ton s Method: 2/13

The line through (a, f(a)) with slope f’(a) is The line through (a, f(a)) with slope f (a) is

y f a x a f a − = ( ) ( )

'

Its intersection point b with the X-axis can be

x a −

found by setting y to zero and solving for x:

b a f a = − ( )

(a f(a))

b a f a = − ( )

'

y=f(x) (a,f(a))

14

a b x*

New ton’s Method: 3/13 New ton s Method: 3/13

Starting with a x0 that is close to a root x*, Starting with a x0 that is close to a root x , Newton’s method uses the following to compute x1, x2, … until some xk converges to x*. x1, x2, … until some xk converges to x . x x f x

i i i

+ =

−

1

( )

'

f x

i i i

+1

( )

'

x*

15

y=f(x) x1 x0 x2 x3

New ton’s Method: 4/13 New ton s Method: 4/13

Conver Convergence 1/2: ence 1/2: g If there is a constant ρ, 0 < ρ < 1, such that

lim | |

* k

x x

+ −

1

ρ

the sequence x0, x1, …, is said to converge

lim | | | |

* k k k

x x

→∞ +

− =

1

ρ

1

linearly with ratio (or rate) ρ. If ρ is zero, the convergence is superlinear. If xk+1 and xk are close to x*, the above expression means | *| ρ| *| |xk+1 – x*| ≈ ρ|xk – x*| The error at xk+1 is proportional to (i.e., linear) and smaller than the error at x since

16

linear) and smaller than the error at xk since ρ < 1.

New ton’s Method: 5/13 New ton s Method: 5/13

Convergence Convergence 2/2: 2/2: Convergence Convergence 2/2: 2/2: If there is a p > 1 and a C > 0 such that | |

*

lim | | | |

* k k k p

x x x x C

→∞ + −

− =

1

the sequence is said to converge with order p. If xk+1 and xk are close to x*, the above If xk+1 and xk are close to x , the above expression yields |xk+1- x*| ≈ C|xk - x*|p. Since |x - x*| is close to 0 and p > 1 |x

• x*|

Since |xk - x | is close to 0 and p > 1, |xk+1-x | is even smaller, and has p more 0 digits.

17

New ton’s Method: 6/13 New ton s Method: 6/13

The right shows a The right shows a possible algorithm that implements Newton’s

x = initial value Fx = f(x) Dx = f’(x)

implements Newton s method. Newton’s method

f (x) DO IF (ABS(Fx) < ε) EXIT New_X = x – Fx/Dx

Newton s method converges with order 2. H it t

Fx = f(New_X) Dx = f’(New_X) x = New_X

However, it may not converge at all, and a i b f

END DO

maximum number of iterations may be needed t t l

18

to stop early.

New ton’s Method: 7/13 New ton s Method: 7/13

Problem Problem 1/6: 1/6: Problem Problem 1/6: 1/6: While Newton’s method has a fast convergence rate it may not converge at all no convergence rate, it may not converge at all no matter how close the initial guess is.

1/(2 1) n

y x

+

=

verify this fact yourself

y

x1 x2 x3 x4

19

New ton’s Method: 8/13 New ton s Method: 8/13

Problem 2/6: Problem 2/6: Newton’s method can oscillate. If y=f(x)=x3-3x2+x+3 then y’=f’(x)=3x2-6x+1 If y f(x) x 3x +x+3, then y f (x) 3x 6x+1. If the initial value is x0=1, then we will have x1 = 2, x2 = 1, x3 = 2, …. (see textbook page 78 for 2, x2 1, x3 2, …. (see textbook page 78 for a function plot). Note that if x0 = -1, in 3 iterations the root is Note that if x0 1, in 3 iterations the root is found at -0.76929265! Therefore, carefully choosing an initial guess Therefore, carefully choosing an initial guess is important.

20

New ton’s Method: 9/13 New ton s Method: 9/13

Problem 3/6: Problem 3/6: Newton’s method can oscillate around a minimal point, eventually causing f’(x) to be zero. In the following, x1, x3, x5, … approach the i i hil h i fi it minimum while x2, x4, …, approach infinity.

21

x1 x2 x3 x4 x5 x*

New ton’s Method: 10/13 New ton s Method: 10/13

Problem 4/6: Problem 4/6: Since Newton’s method uses xk+1 = xk - f(xk)/f’(xk), it requires f’(xk) ≠ 0. As a result, multiple roots can be a problem because f’(x)=0 at a multiple root. Inflection points may also cause problems.

f’(x) = 0 at this inflection point

f’(x )=0

f (x) = 0 at this inflection point

f (x2)=0

22

x1 x2 x*

f’(x) = 0 at a multiple root

New ton’s Method: 11/13 New ton s Method: 11/13

Problem 5/6: Problem 5/6: Newton’s method can converge to a remote root even though the initial guess is close to the ti i t d anticipated one.

x1 x2 x* x3 x4

1 3 4

23

This is what we want! But, we get this one

New ton’s Method: 12/13 New ton s Method: 12/13

Problem 6/6: Problem 6/6: Even though Newton’s method has a convergence rate of 2, it can be very slow if the i iti l i t i ht Th f ll i l initial guess is not right. The following solves x10-1=0. Note that this x* = 1 is a multiple root!

It ti 0 0 5 f( ) 0 99902343 Iteration 0 x = 0.5 f(x) = -0.99902343 Iteration 1 x = 51.65 f(x) = 1.3511494E+17 Iteration 2 x = 46.485 f(x) = 4.711166E+16 Iteration 3 x = 41.836502 f(x) = 1.6426823E+16 ( ) Iteration 4 x = 37.65285 f(x) = 5.727679E+15 ...... Other output ...... Iteration 10 x = 20.01027 f(x) = 1.0292698E+13 Oth t t ...... Other output ...... Iteration 20 x = 6.9771494 f(x) = 273388500.0 Iteration 30 x = 2.4328012 f(x) = 7261.167 Iteration 40 x = 1 002316 f(x) = 0 02340281

slow!

24

Iteration 40 x = 1.002316 f(x) = 0.02340281 Iteration 41 x = 1.000024 f(x) = 2.3961067E-4 Iteration 42 x = 1.0 f(x) = 0.0E+0

New ton’s Method: 13/13 New ton s Method: 13/13

A few important notes: Newton’s method is also referred to as the Newton-Raphson method. p Plot the function to find a good initial guess. When the computation converges plug the x* When the computation converges, plug the x* back to the original equation to make sure it is indeed a root is indeed a root. Check for f’(x) = 0. The use of a maximum number of iterations would help prevent infinite loop.

25

New ton’s Method Example: 1/3 New ton s Method Example: 1/3

Suppose we wish to find the root closest to 4 of the following equation:

f x e x

x

( ) sin( cos( )) = +

−

The following is a plot in [0,10]:

f x e x ( ) sin( cos( )) +

f x e x

x

( ) sin( cos( )) = +

−

26

New ton’s Method Example: 2/3 New ton s Method Example: 2/3

A plot in [4,6] shows that the desired root is approximately 4.7, and f(x) is monotonically increasing in [4,6]:

f x e x

x

( ) sin( cos( )) = +

−

27

New ton’s Method Example: 3/3 New ton s Method Example: 3/3

The following has f(x) and f’(x):

f x e x f x e x e x

x x x

( ) sin( cos( )) ( ) cos( cos( ))( sin( ))

'

= + = − + +

− − −

Initial x0 = 4. 3 iterations

x f(x) f’(x)

3 iterations x* ≈ 4.703324 f(x*) ≈ 0 8102506×10-7

4

• 0.59344154

0.5943911 1 4.998402 0.2848776 0.9131543

f(x*) ≈ 0.8102506×10 7 f’(x*) ≈ 0.99089384

2 4.686432

• 0.016733877

0.99030494 3 4.703329 0.5750917×10-5 0.99089395 4 4.703324 0.8102506×10-7 0.99089384

28

The Secant Method: 1/8 The Secant Method: 1/8

Idea Idea :Newton’s method requires f’(x). What if f’(x) is difficult to compute or even not available? Methods that use an approximation g of f’(x) Methods that use an approximation gk of f (x) are referred to as quasi-Newton methods. Thus, xk 1 is computed (from xk) as follows: xk+1 is computed (from xk) as follows:

x x f x g

k k k

+ =

−

1

( )

The secant method offers a simple way for

gk

estimating gk.

29

The Secant Method: 2/8 The Secant Method: 2/8

The secant method uses the slope of the chord b t d i t f f’( ) between xk-1 and xk as an approximate of f’(xk) . Therefore, xk+1 is the intersection point of this h d d th X i chord and the X-axis. Since the secant method uses two points, it is ll f d t i h d usually referred to as a two-point method.

(xk-1,f(xk-1)) (xk,f(xk))

30

x* xk-1 xk xk+1

Newton’s method

The Secant Method: 3/8 The Secant Method: 3/8

Since the slope gk is Since the slope gk is

g f x f x

k k k

= −

−

( ) ( )

1

Some algebraic manipulations yield the secant

g x x

k k k

−

−1

method formula:

1

( )

k k

x x x x f x

−

− = − ×

It is more complex than that of Newton’s method;

1 1

( ) ( ) ( )

k k k k k

x x f x f x f x

+ −

= − × −

It is more complex than that of Newton s method; but, it may be cheaper than evaluating f’(xk) .

31

The Secant Method: 4/8 The Secant Method: 4/8

Convergence: Convergence: Convergence: Convergence: The secant method, in fact all two-point methods are superlinear! methods, are superlinear! More precisely, one can show that the f ll i h ld ith 1 < < 2 d C > 0 following holds with a 1 < p < 2 and C > 0

lim | |

* k

x x C

+ −

=

1

Moreover, p is the root of p2 – p – 1 =0:

lim | |

* k k p

x x C

→∞

−

Moreover, p is the root of p p 1 0:

p = + = 1 5 2 1618 . .....

32

p 2

The Secant Method: 5/8 The Secant Method: 5/8

Algorithm: Algorithm:

a = initial value #1

Algorithm: Algorithm:

The right is the secant method.

b = initial value #2 Fa = f(a) Fb = f(b)

Note that Fb and Fa may be equal, causing

DO c = b – Fb*(b-a)/(Fb-Fa) Fc = f(c)

either overflow or division by zero!

IF (ABS(Fc) < ε) EXIT a = b Fa = Fb b

0/0 is also possible! Since there are three bt ti it

b = c Fb = Fc END DO

subtractions, rewrite the expression to avoid cancellation if possible

33

cancellation if possible.

The Secant Method: 6/8 The Secant Method: 6/8

The following shows the output of solving f(x) = The following shows the output of solving f(x) x10-1 with initial points x0=2 and x1=1.5. It is faster than Newton’s method (42 iterations) It is faster than Newton s method (42 iterations).

0 : 2.0 1023.0 1.5 56.66504

a f(a) b f(b)

1 : 1.5 56.66504 1.4706805 46.33511 2 : 1.4706805 46.33511 1.3391671 17.550165 3 : 1.3391671 17.550165 1.2589835 9.004613 4 : 1 2589835 9 004613 1 1744925 3 994619 4 : 1.2589835 9.004613 1.1744925 3.994619 5 : 1.1744925 3.994619 1.1071253 1.7667356 6 : 1.1071253 1.7667356 1.0537024 0.68725025 7 : 1.0537024 0.68725025 1.0196909 0.21530557 8 : 1.0196909 0.21530557 1.0041745 0.04253745 9 : 1.0041745 0.04253745 1.0003542 3.5473108E-3 10 : 1.0003542 3.5473108E-3 1.0000066 6.556511E-5

34

Converged after 11 iterations x* = 1.0 f(x*) = 0.0E+0

The Secant Method: 7/8 The Secant Method: 7/8

Due to subtractions, cancellation and loss of Due to subtractions, cancellation and loss of significant digits are possible. If we start with x = 0 and x = 1 5 division by If we start with x0 = 0 and x1 = 1.5, division by zero can occur! We hit a flat area of the curve.

l i h

Initially:

a = 0.0E+0 f(a) = -1.0 b = 1.5 f(b) = 56.66504

explain why

a 0.0E 0 f(a) 1.0 b 1.5 f(b) 56.66504 c = 0.026012301 f(c) = -1.0

Iteration 1 :

a = 1.5 f(a) = 56.66504 b = 0.026012301 f(b) = -1.0 c = 0.051573503 f(c) = -1.0

Iteration 2 :

a = 0.026012301 f(a) = -1.0 b = 0.051573503 f(b) = -1.0

Floa

• ating

ting exception ception (Oops!!!) (Oops!!!)

35

Fl Floating exception ception (Oops!!!) (Oops!!!)

The Secant Method: 8/8 The Secant Method: 8/8

A few important notes: A few important notes: The secant method evaluates the function

• nce per iteration
• nce per iteration.

Plot the function to find a good initial guess. Check if the computed result is indeed a root. Check for f(xk) – f(xk-1) ≈ 0 (i.e., flat area) as

k k 1

this can cause overflow or division by zero. The secant method shares essentially the The secant method shares essentially the same pitfalls with Newton’s method. The use of a maximum number of iterations

36

The use of a maximum number of iterations can help prevent infinite loop.

Handling Multiple Roots: 1/2 Handling Multiple Roots: 1/2

A polynomial f(x) has a root x* of multiplicity p p y f( ) p y p > 1, where p is an integer, if f(x) = (x – x*)pg(x) for some polynomial g(x) and g(x*)≠0. Then we have f(x*) = f’(x*) = f”(x*) = … f[p-1](x*) = f[p](x*) = 0 f( ) f ( ) f ( ) f ( ) f ( ) For example, f(x) = (x-1)(x-1)(x-2)(x-3) has a double root (i.e., p = 2) and f(x)=(x-1)2g(x) and ( , p ) f( ) ( ) g( ) g(x) = (x-2)(x-3). Newton’s method and the secant method, , converge only linearly for multiple roots.

37

Handling Multiple Roots: 2/2 Handling Multiple Roots: 2/2

Two modifications to Newton’s method can still Two modifications to Newton s method can still maintain quadratic convergence. If x* is a root of multiplicity p then use the If x is a root of multiplicity p, then use the following instead of the original:

( ) f

1 '

( ) ( )

k k k k

f x x x p f x

+ =

− ×

Or, let g(x) = f(x)/f’(x) and use the following:

( )

k

g x

1

( ) '( )

k k k k

g x x g x

+ =

−

38

The Regula Falsi Method: 1/5 The Regula Falsi Method: 1/5

Idea Idea: The Regula Falsi, aka false- positions, Idea Idea: The Regula Falsi, aka false positions, method is a variation of the bisection method. The bisection method does not care about the The bisection method does not care about the location of the root. It just cuts the interval at the mid-point and as a result it is slower the mid-point, and, as a result, it is slower. We may use the intersection point of the X-axis d th h d b t ( f( )) d (b f(b)) and the chord between (a, f(a)) and (b, f(b)) as an

• approximation. Hopefully, this would be faster.

39

The Regula Falsi Method: 2/5 The Regula Falsi Method: 2/5

The intersection point c

( f( ))

The intersection point c may not be closer to x* than the mid-point (a+b)/2 is.

(a,f(a))

the mid point (a b)/2 is. The line of (a, f(a)) and (b, f(b)) is

a b c (a+b)/2 x*

f(b)) is

a b c

y f a x a f b f a b a − − = − − ( ) ( ) ( )

Setting y=0 and solving for x yield c as follows:

(b,f(b))

c a b a f b f a f a b b a f b f a f b = − − − × = − − − × ( ) ( ) ( ) ( ) ( ) ( )

40

( f( ))

The Regula Falsi Method: 3/5 The Regula Falsi Method: 3/5

The Regula Falsi

Fa = f(a)

g method is a variation

• f the bisection

Fa = f(a) Fb = f(b) ! F ! Fa* Fb * Fb must must be < 0 be < 0 DO c = a-Fa*(b-a)/(Fb-Fa)

method. Instead of computing

c = a-Fa*(b-a)/(Fb-Fa) Fc = f(c) IF (ABS(Fc) < ε) EXIT IF (Fa*Fc < 0) THEN

the mid-point, it computes the i t ti i t f

IF (Fa Fc < 0) THEN b = c Fb = Fc ELSE

intersection point of the X-axis and the chord

a = c Fa = Fc END IF

chord. Otherwise, they are identical

END DO

41

identical.

The Regula Falsi Method: 4/5 The Regula Falsi Method: 4/5

The Regula Falsi method can be very slow when g y hitting a flat area. One of the two end points may be fixed, while p y , the other end slowly moves the root. Near a very flat area f(b) – f(a) ≈ 0 holds! y f( ) f( )

x2 x3 x4 x5 x* x1 x0

42

The Regula Falsi Method: 5/5 The Regula Falsi Method: 5/5

Solving x10-1=0 with a = 0 and b = 1.3. Very slow! Solving x 1 0 with a 0 and b 1.3. Very slow!

1 0.0E+0 -1.0 1.3 12.785842 0.094299644 -1.0

a f(a) b c f(b) f(c)

2 0.094299644 -1.0 1.3 12.785842 0.18175896 -0.99999994 3 0.18175896 -0.99999994 1.3 12.785842 0.26287412 -0.99999845 4 0.26287412 -0.99999845 1.3 12.785842 0.33810526 -0.99998044 5 0.33810526 -0.99998044 1.3 12.785842 0.4078781 -0.99987256 10 0.6395442 -0.98855262 1.3 12.785842 0.6869434 -0.97660034 20 0.94343614 -0.44136887 1.3 12.785842 0.95533406 -0.36678296 30 0.99553364 -0.043776333 1.3 12.785842 0.9965725 -0.03375119 40 0.9996887 -3.1086206E-3 1.3 12.785842 0.9997617 -2.3804306E-3 50 0.99997854 -2.1457672E-4 1.3 12.785842 0.99998354 -1.6450882E-4 60 0.99999856 -1.4305115E-5 1.3 12.785842 0.9999989 -1.0728836E-5 43 60 0.99999856 1.4305115E 5 1.3 12.785842 0.9999989 1.0728836E 5 63 0.9999995 -4.7683715E-6

Fixed-Point Iteration: 1/19 Fixed Point Iteration: 1/19

Rewrite f(x) = 0 to a new form x = g(x). Rewrite f(x) 0 to a new form x g(x).

f(x) = sin(x)-cos(x) = 0 becomes sin(x)=cos(x). Hence, x = sin-1(cos(x)) or x = cos-1(sin(x)). ( ( )) ( ( )) f(x) = e-x – x2 = 0 becomes e-x = x2. Hence, x = -LOG(x2)

• r x = SQRT(e-x).

The fixed-point iteration starts with an initial value x0 and computes xi+1 = g(xi) until xk+1 ≈ g(xk) for p

i+1

g( i)

k+1 g( k)

some k. We are seeking a x* such that x*=g(x*) holds We are seeking a x such that x g(x ) holds. This x* is a fixed-point of g() since function g() maps x* to x* itself and hence a root of f(x)

44

maps x to x itself, and, hence, a root of f(x).

Fixed-Point Iteration: 2/19 Fixed Point Iteration: 2/19

The fixed-point iteration x=g(x) can be The fixed point iteration x g(x) can be considered as computing the intersection point of the curve y=g(x) and the straight line y = x. the curve y g(x) and the straight line y x.

y = x y = g(x) x* = g(x*)

45

Fixed-Point Iteration: 3/19 Fixed Point Iteration: 3/19

Let us rewrite xk+1= g(xk)

k+1

g( k) as y = g(xk) and xk+1 = y . Thus, y = g(xk) means , y g( k) from xk find the y-value of g(xk).

x = y ( ) y = x Y

Then, xk+1 = y means use the y-value to find the

xk+1= y xk+2= y y = g(x)

new x-value xk+1 on y = x. Repeat this procedure as

y = g(xk) y = g(xk+1) X

shown until, hopefully, it converges.

xk xk+1 X

46

Fixed-Point Iteration: 4/19 Fixed Point Iteration: 4/19

Newton’s method uses xk+1 = xk – f(xk)/f’(xk). Let Newton s method uses xk+1 xk f(xk)/f (xk). Let g(x) = x – f(x)/f’(x) and Newton’s method becomes the fixed-point iteration: xk+1 = g(xk) . becomes the fixed point iteration: xk+1 g(xk) . The fixed-point iteration converges linearly if |g’(x)| < 1 in the neighborhood of the correct root |g (x)| < 1 in the neighborhood of the correct root. Observation: if g’(x)∈(0,1) and g’(x)∈(-1,0), the i t ti d ill t convergence is asymptotic and oscillatory, respectively, and it is faster if g’(x) approaches 0.

47

Fixed-Point Iteration: 5/19 Fixed Point Iteration: 5/19

Two convergence cases: Two convergence cases: Left g’(x) ∈ (-1,0) and hence oscillatory Ri ht ’( ) (0 1) d t ti Right g’(x) ∈ (0,1) and asymptotic

’( ) ∈ ( 1 0) ’( ) ∈ (0 1) g’(x) ∈ (-1,0) g’(x) ∈ (0,1)

x1 x2 x3 x1 x2 x3

48

x1 x2 x3 x1 x2 x3

• scillatory

asymptotic

Fixed-Point Iteration: 6/19 Fixed Point Iteration: 6/19

Two divergence cases: Two divergence cases: Left g’(x) < -1 Ri ht ’( ) > 1 Right g’(x) > 1

g’(x) > 1 g’(x) < -1

49

x1 x2 x3 x4 x1 x2 x3

Fixed-Point Iteration: 7/19 Fixed Point Iteration: 7/19

Since f(x) = 0 has to be transformed to a form of Since f(x) 0 has to be transformed to a form of x = g(x) in order to use the fixed-point iteration, algebraic manipulation is necessary. algebraic manipulation is necessary. Incorrect transformations can yield Incorrect transformations can yield roots roots w e don don’t w ant ant diverge diverge or roots roots w e w e don don t w ant w ant, , diverge diverge, , or

• r

cause runtime errors. cause runtime errors. A lt ft t f ti it i As a result, after transformation, it is very helpful to plot the curve of y = g(x) and the line y d i t th l ’( ) th d i d = x, and inspect the slope g’(x) near the desired fixed point to find a good initial guess.

50

Fixed-Point Iteration: 8/19 Fixed Point Iteration: 8/19

Consider f(x) = cos(x) – x1/2. Consider f(x) cos(x) x .

a root in the range of [0,2] f(x) = cos(x) – x1/2

51

Fixed-Point Iteration: 9/19 Fixed Point Iteration: 9/19

Setting f(x) to zero yields Setting f(x) to zero yields R i t i ld

cos( ) x x − =

Rearranging terms yields cos2(x) = x.

2

y = cos2(x)

We use g(x) = cos2(x) for the fixed-point iteration. From curve y = cos2(x) and line y = x, we see only

y = x

0 6 0 8

• ne root at .64171…..

approximately .

0.6 0.8

52

pp y

Fixed-Point Iteration: 10/19 Fixed Point Iteration: 10/19

Setting f(x) to zero yields g f( ) y Rearranging terms yields

cos( ) x x − =

Rearranging terms yields cos(x) = x1/2. Taking cos-1 yields g(x) =

y = cos-1(x1/2)

Taking cos yields g(x) cos-1(x1/2). From y = cos-1(x1/2) and From y cos (x ) and line y = x, we see only one root at 0.64171…..

y = x

0 6 0 8

approximately.

0.6 0.8

53

Fixed-Point Iteration: 11/19 Fixed Point Iteration: 11/19

One may transform f(x) in many ways, each of One may transform f(x) in many ways, each of which could yield a different root. Some may diverge or cause runtime errors. diverge or cause runtime errors. Again, plotting y = g(x) and y = x and inspecting the slope g’(x) can help determine the the slope g (x) can help determine the transformation that can yield the best result. F l th f ll i f( ) h t b i For example, the following f(x) has two obvious transformations:

2

( ) log( 1) 1 f x x x = + − +

54

Fixed-Point Iteration: 12/19 Fixed Point Iteration: 12/19

Plotting this function f(x) shows that there are Plotting this function f(x) shows that there are roots in [-0.8, -0.7], [2.0, 3.0] and [70, 80].

55

Fixed-Point Iteration: 13/19 Fixed Point Iteration: 13/19

If we transform f(x) in the following way: If we transform f(x) in the following way:

2

log( 1) 1 x x + +

2 1

log( 1) 1 1

x

x x x e

+

+ = + + =

2 1

1 1

x

x e x e

+

+ = −

1

1

x

x e

+

= −

1

( ) 1

x

g x e

+

= −

56

Fixed-Point Iteration: 14/19 Fixed Point Iteration: 14/19

The fixed-point iteration converges to x* = The fixed point iteration converges to x 2.2513…… easily since 0 < g’(x) < 1; but, it is difficult to get the other two roots. difficult to get the other two roots.

This is the asymptotic case y = x

1

( ) 1

x

g x e

+

= −

57

Fixed-Point Iteration: 15/19 Fixed Point Iteration: 15/19

If f(x) is transformed in a different way: If f(x) is transformed in a different way:

2

log( 1) 1 x x + = +

( )

2 2

log( 1) 1 x x + = +

( )

2 2 2

log( 1) 1 x x = + −

( )

2 2

( ) log( 1) 1 g x x = + −

58

Fixed-Point Iteration: 16/19 Fixed Point Iteration: 16/19

This is also the asymptotic case

( )

2

( )

( )

2 2

( ) log 1 1 g x x = + −

If the initial guess is larger than x* = 2.2513… the fixed-point iteration converges to 72.309…

59

Fixed-Point Iteration: 17/19 Fixed Point Iteration: 17/19

This is the oscillatory case If th i iti l i l th * 2 2513 If the initial guess is less than x* = 2.2513… the fixed-point iteration converges to -0.7769…

( )

( )

2 2

( ) log 1 1 g x x = + −

60

( )

( )

( ) g g

Fixed-Point Iteration: 18/19 Fixed Point Iteration: 18/19

( )

( )

2 2

( ) log 1 1 g x x = + −

no way to reach x* = 2.2513… converge to 72.309… converge to -0.7769…

61

Fixed-Point Iteration: 19/19 Fixed Point Iteration: 19/19

Conclusion: The fixed-point iteration is easy to Conclusion: The fixed point iteration is easy to

• use. But, different transformations

But, different transformations may may yield ield different different roots roots and and may may yield yield different different roots roots and and sometimes sometimes may may not converge. not converge. Moreo eover, so some transform rmat ations m may ay

• eo e , so

e t a s o at at o s ay ay lead to runtime errors. lead to runtime errors. One may try different transformations to find One may try different transformations to find the desired root. Only those transformations that can compute the desired root are considered that can compute the desired root are considered correct. U th fi d i t it ti ith

62

Use the fixed-point iteration with care.

Fixed-Point Iteration Example: 1/9 Fixed Point Iteration Example: 1/9

Consider the following function f(x): Consider the following function f(x): Th f ll i l t h th t f( ) 0 h f

f x e x

x

( ) sin( ) = −

−

The following plot shows that f(x) = 0 has four roots in [0,10]. In fact, f(x) = 0 has an infinite b f t number of roots.

f x e x

x

( ) sin( ) = − =

−

63

Fixed-Point Iteration Example: 2/9 Fixed Point Iteration Example: 2/9

Since e-x – sin(x) = 0, we have sin(x) = e-x, and x = Since e sin(x) 0, we have sin(x) e , and x sin-1(e-x). Therefore, we may use g(x) = sin-1(e-x). We have a problem We have a problem.

y = g(x) and and y=x ha have onl

• nly

y one

• ne

y x =

inter intersection point ection point ar around

• und 0.5.

0.5. As a r As a result, sult, onl

• nly

y the r the root in

• ot in

[0,1] can be f can be found w ith

• und w ith the

the [ , ] fix fixed-point iter d-point iteration. ion.

y g x e x = =

− −

( ) sin ( )

1 64

0.5885…

Fixed-Point Iteration Example: 3/9 Fixed Point Iteration Example: 3/9

Since 0 ≤ e-x ≤ 1 g(x) = sin-1(e-x) can be Since 0 ≤ e ≤ 1, g(x) sin (e ) can be computed without problems. Since e-x decreases monotonically from 1 = e-0 to 0 as x→∞, sin-1(e-x) also decreases monotonically from π/2 = sin-1(e-0) = sin-1(1) to 0 as x→∞. to 0 as x→ . Therefore, g(x) = sin-1(e-x) and y = x has l i t ti i t d th fi d

• nly one intersection point, and the fixed-

point iteration can only find one root of f(x)

65

= 0.

Fixed-Point Iteration Example: 4/9 Fixed Point Iteration Example: 4/9

Left is a plot of g(x) in [0,1]. p g( ) [ , ] g’(x) is calculated below:

g x e x ( ) sin ( ) =

− −

1

x x (sin ( ))' = −

−1

2

1 1

Si ’(0 4) 0 90 ’(0 5)

g x e e

x x

( )

'

= − −

− −2

1

Since g’(0.4) = -0.90…, g’(0.5) = -0.76… , g’(0.6) = -0.66…, |g’(x)| < 1 around the root |g’(x)| < 1 around the root 0.58855… and the fixed- point iteration converges.

66

point iteration converges.

Fixed-Point Iteration Example: 5/9 Fixed Point Iteration Example: 5/9

xi 0.5 1 0.6516896 2 0 5482148 2 0.5482148 3 0.616252 4 0 5703949 4 0.5703949 5 0.6007996

26 it ti 26 iterations x* ≈ 0.58853024

67

x0 x1 x2 x3 x4 x5

Fixed-Point Iteration Example: 6/9 Fixed Point Iteration Example: 6/9

We may transform e-x – sin(x) = 0 differently. y ( ) y Since e-x – sin(x) = 0, we have e-x = sin(x). Taking natural logarithm yields -x = ln(sin(x)) g g y and hence x = -ln(sin(x)) . We may use g(x) = -ln(sin(x)). However, this approach has two problems: If sin(x) = 0, ln(sin(x)) is -∞. If sin(x) < 0, ln(sin(x)) is not defined. Since sin(x) is periodic, ln(sin(x)) is undefined on i fi i b f i l (i [ 2 ] an infinite number of intervals (i.e., [π,2π], [3π,4π], [5π,6π], …). In general, ln(sin(x)) is undefined on [(2n 1)π 2nπ] where n ≥ 1

68

undefined on [(2n-1)π, 2nπ], where n ≥ 1.

Fixed-Point Iteration Example: 7/9 Fixed Point Iteration Example: 7/9

A plot shows y = g(x) and y = x at two points in (0,π). There are other roots not in (0,π).

y x = y x =

Ob i l th d i ti Obviously, the derivative at this larger root is larger than 1. Hence, the fixed point iteration won’t

g x x ( ) ln(sin( )) = −

fixed-point iteration won’t be able to find it!

x = π x = 0

g( ) ( ( ))

Ex Exer ercis cise: Plot g(x) in [0 20] to see other part

69

[0,20] to see other part.

Fixed-Point Iteration Example: 8/9 Fixed Point Iteration Example: 8/9

Let us examine the smaller root in [0.4,0.6] . Let us examine the smaller root in [0.4,0.6] . Since g(x) = -ln(sin(x)), g’(x) is (use chain rule):

( ) g x x x x x

' '

( ) ( ln(sin( )) cos( ) sin( ) cot( ) = − = − = −

We have g’(0.4) = -2.36…, g’(0.5) = -1.83…, g’(0.6) = -1.46…. Since |g’(x)| > 1 around the root, the fixed-point iteration diverges and is not be able to find the iteration diverges and is not be able to find the desired root.

70

Fixed-Point Iteration Example: 9/9 Fixed Point Iteration Example: 9/9

If we start with x0 = 0.4, we have x1 = 0.9431011, If we start with x0 0.4, we have x1 0.9431011, x2 = 0.2114828, x3 = 1.561077, x4 = 0.0000472799, x5 = 9.960947 x5 9.960947 Since sin(x5) = sin(9.960947) = -0.51084643, g(x5) = g(9 960947) = ln(sin(-0 51084643)) is not well- = g(9.960947) = ln(sin(-0.51084643)) is not well- defined, and the fixed-point iteration fails. Hence ( ) ln(sin( )) is NOT NOT a correct Hence, g(x) = -ln(sin(x)) is NOT NOT a correct transformation. You may follow this procedure to do a simple analysis before using the fixed-point iteration.

71

The End The End

72