Processing math: 100%



Fixed-Point Iteration


By Prof. Seungchul Lee
http://iailab.kaist.ac.kr/
Industrial AI Lab at KAIST

Source

  • By Prof. David J. Malan from Harvard University
  • By Prof. Erik Demaine from MIT online lecture

1. Fixed-Point Iteration

1.1. Numerical approach

For the given equation:

f(x)=0x=g(x)

Remark: always achievable

x=f(x)+x=g(x)x=f(x)+x=g(x)

Goal: numerically find the solution of x=g(x)

Main idea:

  • Make a guess of the solution, xk

  • If g(xk) is 'nice', then hopefully, g(xk) will be closer to the answer. If so, we can iterate

Iteration algorithm

  1. choose an initial point x0
  2. Do the iteration xk+1=g(xk) until meeting stopping criteria

Convergence check (or analysis)

Let r be the exact solution, r=g(r)

xk+1=g(xk)error:ek=xkrek+1=xk+1r=g(xk)g(r)=g(η)(xkr)=g(η)ek,η(xk,r)|ek+1||g(η)||ek|


If |g(η)|<1,error decreases (iteration converges)

1) Example of x=cos(x)

In [1]:
# Computational Thinking on how to calculate cos(x) = x

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
In [2]:
x = np.linspace(-2*np.pi, 2*np.pi, 100)
y = np.cos(x)

plt.figure(figsize = (8,6))
plt.plot(x, y, linewidth = 2)
plt.plot(x, x, linewidth = 2)
# plt.xlim(-2*np.pi, 2*np.pi)
plt.axvline(x=0, color = 'k', linestyle = '--')
plt.axhline(y=0, color = 'k', linestyle = '--')
plt.legend(['cos(x)','x'])
plt.axis('equal')
plt.ylim([-1,1])
plt.show()
No description has been provided for this image
In [3]:
# naive approach

x = 0.3
print (np.cos(x))
print (np.cos(np.cos(x)))
print (np.cos(np.cos(np.cos(x))))
print (np.cos(np.cos(np.cos(np.cos(x)))))
print (np.cos(np.cos(np.cos(np.cos(np.cos(x))))))
print (np.cos(np.cos(np.cos(np.cos(np.cos(np.cos(x)))))))
print (np.cos(np.cos(np.cos(np.cos(np.cos(np.cos(np.cos(x))))))))
0.955336489125606
0.5773340444711864
0.8379206831271269
0.6690097308223832
0.7844362247423562
0.7077866472756374
0.7598027552852303
In [4]:
# better way

x = 0.3

for i in range(24):
    tmp = np.cos(x)
    x = np.cos(tmp)
    
print (x)
0.7390851312923253
In [5]:
# better way

x = np.zeros((24, 1))
x[0] = 0.3

for i in range(23):
    x[i+1] = np.cos(x[i])
    
print (x)
[[0.3       ]
 [0.95533649]
 [0.57733404]
 [0.83792068]
 [0.66900973]
 [0.78443622]
 [0.70778665]
 [0.75980276]
 [0.72497188]
 [0.74851807]
 [0.73269821]
 [0.74337234]
 [0.73619044]
 [0.74103194]
 [0.73777234]
 [0.73996881]
 [0.73848959]
 [0.73948617]
 [0.73881493]
 [0.73926712]
 [0.73896253]
 [0.73916771]
 [0.73902951]
 [0.7391226 ]]
In [6]:
# better way

x = 10
for i in range(24):
    x = np.cos(x)

print (x)
0.7390735444682907

2) example of 2(x2=2)

In [7]:
# Use an idea of a fixed point

x = 2
for i in range(10):
    x = 2/x
    print (x)
1.0
2.0
1.0
2.0
1.0
2.0
1.0
2.0
1.0
2.0
In [8]:
x = np.linspace(-3, 3, 100)
y = 2/x

plt.figure(figsize = (8,6))
plt.plot(x, y, linewidth = 2)
plt.plot(x, x, linewidth = 2)
plt.axvline(x=0, color = 'k', linestyle = '--')
plt.axhline(y=0, color = 'k', linestyle = '--')
plt.legend(['2/x','x'])
plt.axis('equal')
plt.ylim([-1,1])
plt.show()
No description has been provided for this image
In [9]:
# How to overcome 
# Use an idea of a fixed point +   kind of *|damping|*

x = 3
for i in range(10):
    x = (x + 2/x)/2
    print (x)
1.8333333333333333
1.4621212121212122
1.414998429894803
1.4142137800471977
1.4142135623731118
1.414213562373095
1.414213562373095
1.414213562373095
1.414213562373095
1.414213562373095
In [10]:
x = np.linspace(-4, 4, 100)
y = (x + 2/x)/2

plt.figure(figsize = (8,6))
plt.plot(x, y, linewidth = 2)
plt.plot(x, x, linewidth = 2)
plt.axvline(x=0, color = 'k', linestyle = '--')
plt.axhline(y=0, color = 'k', linestyle = '--')
plt.axis('equal')
plt.ylim([-1,1])
plt.show()
No description has been provided for this image

Think about why it gives different results.

1.2. System of Linear Equations


4x1x2+x3=74x18x2+x3=212x1+x2+5x3=15

In [11]:
import numpy as np
import matplotlib.pyplot as plt

# matrix inverse

A = np.array([[4, -1, 1], [4, -8, 1], [-2, 1, 5]])
b = np.array([[7, -21, 15]]).T

x = np.linalg.inv(A).dot(b)

print (x)
[[2.]
 [4.]
 [3.]]

This solution only possible for small size problems. There are many iterative methods for large problems.

4x1x2+x3=7x1=14x214x3+744x18x2+x3=21x2=12x1+18x3+2182x1+x2+5x3=15x3=25x115x2+155


  • In a matrix form

[x1x2x3]=[014141201825150][x1x2x3]+[742183]


  • iteration

[x(k+1)1x(k+1)2x(k+1)3]=[014141201825150][x(k)1x(k)2x(k)3]+[742183]

In [12]:
# Iterative way

A = np.array(([[0, 1/4, -1/4 ], 
               [4/8, 0, 1/8], 
               [2/5, -1/5, 0]]))
b = np.array([[7/4, 21/8, 15/5]]).T

# initial point
x = np.array([[1, 1, 2]]).T

A = np.asmatrix(A)
b = np.asmatrix(b)
x = np.asmatrix(x)

for i in range(20):
    x = A*x + b
    
print (x)
[[2.]
 [4.]
 [3.]]

Remark) try this one

x1=3x1+x2x3+7x2=4x17x2+x3+21x3=2x1x2+4x3+15


[x1x2x3][311471214][x1x2x3]+[72115]


Convergence check


xAx+b

xk+1=Axk+b=A(Axk1+b)+b=A2Xk1+Ab+b=Ak+1x0+Akb++Ab+b

In [13]:
# think about why this one does not work

A = np.array(([[3, 1, -1 ], 
               [4, 7, 1], 
               [2, -1, -4]]))
b = np.array([[7, 21, 15]]).T

# initial point
x = np.array([[1, 2, 2]]).T

for i in range(10):
    x = A.dot(x) + b
    
print (x)
[[1076845340]
 [ 660465147]
 [-237408283]]
In [14]:
# stability, check eigenvalue of A

A = np.array(([[3, 1, -1 ], 
               [4, 7, 1], 
               [2, -1, -4]]))

np.linalg.eig(A)
Out[14]:
(array([ 7.82147369,  1.65496639, -3.47644008]),
 array([[-0.21213573, -0.71387137,  0.17464301],
        [-0.97612458,  0.60136269, -0.15942549],
        [ 0.04668226, -0.3588183 ,  0.97163951]]))
In [15]:
# stability, check eigenvalue of A

A = np.array(([[0, 1/4, -1/4 ], 
               [4/8, 0, 1/8], 
               [2/5, -1/5, 0]]))

np.linalg.eig(A)
Out[15]:
(array([ 0.33471648+0.j        , -0.16735824+0.28987297j,
        -0.16735824-0.28987297j]),
 array([[-0.52873647+0.j        , -0.10993842+0.35839967j,
         -0.10993842-0.35839967j],
        [-0.83865529+0.j        ,  0.40608435-0.36739725j,
          0.40608435+0.36739725j],
        [-0.13074806+0.j        ,  0.74804945+0.j        ,
          0.74804945-0.j        ]]))
In [16]:
%%javascript
$.getScript('https://kmahelona.github.io/ipython_notebook_goodies/ipython_notebook_toc.js')