Fixed Point Iteration

Fixed Point Iteration

edit

What is Fixed Point Iteration?

edit

http://en.wikipedia.org/wiki/Fixed_point_iteration

Algorithm

edit

To find a solution to p = g(p) given an initial approximation p0.

INPUT: Initial approximation p0; tolerance TOL; maximum number of iterations N0. OUTPUT: approximate solution p or message of failure.

Step 1: Set i =1

Step 2: While   Steps 3 - 6

Step 3: Set p = g(p0)

Step 4: If   then OUTPUT (p); (The procedure was successful). STOP.

Step 5: Set i = i + 1

Step 6: Set p0 = p. (Update p0)

Step 7: OUTPUT ('The method failed after N0 iterations'). STOP.

Examples

edit

1. Find square root of 2 accurate till third decimal (10-3). The initial point is 2 .

Answer.

We know that,
x0 = 2
xn+1 = 0.5 * ( x0 + (2/ x0))
x1 = 0.5 * (2 + (2/2)) = 1.5
x2 = 0.5 * (1.5 + (2/1.5)) = 1.416
x3 = 0.5 * (1.416 + (2/1.416)) = 1.4142
We found out the square root of 2 accurate till 4th decimal in 3 steps.
As implemented in python
edit

The following python code implements the functionality of this section.

The two statements under "# Original estimate of square root (initial approximation.)" provide a pretty good estimate of initial approximation of   Unfortunately they are computationally intensive and add significant time to completion of function simpleSquareRoot().

The code within the while loop is both very simple and very fast.

Tolerance TOL is determined by the desired precision of the sqRoot, decimal.getcontext().prec, set by the code that invokes this function.

# python code.

import decimal
D = decimal.Decimal

def simpleSquareRoot (N) :
    '''
sqRoot = simpleSquareRoot (N)
Input N must be type Decimal or convertible to type Decimal.
output sqRoot is normally type Decimal.
output may be type None.
'''
    if type(N) != D :
        status = 0
        try : N = D(str(N))
        except: status = 1
        if status :
            print ('simpleSquareRoot (N): Error converting input to type Decimal.')
            return None

    if N in (0,1) : return N
    if N < 0 :
        print ('simpleSquareRoot (N) : input must be >= 0.')
        return None

    originalPrecision = decimal.getcontext().prec
    decimal.getcontext().prec += 2

    # Original estimate of square root (initial approximation.)
    # These next 2 statements are computationally intensive.
    # They add significant computational time to the execution of this function.
    sign, digits, exponent = N.as_tuple()
    n_ = n = D(10) ** ( (len(digits) + exponent) >> 1 )

    last = 0; half = D('0.5')

    if not simpleSquareRootDebug :
        while last != n :
            last = n
            n = half * (n + (N/n))

        decimal.getcontext().prec = originalPrecision
        sqRoot = (n+0).normalize()
        # Normal exit.
        return sqRoot

The following optional code is valid in debugging mode. It provides information about what happens during the while loop and verifies that the calculated square root is in fact within tolerance.

    # Here simpleSquareRootDebug is True.
    print ('simpleSquareRoot (N): N =',N)
    count = 0
    while last != n :
        count += 1
        last = n
        n = half * (n + (N/n))
        print ('n =',n)
    # Function simpleSquareRoot () checks itself.
    # Initial estimate compared with calculated value.
    ratio = decimal.Context(prec=5).create_decimal(n_/n)
    if 10 > ratio > D('0.1') :
        print ('    Original estimate =',n_,', ratio =',  ratio )
    else :
        print ('    Original estimate =',n_,', ratio =',  ratio, '*excessive*.' )
    print ('    count =',count)
    tolerance = D('1e-' + str(originalPrecision + 1))
    N_ = n * n
    v1 = abs(N-N_) ; v2 = abs(N+N_)/2
    status = (v1 > tolerance*v2)
    decimal.getcontext().prec = originalPrecision
    if status :
        print ('    Internal error 1.')
        return None
    sqRoot = (n+0).normalize()
    print ('    sqRoot =',sqRoot)
    return sqRoot

Compared with python's method decimal.Decimal.sqrt(), this code is slow for small values of precision. However, for numbers containing 550 decimal digits or more, this code (without debugging) is a little faster than python's method.

sqrt(2) with precision of 5
edit
simpleSquareRootDebug = 1

precision = decimal.getcontext().prec = 5

sqRoot = simpleSquareRoot(2)
simpleSquareRoot (N): N = 2
n = 1.5
n = 1.416666
n = 1.414216
n = 1.414214
n = 1.414214
    Original estimate = 1 , ratio = 0.70711
    count = 5
    sqRoot = 1.4142

sqrt(2) with 1,000 places of decimals
edit
precision = decimal.getcontext().prec =	1001

sqRoot = simpleSquareRoot(2)
simpleSquareRoot (N): N = 2
n = 1.5
n = 1.41666666666666......66666666
n = 1.41421568627450......92156862
n = 1.41421356237468......91918986
n = 1.41421356237309......93557326
n = 1.41421356237309......39275620
n = 1.41421356237309......87665231
n = 1.41421356237309......76778142
n = 1.41421356237309......02485270
n = 1.41421356237309......23938856
n = 1.41421356237309......48847209
n = 1.41421356237309......48847209
    Original estimate = 1 , ratio = 0.70711
    count = 12
    sqRoot = 
1.414213562373095048801688724209698078569671875376948073176679737990732478462
  107038850387534327641572735013846230912297024924836055850737212644121497099
  935831413222665927505592755799950501152782060571470109559971605970274534596
  862014728517418640889198609552329230484308714321450839762603627995251407989
  687253396546331808829640620615258352395054745750287759961729835575220337531
  857011354374603408498847160386899970699004815030544027790316454247823068492
  936918621580578463111596668713013015618568987237235288509264861249497715421
  833420428568606014682472077143585487415565706967765372022648544701585880162
  075847492265722600208558446652145839889394437092659180031138824646815708263
  010059485870400318648034219489727829064104507263688131373985525611732204024
  509122770022694112757362728049573810896750401836986836845072579936472906076
  299694138047565482372899718032680247442062926912485905218100445984215059112
  024944134172853147810580360337107730918286931471017111168391658172688941975
  8716582152128229518488472

Result was produced with only 12 passes through loop.

2. Consider the iteration pn+1 = g(p0) when the function g(x) = 1 + x - x2/4 is used. The fixed points can be found by solving equations x = g(x). The two solutions (fixed points of g) are x = -2 and x = 2. The derivative of the function is g'(x) = 1 - x/2, and there are only two cases to consider

Case 1 (P = -2):
Start with P0 = -2.05
We get
P1 = -2.100625
P2 = -2.20378135
.....  
Therefore, if |g'(x)| > 1 then sequence will not converge to P = -2
Case 2 (P = 2):
Start with P0 = 1.6
We get
P1 = 1.96
P2 = 1.9996
..... 2
Therefore, if |g'(x)| < 1 then sequence will converge to P = 2

Applications

edit

Two methods in which Fixed point technique is used:

1. Newton Raphson Method

Formula
 
where,
xn - initial point
f(xn) is the value of the function at that point
f'(xn) is the value of the differentiated function at that point.
Plug all these values into the above equation to get xn+1. It becomes the next initial point. Repeat until you get a point within an acceptable degree of error

2. Secant Method

Used to avoid differentiated form in Newton Raphson's method. Only problem is you need two initial points for this method (xn and xn-1)
Formula
 
Similar to Newton Raphson's method plug in all values to generate next approximation.

Exercises

edit

If 'f' is continuous and 'x' is a fixed point on 'f' then what is 'f(x)'?


Find the square root of 0.5 using fixed point iteration? Initial point 0.1 and tolerance 10-3

Calculating x1 i.e. 1st Iteration

Calculating x2 i.e. 2nd Iteration

Calculating x3 i.e. 3rd Iteration

Calculating x4 i.e. 4th Iteration

Calculating x5 i.e. 5th Iteration

In above problem, how many iterations or steps are needed to achieve an accuracy of 10-8

Quizzes

edit

<quiz display=simple> {Is there any possibility that the fixed point isn't unique |type="()"} -Yes +No

{ Assuming g ε C[a,b], If the range of the mapping y = g(x) satisfies y ε [a,b], then} +g has a fixed point in [a,b] -g has no fixed point in [a,b] -Depends on some other condition

{ If |g'(x)| > 1, then the iteration xn+1 = g (x) produces a sequence } +that diverges away from P -that converges towards P -Depends

{The fixed point iteration will diverge unless x0 = 0. |type="()"} +True -False

{Consider the following fixed point iteration;  . Rate (order) of convergence of this iteration is |type="()"} +Quadratic -Linear -Cubic