Disclaimer: This code is far from beeing efficient or optimal. It's intention is be simple, while using numpy and matplotlib for convenience.
In [ ]:
%matplotlib inline
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from IPython.display import set_matplotlib_formats
set_matplotlib_formats('png')
mpl.rcParams['savefig.dpi'] = 140
logistic map: $$x_{n+1} = rx_n(1-x_n)$$
In [ ]:
def plot_feigenbaum(r_intervall = [2.4, 4.0],
n_steps = 800,
n_pre_iterations = 200,
n_show_iterations = 200):
r = np.linspace(*r_intervall, num=n_steps)
x = np.random.random(n_steps)
x_results = np.zeros([n_steps, n_show_iterations])
for i in range(n_pre_iterations + n_show_iterations):
x = r * x * (1 - x)
if i >= n_pre_iterations:
x_results[:, i - n_pre_iterations] = x
plt.plot(r[np.newaxis,:].repeat(n_show_iterations, axis=1).flatten(),
x_results.flatten(),
marker='.', linestyle='None', color='black', markersize=0.2)
plt.title('$x$ values over {} iterations (after first {})'.
format(n_show_iterations, n_pre_iterations))
plt.xlabel('$r$')
plt.ylabel('$x$')
plot_feigenbaum()
In [ ]:
plot_feigenbaum(r_intervall = [3.6, 3.8], n_steps=1000, n_pre_iterations = 1000)
In [ ]:
plot_feigenbaum(r_intervall = [3.74, 3.745], n_steps=1000, n_pre_iterations = 1000, n_show_iterations=500)
plt.ylim((0.46, 0.53))
In [ ]:
def plot_barnsley(n_points=10**5):
# matrices for the linear part of the transformations
A = np.array([
[[.0, .0],
[.0, .16]],
[[.85, .04],
[-.04, .85]],
[[.2, -.26],
[.23, .22]],
[[-.15, .28],
[.26, .24]]
])
# translation vectors of the transformations
b = np.array([
[.0, .0],
[.0, 1.6],
[.0, 1.6],
[.0, .44]
])
# probabilities with which each transformation is chosen
p = [.01, .85, .07, .07]
p_aggr = [p[0], p[0]+p[1], p[0]+p[1]+p[2]]
points = np.zeros((n_points, 2)) # stores the result of each iteration
x = np.array([.0, .0])
for i in range(n_points):
# pick a transformation
s = np.random.random()
if s < p_aggr[0]:
i_trafo = 0
elif s < p_aggr[1]:
i_trafo = 1
elif s < p_aggr[2]:
i_trafo = 2
else:
i_trafo = 3
x = np.dot(A[i_trafo], x) + b[i_trafo]
points[i] = x
plt.figure(figsize=(6,6))
plt.plot(points[:,0], points[:,1],
marker='.', linestyle='None', color='green', markersize=0.5)
plt.title('Barnsley fern')
plt.xlabel('$x_1$')
plt.ylabel('$x_2$')
plot_barnsley()
In [ ]:
def plot_julia(c = -0.7 + 0.3j,
z_min = -2 - 2j,
z_max = 2 + 2j,
max_iter = 160,
breakout = 100.0,
resolution = [800, 600]):
iteration = np.zeros(resolution, dtype='int')
z = (np.linspace(z_min.real, z_max.real, resolution[0])[:, np.newaxis] +
1j * np.linspace(z_min.imag, z_max.imag, resolution[1])[np.newaxis, :])
for i in range(max_iter):
z = z**2 + c
iteration[(np.abs(z) > breakout) & ~iteration.astype('bool')] = i
plt.figure()
plt.imshow(iteration.T,
extent=[z_min.real, z_max.real, z_min.imag, z_max.imag],
interpolation='nearest', origin='lower')
plt.xlabel('$\Re(z)$')
plt.ylabel('$\Im(z)$')
plt.title('Julia set for $c={}$'.format(c))
plt.colorbar(label='breakout iteration')
plot_julia()
In [ ]:
plot_julia(z_min = -1.0 - 0.5j, z_max= 0.0 + 0.5j, max_iter = 260)
In [ ]:
plot_julia(c = 0.41 + 0.2j, max_iter = 80)
In [ ]:
def plot_mandelbrot(c_min = -2.4 - 1.5j,
c_max = 1.2 + 1.5j,
max_iter = 30,
breakout = 100.0,
resolution = [800, 600]):
iteration = np.zeros(resolution, dtype='int')
z = 0.0
c = (np.linspace(c_min.real, c_max.real, resolution[0])[:, np.newaxis] +
1j * np.linspace(c_min.imag, c_max.imag, resolution[1])[np.newaxis, :])
for i in range(max_iter):
z = z**2 + c
iteration[(np.abs(z) > breakout) & ~iteration.astype('bool')] = i
plt.figure()
plt.imshow(iteration.T,
extent=[c_min.real, c_max.real, c_min.imag, c_max.imag],
interpolation='nearest', origin='lower')
plt.xlabel('$\Re(c)$')
plt.ylabel('$\Im(c)$')
plt.title('Mandelbrot set')
plt.colorbar(label='breakout iteration')
plot_mandelbrot()
In [ ]:
plot_mandelbrot(c_min=-0.4+0.6j, c_max=0.2+1.1j, max_iter=80)
In [ ]:
def plot_newton(z_min = -1.2 - 1j,
z_max = 1.2 + 1j,
p = [3, 0, 0, -1], # cofficients of the polynomial for numpy.polyval, e.g. p[0] * z*1 + p[1] * z*0
eps = 0.001,
a = 1,
max_iter = 30,
resolution = [800, 600]):
f = np.poly1d(p)
df = f.deriv()
roots = f.r
print('roots:')
for root in roots:
print(root)
iteration = np.zeros(resolution, dtype='int') # iteration in which the distance to a root is less than eps
converged = np.zeros(resolution, dtype='int') # index of the root to which each point converges
z = (np.linspace(z_min.real, z_max.real, resolution[0])[:, np.newaxis] +
1j * np.linspace(z_min.imag, z_max.imag, resolution[1])[np.newaxis, :])
for i in range(max_iter):
z = z - a * f(z) / df(z)
for i_root, root in enumerate(roots):
mask = np.abs(z - root) < eps
iteration[mask & ~converged.astype('bool')] = i
converged[mask] = i_root + 1
plt.figure()
plt.imshow(iteration.T,
extent=[z_min.real, z_max.real, z_min.imag, z_max.imag],
interpolation='nearest', cmap=plt.cm.gray, origin='lower')
plt.colorbar(label='convergence iteration')
plt.imshow(converged.T,
extent=[z_min.real, z_max.real, z_min.imag, z_max.imag],
interpolation='nearest', cmap=plt.cm.hsv, alpha=0.4, origin='lower')
plt.xlabel('$\Re(z)$')
plt.ylabel('$\Im(z)$')
plt.title('Newton basin for ${}$'.format(p))
plot_newton(resolution=[1600, 1200])
In [ ]:
plot_newton(p=[1, 0, 0, 1, 0, 0, -1])
In [ ]:
plot_newton(p=[3, 0, 0, 2, 0, 0, 1])
In [ ]: