""" Sam Dolan. Registration number: xxxxxxxxx """ ##### # (1) Import modules import numpy as np import matplotlib.pyplot as plt import time ##### # (2) Define all the functions we will use here # (Add comments and a docstring). def nr_step_f1(x): """This function takes one step of the NR method to find the root of the function f_1(x) = x e^{-x} - 1/10.""" f1 = x*np.exp(-x) - 0.1 f1x = (1-x)*np.exp(-x) return x - f1/f1x def nr_step_f2(x): f2 = x + np.sin(x) - np.log(x) - 2.0 f2x = 1 + np.cos(x) - 1.0/x return x - f2/f2x def nr_step_f3(z): f3 = z**3 - 1.0 f3z = 3*z**2 return z - f3/f3z def findroot(x0, step_func, max_steps=25, tol=1e-12): """ Seeks a root of f1(x) by applying the Newton-Raphson step repeatedly to generate a sequence. Here "x0" is the initial guess, and "step_func" is a function (such as nr_step_f1) that carries out a single iterative step. The method should stop once either of two conditions is met: (i) the absolute difference between successive estimates is less than "tol" (indicating that the sequence has converged sufficiently), or, (ii) the number of steps taken exceeds "max_steps". The function returns: A list of the sequence of estimates. """ xs = [x0] for i in range(max_steps): #print("%d:\t%.14f" % (i, xs[-1])) xs.append(step_func(xs[-1])) if abs(xs[-1]-xs[-2]) < tol: break #print("------------------------") #print("Result:\t%.12f" % xs[-1]) return xs print("Part 1(a)") print("See code: I have added comments and a docstring to the function.") print("-----") print("Part 1(b)") x0 = 0.2 seq = findroot(x0, nr_step_f1) print("With a starting value of x0=%.2f the NR method converged after %i iterations to:" % (x0, len(seq)-1)) print(" %.10f" % seq[-1]) x0 = 2.0 seq = findroot(x0, nr_step_f1) print("With a starting value of x0=%.2f the NR method converged after %i iterations to:" % (x0, len(seq)-1)) print(" %.10f" % seq[-1]) print("-----") print("Part 1(c)") print("The function f2(x) has four roots in the domain (0,5]:") print(" %.10f" % findroot(0.001, nr_step_f2)[-1]) print(" %.10f" % findroot(1.0, nr_step_f2)[-1]) print(" %.10f" % findroot(3.0, nr_step_f2)[-1]) print(" %.10f" % findroot(4.0, nr_step_f2)[-1]) print("-----") print("Part 1(d)") print(" This part is for discussion, and is not assessed.") print("-----") print("Part 2(a)") print("The NR method fails for a starting guess of 0, and for the three values that ") print("map onto zero, and the 9 values that map onto those 3, and so on.") print("Solving z_{1} = 0 gives 2 z^3 + 1 = 0. The three roots to this equation can be") print("written in closed form, or found numerically using the NR method again:") def nr_step_f4(z): f4 = 2*z**3 + 1.0 f4z = 6*z**2 return z - f4/f4z ics = [1.0, 1.0+1j, 1.0-1j] zbad = [findroot(ic, nr_step_f4)[-1] for ic in ics] def printcomplex(z): print(" (%.6f, %.6f)" % (z.real, z.imag)) # without carriage return for z in zbad: printcomplex(z) print("\n-----") print("Part 2(b)") print("Testing the NR method for finding the cubic roots of unity.") zroots = [1.0, 0.5*(-1.0+np.sqrt(3.0)*1j), 0.5*(-1.0-np.sqrt(3.0)*1j)] z0s = [0.02+0*1j, 0.02*1j, -0.01*(1+1j), -2**(-1/3.0)+0.001*1j] for z0 in z0s: seq = findroot(z0, nr_step_f3, max_steps=50) print("\nWith a starting value of z0=(%.4f, %.4f) the NR method converged after %i iterations to:" % (z0.real, z0.imag, len(seq)-1)) print(" %.6f + %.6f i" % (seq[-1].real, seq[-1].imag)) print("-----") print("Part 2(c)") print("Creating a 2D image of the basins of attraction of the 3 roots in the complex plane.") N1 = 200 N2 = 200 xwid = 1.4 xs = np.linspace(-xwid,xwid,N1) ys = np.linspace(-xwid*N2/N1,xwid*N2/N1,N2) arr = np.zeros((N1,N2)) def iterate(z): z2 = z**2 z3 = z2*z return z - (z3 - 1.0)/(3.0*z2) eps = 1.e-12 maxits = 25 def intval(z): return (abs(z-zroots[0]) < eps)*1 + (abs(z-zroots[1]) < eps)*2 + (abs(z-zroots[2]) < eps)*3 print("Calculating ...") t1 = time.time() for i in range(N1): for j in range(N2): z = xs[j] + 1j*ys[i] for k in range(maxits): z = iterate(z) arr[i,j] = intval(z) t2 = time.time() print("Time taken for calculation : %f sec" % (t2-t1)) print("Plotting ...") plt.imshow(arr) # plt.show() plt.savefig("nrbasins.png")