In numerical analysis, Newton's method, also known as the Newton–Raphson method is a root-finding algorithm which produces successively better approximations to the roots (or zeroes) of a real-valued function. The most basic version starts with a single-variable function f defined for a real variable x, the function's derivative f′, and an initial guess x0 for a root of f. If the function satisfies sufficient assumptions and the initial guess is close, then
$$ x_{1}= x_{0}-\frac {f(x_{0})}{f'(x_{0})}$$is a better approximation of the root than x0. Geometrically, (x1, 0) is the intersection of the x-axis and the tangent of the graph of f at (x0, f(x0)): that is, the improved guess is the unique root of the linear approximation at the initial point. The process is repeated as
$$x_{n+1}=x_{n}-\frac {f(x_{n})}{f'(x_{n})}$$using numpy below is another implementation that I found on GitHub.
The %%timeit command lets you test how fast code runs: this is called benchmarking and it's very easy to do in Jupyter.
#%%timeit
# zeroes of z -> z^3 - 1
npts = 512
X = np.linspace(-2, 2, npts)
Y = np.linspace(-2, 2, npts)
#broadcast X,Y to a square array
C = X[:, None] + 1J*Y
#Newton for z -> z^3 - 1
max_iter = 20
for k in range(max_iter):
C = C - (C**3 - 1)/ (3 * C**2)
#make a blank rgb image
im = np.zeros((C.shape[0],
C.shape[1],
3))
# set colors according to the basin of attraction
im[C.imag > 0,0] = 1
im[C.imag < 0,2] = 1
yellow = C.real > 0
im[yellow,1] = 1
im[yellow,0] = 1
fig = plt.figure(figsize=(8,8))
ax = fig.add_subplot(111)
ax.imshow(im);
#plt.savefig("newton.png")
#%%timeit
# zeroes of z -> (z^2 - 1)(z^2 -t^2)
def mk_grid_cplx(npts=2**9):
X = np.linspace(-2, 2, npts)
Y = np.linspace(-2, 2, npts)
#broadcast X,Y to a square array
return X[:, None] + 1J*Y
def do_frame(t=1):
C = mk_grid_cplx(npts=2**9)
max_iter = 40
for k in range(max_iter):
U = C**2
C = C - (U - 1)*(U - t**2)/ ( 2*C*(2*U - (1 + t**2)) )
#make a blank rgb image
im = np.zeros((C.shape[0],
C.shape[1],
3),
dtype=np.uint8)
# set colors according to the basin of attraction
im[np.abs(C - 1) < .1 , 0] = 255 #R
im[np.abs(C + 1) < .1, 2] = 255 #B
yellow = np.abs(C - t) < .1
im[yellow, 1] = 255
im[yellow, 0] = 255
return im
#frames = [do_frame(t=s) for s in np.linspace(4,.1,40)]
frames = [do_frame(t= .9*np.exp(s*1J))
for s in np.linspace(0, 2*np.pi, 160)]
from matplotlib import animation
from IPython.display import HTML
%matplotlib inline
#frames = imgs
def plot_images(img_list):
def init():
img.set_data(img_list[0])
return (img,)
def animate(i):
img.set_data(img_list[i])
return (img,)
fig = plt.figure()
ax = fig.gca()
plt.axis('off')
img = ax.imshow(img_list[0])
anim = animation.FuncAnimation(fig, animate,
init_func=init,
frames=len(img_list),
interval=200,
blit=True)
return anim
anim = plot_images(frames).to_html5_video()
HTML(anim)
with open("newton_frac.html", "w") as fp:
fp.write(anim)
HTML(anim)
if you have python 3.9 this isn't easy
pip3 install opencv-python
import cv2
size = frames[0].shape
fps = 5
out = cv2.VideoWriter('xx.mp4',
cv2.VideoWriter_fourcc( * 'mp4v'),
fps,
(size[0], size[1]))
for i, ff in enumerate(frames):
out.write(ff)
out.release()
sucks this doesn't work
should use but doesn't really work
array([[0., 0., 1.], [0., 0., 1.], [0., 0., 1.], ..., [0., 0., 1.], [0., 0., 1.], [0., 0., 1.]])
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap
# A list of colors to distinguish the roots.
colors = ['b', 'r', 'g', 'y']
TOL = 1.e-8
def newton(z0, f, fprime, MAX_IT=1000):
"""The Newton-Raphson method applied to f(z).
Returns the root found, starting with an initial guess, z0, or False
if no convergence to tolerance TOL was reached within MAX_IT iterations.
"""
z = z0
for i in range(MAX_IT):
dz = f(z)/fprime(z)
if abs(dz) < TOL:
return z
z -= dz
return False
def plot_newton_fractal(f, fprime, n=200, domain=(-1, 1, -1, 1)):
"""Plot a Newton Fractal by finding the roots of f(z).
The domain used for the fractal image is the region of the complex plane
(xmin, xmax, ymin, ymax) where z = x + iy, discretized into n values along
each axis.
"""
roots = []
m = np.zeros((n, n))
def get_root_index(roots, r):
"""Get the index of r in the list roots.
If r is not in roots, append it to the list.
"""
try:
return np.where(np.isclose(roots, r, atol=TOL))[0][0]
except IndexError:
roots.append(r)
return len(roots) - 1
xmin, xmax, ymin, ymax = domain
for ix, x in enumerate(np.linspace(xmin, xmax, n)):
for iy, y in enumerate(np.linspace(ymin, ymax, n)):
z0 = x + y*1j
r = newton(z0, f, fprime)
if r is not False:
ir = get_root_index(roots, r)
m[iy, ix] = ir
nroots = len(roots)
if nroots > len(colors):
# Use a "continuous" colormap if there are too many roots.
cmap = 'hsv'
else:
# Use a list of colors for the colormap: one for each root.
cmap = ListedColormap(colors[:nroots])
plt.imshow(m, cmap=cmap, origin='lower')
plt.axis('off')
plt.show()
f = lambda z: z**4 - 1
fprime = lambda z: 4*z**3
plot_newton_fractal(f, fprime, n=500)
(512, 512)