<jupyter_start><jupyter_text>Introduction to Newton's MethodThis is the **master** version. See the **student** branch for a open ended demo.In the [Ch2_IntroToBisection.ipynb](https://github.com/stephenbeckr/numerical-analysis-class/blob/master/Demos/Ch2_IntroToBisection.ipynb) demo, we looked at [Lambert's W function](https://en.wikipedia.org/wiki/Lambert_W_function). Let's use this example again for root finding, but this time via Newton's method<jupyter_code>import numpy as np

x = 20
f = lambda y : y*np.exp(y) - x


# Let's also get a "true" answer, or at least a fairly good one
import scipy.special.lambertw as W

trueAnswer = np.real( W(20) )
relError = lambda guess : abs( guess - trueAnswer )/abs(trueAnswer)
correctDigits = lambda guess : -np.log10( relError(guess) + 1e-17)

# Plotting

yGrid = np.linspace(0,3,100)

%matplotlib inline
import matplotlib.pyplot as plt
plt.style.use('ggplot')
#plt.style.use('seaborn-notebook')

fig, ax = plt.subplots()
plt.plot(yGrid,f(yGrid))
plt.plot( trueAnswer,f(trueAnswer),'o',markersize=15)
ax.spines['left'].set_position('zero')
ax.spines['right'].set_color('none')
ax.spines['bottom'].set_position('zero')
ax.spines['top'].set_color('none')
ax.minorticks_on()<jupyter_output><empty_output><jupyter_text>Task: implement Newton's Method1. Start with a simple implementation that runs the iterations, maybe printing out the value of x and f(x) at every iteration.  Run the method for a fixed number of iterations (like 20)2. You can make fancier implementations that **save** a list of x and f(x) in case you want to display or plot them later3. You can implement any of the stopping criteria discussed in our book, e.g.,  - $|x_n - x_{n-1}| \le \epsilon$  - $\frac{|x_n - x_{n-1}|}{\max\{|x_,10^{-15}\}} \le \epsilon$  - $|f(x_n)| \le \epsilon$ 4. Try different starting conditions, like  - $x_0 = 2$  - $x_0 = 200$  - $x_0 = -2$  Do all starting conditions lead to convergence? Are they all quadratic convergence right away?<jupyter_code># Basic code
def Newton_simple(f,fprime,x0,maxIter = 100):
  
  x   = np.asarray(x0).copy() # in case it is np.array, we want to do a copy
  # (np.assarray will do nothing if input is already an np.array, and if it isn't,
  #  it'll try to convert it to one, e.g., convert a List to np.array)
  
  for n in range(1,maxIter):
    x  = x - f(x)/fprime(x)
  return x

# Fancier code (with saved output and nicer stopping conditions)
def Newton(f,fprime,x0,maxIter = 100, fTol = 1e-8, relTol = 1e-8,Verbose=False):
  history_x  = np.zeros(maxIter)
  history_fx = np.zeros(maxIter)
  x   = np.asarray(x0,dtype=np.double).copy()
  fx  = f(x)
  history_x[0]  = x
  history_fx[0] = fx
  for n in range(1,maxIter):
    try:
      x -= fx / fprime(x)
    except ZeroDivisionError:
      return x, history_x, history_fx
    #print(x,fprime(x)) # for debugging
    if Verbose:
      print("Iteration {:4d}, x is {:+14.8e}, f(x) is {:+14.8e},  f'(x) is {:+14.8e}".format(n,x,fx,fprime(x)))
    fx = f(x)
    history_x[n] = x
    history_fx[n] = fx
    if np.abs(fx) <= fTol or np.abs(x-history_x[n-1])/max(1e-15,x) <= relTol:
      print('Reached stopping criteria at iteration',n)
      history_x  = history_x[:n]
      history_fx = history_fx[:n]
      break
  return x, history_x, history_fx
# Now, run Newton's method on the function f
fprime = lambda y : np.exp(y) + y*np.exp(y)

x0  = 20
p, history_x, history_fx = Newton(f,fprime,x0,maxIter=50,fTol = 1e-15,relTol = 1e-15)
p
plt.plot(abs(history_fx),'o-',MarkerSize=12);
plt.yscale("log")
plt.xlabel("Iterations n")
plt.ylabel("f(x_n)");<jupyter_output>Reached stopping criteria at iteration 26<jupyter_text>Discussion of results:- for $x_0=2$ we get very rapid convergence right away- for $x_0=20$, we're not close enough for the theorems to apply. We still convergen, just at a linear rate, not quadratic. Then, once we get close enough, the theorems do apply, and we get very fast quadratic convergence.  If we keep going for further iterations, we don't get any lower error because we've reached the limit of floating point precisino- for $x_0 = -2$ we diverge to $-\infty$  Task: implement the Secant MethodYou'll notice one important detail: we need to specify **two** initial points, $x_0$ and $x_1$<jupyter_code>def Secant(f,x0,x1,maxIter = 100, fTol = 1e-8, relTol = 1e-8):
  history_x  = np.zeros(maxIter)
  history_fx = np.zeros(maxIter)

  
  x_old= np.asarray(x0).copy() # .copy() probably not necessary
  x    = np.asarray(x1).copy()
  fx_old = f(x_old)
  fx     = f(x)
  history_x[0]  = x
  history_fx[0] = fx

  for n in range(1,maxIter):
    xNew  = x - fx*(x-x_old)/(fx-fx_old)
    x_old = x
    fx_old= fx
    x     = xNew

    fx = f(x)
    # print(x,x_old,fx,fx_old) # debugging
    history_x[n] = x
    history_fx[n] = fx
    if np.abs(fx) <= fTol or np.abs(x-history_x[n-1])/max(1e-15,x) <= relTol:
      print('Reached stopping criteria at iteration',n)
      history_x  = history_x[:n]
      history_fx = history_fx[:n]
      break
  return x, history_x, history_fx
# Now, run the Secant method on the function f
x0  = 2
x1  = x0 + .1
p, history_x, history_fx = Secant(f,x0,x1,maxIter=50,fTol = 1e-15,relTol = 1e-15)
p, Newton_history_x, Newton_history_fx = Newton(f,fprime,x1,maxIter=50,fTol = 1e-15,relTol = 1e-15)

plt.plot(abs(history_fx),'o-',MarkerSize=12,label="Secant method");
plt.plot(abs(Newton_history_fx),'*-',MarkerSize=12,label="Newton's method");
plt.yscale("log")
plt.xlabel("Iterations n");
plt.ylabel("f(x_n)");
plt.legend();<jupyter_output>Reached stopping criteria at iteration 7
Reached stopping criteria at iteration 5<jupyter_text>Task: compare Newton's method, the Secant method, and bisection method<jupyter_code><jupyter_output><empty_output><jupyter_text>Task: "break" Newton's methodFind a function $f(x)$ that satisfies all of the following properties:* $f$ is *strictly monotonically increasing** $f$ has a *simple root* at $x=0$* $f'(x)$ exists for all real $x$* Newton's method diverges when started with an initial guess $x_0 = 1$. Definition* A function $f$ is *strictly monotonically increasing* if $f(a)  0$ for all $x$. One answer: $f(x)=\arctan(c\cdot x)$There are lots of possiblities.  It helps to draw a sketch of a monotonically increasing function with a root at $x=0$, and then start drawing tangent lines to see what can go wrong.  One choice is $f(x)=\arctan(c\cdot x)$ for a large enough constant $c$ (which has the effect of using $\arctan(x)$ but for $x_0 = c$); in fact, Newton's method is "affine invariant" which means scaling the inside by $x \mapsto c\cdot x + b$ will not affect the iterations, though it *does* affect the starting guess.<jupyter_code>c      = 5
f      = lambda x : np.arctan(c*x)
fprime = lambda x : c/(1+(c*x)**2)

x0    = .1
print("Running Newton's method starting at x0 = {}\n".format(x0))
p, Newton_history_x, Newton_history_fx = Newton(f,fprime,x0,maxIter=20,fTol = 1e-15,relTol = 1e-15,Verbose=True)

# Now, try it again for a starting point farther away
x0    = 1
print("\nRunning Newton's method starting at x0 = {}\n".format(x0))
p, Newton_history_x, Newton_history_fx = Newton(f,fprime,x0,maxIter=12,fTol = 1e-15,relTol = 1e-15,Verbose=True)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13666574356134034, 'fraction_numerical': 0.024646912212683468, 'mean_word_length': 4.30712711241734, 'pattern_counts': {'":': 0, '<': 23, '<?xml version=': 0, '>': 19, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '398694', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 2831, 'n_tokens_neox': 2677, 'n_words': 795, 'path': 'Demos/Ch2_NewtonsMethod.ipynb', 'repo_name': 'jakemcgrath/numerical-analysis-class', 'size': '7222'}