#!/usr/bin/env python # coding: utf-8 # # Root Finding # # ### Table of Contents: # # 1. [Fixed-point iteration](#first-bullet) # 2. [Newton's method](#second-bullet) # 3. [Secant method](#third-bullet) # 4. [Characteristics](#fourth-bullet) # 5. [Implementation and comparison](#fifth-bullet) # # --- # ## 1. Fixed-point iteration # # Suppose we define an iteration # # $$ # x_{k+1} = g(x_k) # $$ # # If $\alpha$ is a point such that $g(\alpha) = \alpha$, we call $\alpha$ a fixed point of $g$. # # A fixed-point iteration terminates once a fixed point is reached, since if $g(x_k)= x_k$ then we get $x_{k+1} = x_k$. # # Also, if $x_{k+1}$ converges as $k \rightarrow \infty$, it must converge to a fixed point: Let $\alpha = \lim_{k \rightarrow \infty} x_k$, then # # $$ # \alpha = \lim_{k \rightarrow \infty} x_{k+1} = \lim_{k \rightarrow \infty} g(x_k) = g\left(\lim_{k \rightarrow \infty} x_k\right) = g(\alpha) # $$ # # Hence, we need to guarantee that a fixed-point iteration will converge. # # --- # We can use **Lipschitz condition** to prove that the fixed point iteration converges to $\alpha$. # # If $g $ satisfies a Lipschitz condition in an interval $[a,b]$ if $\exists L \in \mathbb{R}_{>0}$ such that # # $$ # |g(x) - g(y)| \leq L |x-y|, \quad \forall x,y \in [a,b] # $$ # # $g$ is called a **contraction** if $L < 1$. Simply speaking, contraction takes two numbers and maps them close to each other. # # $$ # |x_k - \alpha| = |g(x_{k-1}) - g(\alpha)| \leq L|x_{k-1} - \alpha| # $$ # # which implies # # $$ # |x_k - \alpha| \leq L^k |x_0 - \alpha| # $$ # # and since $L < 1$, $|x_k - \alpha| \rightarrow 0$ as $k \rightarrow \infty$. This shows that error decreases by factor of $L$ each iteration. # # Therefore, we define the rule of convergence of fixed-point iteration: # # > Suppose that $g(\alpha) = \alpha$ and that $g$ is a contraction on $[\alpha - A, \alpha + A]$. Suppose also that $|x_0 - \alpha| \leq A$. Then the fixed point iteration converges to $\alpha$. # --- # We can verify convergence of a fixed-point iteration by checking the gradient of $g$, since if $g \in C^1[a,b]$, we can obtain a Lipschitz constant based on $g'$: # # $$ # L = \max_{\theta \in (a,b)} |g'(\theta)| # $$ # # We want to check if $|g'(\alpha)|<1$, i.e. if there is a neighborhood of $\alpha$ on which $g$ is a contraction. # # Furthermore, as $k \rightarrow \infty$, # # $$ # {|x_{k+1} - \alpha| \over |x_k - \alpha|} = {| g(x_k) - g(\alpha)| \over |x_k - \alpha|} \rightarrow |g'(\alpha)| # $$ # # Hence, asymptotically, error decreases by a factor of $|g'(\alpha)|$ each iteration. # # $|g'(\alpha)|$ can be used to dictate the **asymptotic convergence rate**: # # We say an iteration converges **linearly** if, for some $\mu \in (0,1)$, # # $$ # \lim_{k \rightarrow \infty} {|x_{k+1} - \alpha| \over |x_k - \alpha|} = \mu # $$ # # An iteration converges **superlinearly** if # # $$ # \lim_{k \rightarrow \infty} {|x_{k+1} - \alpha| \over |x_k - \alpha|} = 0 # $$ # --- # ## 2. Newton's method # # Consider the following fixed-point iteration # # $$ # x_{k+1} = x_k - \lambda(x_k)f(x_k), \quad k = 0,1,2,... # $$ # # corresponding to $g(x) = x - \lambda(x_k) f(x_k)$, for some function $\lambda$. # # We want to achieve a fixed point $\alpha$ of $g$ such that $f(\alpha) = 0$, and we also want $|g'(\alpha)| = 0$ t get superlinear convergence. # # We take the first derivative and evaluate it at $f(\alpha) = 0$: # # $$ # g'(\alpha) = 1 - \lambda'(\alpha) f(\alpha) - \lambda(\alpha) f'(\alpha) = 1 - \lambda(\alpha) f'(\alpha) # $$ # # To satisfy $|g'(\alpha)| = 0$ we choose $\lambda(x) = 1/f'(x)$ to get **Newton's Method**: # # $$ # x_{k+1} = x_k - {f(x_k)\over f'(x_k)}, \quad k = 0,1,2,... # $$ # --- # To understand the superlinear convergence rate of Newton's Method, we apply **Taylor expansion** for $f(\alpha)$ about $f(x_k)$: # # $$ # 0 = f(\alpha) = f(x_k) + (\alpha - x_k) f'(x_k) + {(\alpha - x_k)^2 \over 2} f''(\theta_k) # $$ # # for some $\theta_k \in (\alpha , x_k)$. # # Dividing through by $f'(x_k)$ gives # # $$ # \left(x_{k+1} - {f(x_k)\over f'(x_k)}\right) - \alpha = {f''(\theta_k)\over 2f'(x_k)} (x_k - \alpha)^2 # $$ # # equivalently, # # $$ # x_{k+1} - \alpha = {f''(\theta_k)\over 2f'(x_k)} (x_k - \alpha)^2 # $$ # # Hence, the error at iteration $k+1$ is the square of the error at each iteration $k$. This refers to as **quadratic convergence**. # # Note that this result relies on the second order Taylor expansion near $\alpha$, we need to be **sufficiently close** to $\alpha$ to get quadratic convergence. # --- # ## 3. Secant method # # # Secant method uses finite difference to approximate $f'(x_k)$: # # $$ # f'(x_k) \approx {f(x_k) - f(x_{k-1}) \over x_k - x_{k-1}} # $$ # # Substituting this into the iteration leads to **Secant Method**: # # $$ # x_{k+1} = x_k - f(x_k)\left({x_k - x_{k-1}\over f(x_k) - f(x_{k-1})}\right), \quad k = 1,2,3... # $$ # # # The convergence rate fo Secant Method can be shown as # # $$ # \lim_{k \rightarrow \infty} {|x_{k+1} - \alpha| \over |x_k - \alpha|^q} = \mu # $$ # # where $\mu$ is a positive constant and $q$ is the golden ratio $q = {1+ \sqrt{5} \over 2} \approx 1.6$. # # Thus, Secant Method has a **golden ratio quadratic rate** of convergence. # --- # ## 4. Characteristics # # * **Fixed-point iteration** # # * Sometimes need to re-write the function $f(x)=0$ as $x = g(x)$ to make sure $g(x)$ has an appropriate property. # # * **Newton's Method** # # * Quadratic convergence is very rapic. # * Requires to evaluate both $f(x_k)$ and $f'(x_k)$ per iteration. # # * **Secant Method** # # * Faster than fixted-point iteration, slower than Newton's Method. # * Does not require us to determine $f'(x)$ analytically. # * Requires only one extra function evaluation, $f(x_k)$, per iteration. # * When iterations become so close, we might get division by 0. # # --- # ## 5. Implementation and comparison # # The following section compares how fixed-point interation, Newton's method and Secant method solve $f(x) = e^x - x -2 = 0$. # In[1]: from math import * # f(x) = e^x-x-2 def f(x): return exp(x)-x-2 # f'(x) = e^x-1 def df(x): return exp(x)-1 startpt = 2 # Newton method setup xa=startpt # Secant method setup xb=startpt # Fixed-point iteration setup xc=startpt # Initialize previous step x_{k-1} in secant method xbb=startpt + 0.1 fxbb=f(xbb) # store xas = [xa] xbs = [xb] xcs = [xc] for k in range(30): print("%17.10g %17.10g %17.10g %17.10g %17.10g %17.10g" \ %(xa,f(xa),xb,f(xb),xc,f(xc))) # Newton xa=xa-f(xa)/df(xa) xas.append(xa) # Secant fxb=f(xb) tem=xb-fxb*(xb-xbb)/(fxb-fxbb) xbb=xb;fxbb=fxb xb=tem xbs.append(xb) # Fixed point iteration xc=log(xc+2) xcs.append(xc) if xb-xbb == 0.0 and fxb-fxbb == 0.0: break # In[2]: import matplotlib.pyplot as plt import numpy as np domain = np.arange(0,len(xas), 1) plt.figure(figsize=(15, 8), dpi=80) plt.plot(domain,[f(x) for x in xas], label = "Newton's Method") plt.plot(domain,[f(x) for x in xbs], label = "Secant Method") plt.plot(domain,[f(x) for x in xcs], label = "Fixed-point Iteration") plt.ylabel('f(x)') plt.xlabel('x') plt.legend() plt.show()