5. Iterative Computation#
Many algorithms involve iterative computation until a solution is found with a desired accuracy.
Also many mathematical models are formulated as a mapping from the current state to the next state, which gives discrete-time dynamics.
References:
Python Tutorial chapter 4: Control Flow Tools
Wikipedia: Newton’s method, Logistic map
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
Newton’s method#
For a linear equation \(A x = b\), or \(A x - b = 0\), the solution is given using the inverse matrix as
\[x = A^{-1} b. \]
For a general nonlinear equation \(f(x) = 0\), the solution may require iterative approximation. A typical way is by the Newton’s method:
\[x_{k+1} = x_k - \frac{f(x_k)}{f'(x_k)} \]
starting from an initial guess \(x_0\).
Example#
Let us define a plynomial function and its derivative.
def poly(x, a):
"""Polynomial a[0] + a[1]*x + ... + a[n]*x**n"""
n = len(a) # order+1
xp = np.array([x**i for i in range(n)]) # 1, x,.. x**(n-1)
y = a @ xp
return y
def d_poly(x, a):
"""Derivative of polynomial a[1] + 2*a[2]*x + ... + n*a[n]*x**(n-1)"""
b = np.arange(len(a))*a
return poly(x, b[1:])
a = np.array([3,2,1])
poly(1, a)
np.int64(6)
d_poly(1., a)
np.float64(4.0)
Here is an example polynomial with three zero-crossing points.
# f(x) = (x-3)(x-1)(x+2) = 6 - 5x -2x^2 + x^3
a = np.array([6, -5, -2, 1])
print('a =', a)
x = np.linspace(-3, 5, 50)
plt.plot(x, poly(x,a), x, d_poly(x,a))
plt.grid('on'); plt.xlabel('x'); plt.ylabel('y');
plt.legend(('f(x)','f\'(x)'));
a = [ 6 -5 -2 1]
Here’s a simple implementation of Newton’s method with visualization.
def newton(f, df, x0, a, target=1e-6, maxstep=20):
"""Newton's method.
f, df: function and its derivative
x: initial guess
a: parameter for f(x,a)
target: accuracy target"""
x = np.zeros(maxstep+1)
y = np.zeros(maxstep)
x[0] = x0
for i in range(maxstep):
y[i] = f(x[i],a)
if abs(y[i]) < target:
break # converged!
x[i+1] = x[i] - y[i]/df(x[i],a) # new x
else:
print('did not coverge in', maxstep, 'steps.')
return x[:i+1], y[:i+1]
Let us see how it works.
newton(poly, d_poly, 5, a)
(array([5. , 3.88 , 3.27527854, 3.04063338, 3.00110572,
3.00000085, 3. ]),
array([5.60000000e+01, 1.49022720e+01, 3.30409339e+00, 4.17958407e-01,
1.10657543e-02, 8.54778237e-06, 5.11590770e-12]))
newton(poly, d_poly, 0, a)
(array([0. , 1.2 , 0.98978102, 0.99998301, 1. ]),
array([ 6.00000000e+00, -1.15200000e+00, 6.14172290e-02, 1.01951561e-04,
2.88712165e-10]))
Here is a graphical representation.
def zigsawplot(x, y):
"""zigsaw lines of updates
(x0,0),(x0,y0),(x1,0), (x1,y1),...(xn,0),(xn,yn)"""
x = np.repeat(x, 2) # duplicate items
y = np.c_[np.zeros_like(y),y].ravel() # insert zeros
plt.plot(x, y)
xt, yt = newton(poly, d_poly, 5, a)
print('x =', xt[-1], ' y =', yt[-1])
plt.plot(x, poly(x,a))
zigsawplot(xt, yt)
plt.grid('on'); plt.xlabel('x'); plt.ylabel('y');
x = 3.0000000000005116 y = 5.115907697472721e-12
Try starting from other initial guesses
xt, yt = newton(poly, d_poly, -3, a)
print('x =', xt[-1], ' y =', yt[-1])
plt.plot(x, poly(x,a))
zigsawplot(xt, yt)
plt.grid('on'); plt.xlabel('x'); plt.ylabel('y');
x = -2.00000000000004 y = -6.004086117172847e-13
xt, yt = newton(poly, d_poly, 1.9, a)
print('x =', xt[-1], ' y =', yt[-1])
plt.plot(x, poly(x,a))
zigsawplot(xt, yt)
plt.grid('on'); plt.xlabel('x'); plt.ylabel('y');
x = 0.9999999999963238 y = 2.205702287483291e-11
xt, yt = newton(poly, d_poly, 2.1, a)
print('x =', xt[-1], ' y =', yt[-1])
plt.plot(x, poly(x,a))
zigsawplot(xt, yt)
plt.grid('on'); plt.xlabel('x'); plt.ylabel('y');
x = -2.000000000000004 y = -6.039613253960852e-14
Quasi-Newton Method#
The derivative \(f'(x)\) may not be available or hard to compute. In that case, we can use the slope between two points to approximate the derivative.
\[f'(x) \simeq \frac{f(x_k)-f(x_{k-1})}{x_k-x_{k-1}} \]
Quasi-Neuton method uses pairs of points to find the solution:
\[x_{k+1} = x_k - f(x_k)\frac{x_k-x_{k-1}}{f(x_k)-f(x_{k-1})} \]
starting from an initial guess \(x_0\).
def qnewton(f, x0, x1, a, target=1e-6, maxstep=20):
"""Quasi-Newton's method.
f: return f
x: initial guess
x0, x1: initial guess at two points
*args: parameter for f(x,*args)
target: accuracy target"""
x = np.zeros(maxstep+2)
y = np.zeros(maxstep+1)
x[0], x[1] = x0, x1
y[0] = f(x[0], a)
for i in range(1, maxstep+1):
y[i] = f(x[i], a)
if abs(y[i]) < target:
break # converged!
x[i+1] = x[i] - y[i]*(x[i]-x[i-1])/(y[i] - y[i-1]) # new x
else:
print('did not coverge in', maxstep, 'steps.')
return x[:i+1], y[:i+1]
qnewton(poly, 5, 4, a)
(array([5. , 4. , 3.52631579, 3.19955654, 3.05234137,
3.00641035, 3.00022742, 3.00000102, 3. ]),
array([5.60000000e+01, 1.80000000e+01, 7.34800991e+00, 2.28227200e+00,
5.42734403e-01, 6.43913708e-02, 2.27452245e-03, 1.01671083e-05,
1.61831082e-09]))
def qnplot(x, y):
"""lines for quasi-Newton updates
(x0,0),(x0,y0),(x2,0), (x1,0),(x1,y1),(x3,0),...(xn,0),(xn,yn)"""
for i in range(len(x)-2):
plt.plot([x[i],x[i],x[i+2]], [0,y[i],0])
xt, yt = qnewton(poly, 5, 4, a)
print('x =', xt[-1], ' y =', yt[-1])
plt.plot(x, poly(x,a))
qnplot(xt, yt)
plt.grid('on'); plt.xlabel('x'); plt.ylabel('y');
x = 3.000000000161831 y = 1.6183108186851314e-09
xt, yt = qnewton(poly, 0, 2, a)
print('x =', xt[-1], ' y =', yt[-1])
plt.plot(x, poly(x,a))
qnplot(xt, yt)
plt.grid('on'); plt.xlabel('x'); plt.ylabel('y');
x = 0.999999939072668 y = 3.655639955191248e-07
xt, yt = qnewton(poly, -4, -3, a)
print('x =', xt[-1], ' y =', yt[-1])
plt.plot(x, poly(x,a))
qnplot(xt, yt)
plt.grid('on'); plt.xlabel('x'); plt.ylabel('y');
x = -2.0000000573764902 y = -8.60647380918067e-07
Define your own function and see how Newton’s or Quasi-Newton method works.
Discrete-time Dynamics#
We have already seen the iterative dynamics by multiplication of a matrix, which can cause expansion, shrinkage, and rotation.
With nonlinear mapping, more varieties of behaviors including chaos can be observed.
1D Dynamics#
The simplest case is the logistic map.
\[ x_{k+1} = a x_k(1 - x_k) \]
def logistic(x, a):
"""logistic map: f(x) = a*x*(1-x)"""
return a*x*(1 - x)
x = np.linspace(0, 1, 50)
# plot with different levels of a
A = np.arange(5)
leg = []
for a in A:
y = logistic(x, a)
plt.plot(x, y)
leg.append('a = {}'.format(a))
plt.legend(leg)
plt.plot([0,1], [0,1]) # x=f(x) line
plt.xlabel('x'); plt.ylabel('f(x)')
plt.axis('square');
def iterate(f, x0, a, steps=100):
"""x0: initial value
a: parameter to f(x,a)"""
x = np.zeros(steps+1)
x[0] = x0
for k in range(steps):
x[k+1] = f(x[k], a)
return x
Try iteration with different parameter \(a\).
a = 1
xt = iterate(logistic, 0.1, a, 200)
plt.plot(xt)
xt[-10:] # last 10 points
array([0.00490035, 0.00487633, 0.00485256, 0.00482901, 0.00480569,
0.00478259, 0.00475972, 0.00473707, 0.00471463, 0.0046924 ])
a = 3
xt = iterate(logistic, 0.1, a, 200)
plt.plot(xt)
xt[-10:] # last 10 points
array([0.68272042, 0.64983975, 0.68264415, 0.64992335, 0.68256897,
0.65000571, 0.68249486, 0.65008688, 0.68242179, 0.65016688])
a = 3.5
xt = iterate(logistic, 0.1, a, 200)
plt.plot(xt)
xt[-10:] # last 10 points
array([0.38281968, 0.82694071, 0.50088421, 0.87499726, 0.38281968,
0.82694071, 0.50088421, 0.87499726, 0.38281968, 0.82694071])
Trajectory in \(x_k\)-\(x_{k+1}\) plane, called “cobsplot”.
def cobsplot(x):
"""cobsplot of trajectory x"""
plt.plot([0,1], [0,1]) # x=f(x) line
x2 = np.repeat(x, 2) # duplicate items
plt.plot(x2[:-1], x2[1:], lw=0.5) # (x0,x1), (x1,x1), (x1,x2),...
plt.xlabel('$x_k$'); plt.ylabel('$x_{k+1}$');
plt.axis('square');
plt.plot(x, logistic(x, a)) # plot the map
cobsplot(xt) # plot the trajectory
It is known that \(3.57<a<4\) can cause chaotic oscillation.
a = 3.8
xt = iterate(logistic, 0.1, a, 200)
plt.plot(x, logistic(x, a)) # plot the map
cobsplot(xt) # plot the trajectory
2D Dynamics#
Let us see an example of the Hénon map:
\[ x_{k+1} = 1 - a x_k^2 + y_k \]
\[ y_{k+1} = b x_k \]
def henon(xy, ab=[1.4, 0.3]):
"""Henon map: stretch and fold
x_{k+1} = 1 - a*x_k**2 + y_k
y_{k+1} = b*x_k
xy: state [x, y]
ab: parameters [a, b]
"""
x = 1 - ab[0]*xy[0]**2 + xy[1]
y = ab[1]*xy[0]
return np.array([x, y])
x = np.linspace(-1,1,5)
X, Y = np.meshgrid(x,x)
plt.plot(X.T,Y.T); plt.plot(X,Y,':'); # original gird
fXY = henon([X,Y])
plt.plot(fXY[0].T,fXY[1].T); plt.plot(fXY[0],fXY[1],':') # mapped grid
plt.axis('equal');
henon([1,1])
array([0.6, 0.3])
def iterate_vec(f, x0, *args, steps=100):
"""f: n-dimensional map
x0: initial vector
*args: optional parameter to f(x,*args) """
n = len(x0) # state dimension
x = np.zeros((steps+1, n))
x[0] = x0
for k in range(steps):
x[k+1] = f(x[k], *args)
return x
x = iterate_vec(henon, [1, 0])
plt.plot(x);
x = iterate_vec(henon, [1, 0], [1.4, 0.3], steps=5000)
plt.plot(x[:,0], x[:,1], '.', markersize=0.5)
plt.xlabel('x'); plt.ylabel('y'); plt.axis('equal');
Here’s another example: Gumowski-Mira map, which originates from a model of accelerator beams:
\[ x_{n+1} = y_n + \alpha y_n (1-\sigma y_n^2) + g(x_n) \]
\[ y_{n+1} = -x_n + g(x_{n+1}) \]
\[ g(x) = \mu x + \frac{2(1-\mu)x^2}{1+x^2} \]
def gumowski_mira(xy, asm=[0.009, 0.05, -0.801]):
"""Gumowski-Mira map:
x_{k+1} = y_n + alpha y_n (1-sigma y_n^2) + g(x_n)
y_{k+1} = -x_n + g(x_{n+1})
g(x) = mu x + \frac{2(1-mu)x^2}{1+x^2}
xy: state [x, y]
asm: parameters [a, sigma, mu]
"""
x, y = np.array(xy, dtype=float) # unpack array
alpha, sigma, mu = np.array(asm, dtype=float)
# local function
def g(x):
return mu*x + 2*(1-mu)*x**2/(1+x**2)
x1 = y + alpha*y*(1 - sigma*y**2) + g(x)
y1 = -x + g(x1)
return np.array([x1, y1])
x = np.linspace(-20,20,9)
X, Y = np.meshgrid(x,x)
plt.plot(X.T,Y.T); plt.plot(X,Y,':'); # original gird
fXY = gumowski_mira([X,Y])
plt.plot(fXY[0].T,fXY[1].T); plt.plot(fXY[0],fXY[1],':') # mapped grid
plt.axis('equal');
x = iterate_vec(gumowski_mira, [1, 1], steps=10000)
plt.plot(x[:,0], x[:,1], '.', markersize=0.5)
plt.axis('equal');
Recursive call#
In Python and other modern languages, a function and call itself from inside. This can allow compact coding for complex functions.
def factorial(n):
"""factorial by recursice call"""
if n == 0:
return 1
else:
return factorial(n-1)*n
factorial(4)
24
Drawing fractals#
Recursive calls can be used for drawing fractal figures that have multi-scale self-similar structures.
One example is the Cantor set, which is made by iteratively removing the middle part of a segment.
def cantorplot(x0=0, x1=1, a=1/3, e=1e-4):
"""draw a Cantor set
x0, x1: end points
a: fraction to fill
e: minimal resolution"""
u = x1 - x0 # length
if abs(u) < e: # for a very short segment
plt.plot([x0, x1], [0, 0]) # straight line
else:
cantorplot(x0, x0+a*u, a, e) # left 1/3
cantorplot(x1-a*u, x1, a, e) # right 1/3
# Interactive mode to allow zooming in
%matplotlib widget
plt.figure(figsize=(10,2))
cantorplot(a=0.4)
Here is an example of drawing complex lines like a coastline or a clowd by replacing a segment with zig-zag segments.
def fractaline(x0=[0,0], x1=[1,0], a=[0, 0.3, 0], e=1e-2):
"""draw a fractal line
x0, x1: start, end points
a: shifts of via points
e: minimal resolution"""
n = np.size(a) # number of via points
x = np.zeros((n+2,2))
x[0], x[-1] = x0, x1
u = x[-1] - x[0] # connecting vector
v = np.array([-u[1],u[0]]) # orthogonal vector
for i, ai in enumerate(a):
x[i+1] = x0 + u*(i+1)/(n+1) + ai*v # shift of via points
#print(x)
if sum((u**2)) < e**2: # if the segment is very short
plt.plot(x[:,0],x[:,1]) # draw a straight line
else:
for i in range(n+1): # n segments
fractaline(x[i], x[i+1], a, e)
plt.axis('equal');
plt.figure(figsize=(10,4))
fractaline()
plt.figure(figsize=(10,5))
fractaline(a=[0.2])