Fixed-Point Iteration
Source
For the given equation:
f(x)=0⟹x=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
Convergence check (or analysis)
Let r be the exact solution, r=g(r)
xk+1=g(xk)error:ek=xk−rek+1=xk+1−r=g(xk)−g(r)=g′(η)(xk−r)=g′(η)ek,η∈(xk,r)⟹|ek+1|≤|g′(η)||ek|
If |g′(η)|<1,error decreases (iteration converges)
1) Example of x=cos(x)
# Computational Thinking on how to calculate cos(x) = x
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
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()
# 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))))))))
# better way
x = 10
for i in range(24):
x = np.cos(x)
print (x)
# 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)
2) example of √2(x2=2)
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()
# Use an idea of a fixed point
x = 2
for i in range(10):
x = 2/x
print (x)
n_iter = 10
x = np.linspace(-3, 3, 100)
y = 2/x
def func(x):
return 2/x
x_nu = np.zeros((n_iter,1))
x_nu[0] = 3
for i in range(n_iter-1):
x_nu[i+1] = 2/x_nu[i]
traj_x = []
traj_y = []
for i in range(n_iter):
traj_x.append(x_nu[i])
traj_y.append(func(x_nu[i]))
traj_x.append(func(x_nu[i]))
traj_y.append(func(x_nu[i]))
plt.figure(figsize = (8,6))
plt.plot(x, y, linewidth = 2)
plt.plot(x, x, linewidth = 2)
plt.plot(traj_x,traj_y,'r--')
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()
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()
# 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)
n_iter = 10
x = np.linspace(-4, 4, 100)
y = (x + 2/x)/2
def func(x):
return (x + 2/x)/2
x_nu = np.zeros((n_iter,1))
x_nu[0] = 4
for i in range(n_iter-1):
x_nu[i+1] = (x_nu[i] + 2/x_nu[i])/2
traj_x = []
traj_y = []
for i in range(n_iter):
traj_x.append(x_nu[i])
traj_y.append(func(x_nu[i]))
traj_x.append(func(x_nu[i]))
traj_y.append(func(x_nu[i]))
plt.figure(figsize = (8,6))
plt.plot(x, y, linewidth = 2)
plt.plot(x, x, linewidth = 2)
plt.plot(traj_x,traj_y, 'r--')
plt.axvline(x=0, color = 'k', linestyle = '--')
plt.axhline(y=0, color = 'k', linestyle = '--')
plt.axis('equal')
plt.ylim([-1,1])
plt.show()
Think about why it gives different results.
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)
This solution only possible for small size problems. There are many iterative methods for large problems.
# 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)
Remark) try this one
x1=−3x1+x2−x3+7x2=4x1−7x2+x3+21x3=2x1−x2+−4x3+15Convergence check
x←Ax+b
# 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)
# stability, check eigenvalue of A
A = np.array(([[3, 1, -1 ],
[4, 7, 1],
[2, -1, -4]]))
np.linalg.eig(A)
# 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)
%%html
<center><iframe src="https://www.youtube.com/embed/xD7OCAjJ4TQ?rel=0"
width="420" height="315" frameborder="0" allowfullscreen></iframe></center>
%%html
<center><iframe src="https://www.youtube.com/embed/EPhPj40wBs4?rel=0"
width="420" height="315" frameborder="0" allowfullscreen></iframe></center>
%%javascript
$.getScript('https://kmahelona.github.io/ipython_notebook_goodies/ipython_notebook_toc.js')