Newton's Method

Algorithm

In [ ]:
'''
To find a solution to f(x) = 0 given an initial approximation pâ‚€ :
INPUT initial approximation pâ‚€ ; tolerance TOL; maximum number of iterations Nâ‚€ .
OUTPUT approximate solution p or message of failure.
Step 1 Set i = 1.
Step 2 While i ≤ N₀ do Steps 3–6.
    Step 3 Set p = p₀ − f(p₀ )/f'(p₀). (Compute pᵢ .)
    Step 4 If |p − p₀ | < TOL then
            OUTPUT (p); (The procedure was successful.)
            STOP.
    Step 5 Set i = i + 1.
    Step 6 Set pâ‚€ = p. (Update pâ‚€ .)
Step 7 OUTPUT (‘The method failed after N₀ iterations, N₀ =’, N₀ );
        (The procedure was unsuccessful.)
        STOP.
'''

Code

In [1]:
import pandas as pd
from sympy import *
from plotly.offline import download_plotlyjs, init_notebook_mode, iplot
from plotly.graph_objs import *
import numpy as np

init_notebook_mode(connected=True)

def NewtonsMethod(
    p0, 
    equation, 
    tolerance=0.000001, 
    max_iterations=50, 
    output='answer', 
    plot=False,
    stopping_criterion='tolerance', 
    iteration=50 
):
    '''
    The Newton's Method
    
    Keyword arguments:
    p0 (double) -- initial approximation pâ‚€
    equation -- an equation of 'x' that we need to evaluate
    tolerance (double) -- error tolerance (default = 0.000001)
    max_iterations (int) -- number of maximum iterations (default = 50)
    output (string) -- how to output the result (default = 'answer')
                possible values are:
                  - 'answer' : to show the value of "p"
                  - 'table' : to show all the values in tabular form
    plot (bool) -- whether to show plot or not (default=False)
    stopping_criterion (string) -- when to stop the processing (default = 'tolerance')
                possible vaulse are:
                  - 'tolerance' : to stop at a certain tolerance level
                  - 'iteration' : to stop at a certain iteration (must provide value for 'iteration' argument)
    iteration (int) -- on which iteration to stop the processing (default = 50)
    
    Returns:
    
    It will return the value according to the value in 'output' argument.
    If value of output is 'answer', it will return a float. If it is 'table', it will return a DataFrame.
    '''
    
    f = lambdify(x, equation)
    df = lambdify(x, equation.diff(x))
    
    N = 100
    X = np.linspace(-5, 5, N)
    Y = list(map(f, X))
    
    if (output == "table" or plot):
            P = list()
            FP = list()
            labels = list()
            P.append(p0)
            FP.append(f(p0))
            labels.append('P0')
            
    for i in range(1, max_iterations):
        p = p0 - (f(p0)/df(p0))
        
        if (output == "table" or plot):
            P.append(p)
            FP.append(f(p))
            labels.append("P"+str(i))
        
        if (abs(p - p0) < tolerance or (stopping_criterion == 'iteration' and iteration == i)):
            if (plot):
                trace0 = Scatter(
                x=X,
                y=Y,
                mode = 'lines',
                name = "equation"
                )

                trace1 = Scatter(
                    x=P,
                    y=FP,
                    mode = 'markers',
                    name = "(p, f(p))",
                    text = labels
                )

                trace2 = Scatter(
                    x=P,
                    y=np.zeros((len(P),), dtype=np.int),
                    mode = 'markers',
                    text = labels,
                    name = "p"
                )
                layout = Layout(
                    xaxis=dict(
                        range=[-5, 5]
                    ),
                    yaxis=dict(
                       range=[-20, 20]
                    )
                )
                data = Data([trace0, trace1, trace2])
                fig = Figure(data=data, layout=layout)

                iplot(fig)
            
            if (output == 'answer'):
                return p
            if (output == 'table'):
                # create a table and return it
                d = {'p' : P}
                return pd.DataFrame(data=d).style.format("{:.10}")
        
        p0 = p
    
    return "Method failed after " + max_iterations + " iterations."

Example 1

Consider the function f(x) = cosx−x = 0. Approximate a root of f using Newton’s Method

In [3]:
import math
from sympy import *

x = Symbol('x')
NewtonsMethod(math.pi/4, cos(x) - x, output="table", plot=True)
Out[3]:
p
0 0.7853981634
1 0.7395361335
2 0.7390851781
3 0.7390851332

Exercise Set 2.3

1. Let f(x) = x² − 6 and p₀ = 1. Use Newton’s method to find p₂ .

In [4]:
from sympy import *
NewtonsMethod(1, pow(x, 2) - 6, output="table", stopping_criterion="iteration", iteration=2, plot=True)
Out[4]:
p
0 1.0
1 3.5
2 2.607142857

2. Let f(x) = −x³ − cosx and p₀ = −1. Use Newton’s method to find p₂ . Could p₀ = 0 be used?

In [5]:
from sympy import *

x = Symbol('x')
NewtonsMethod(-1, pow(x, 3) - cos(x), output="table", stopping_criterion="iteration", iteration=2, plot=True)
Out[5]:
p
0 -1.0
1 -0.2864111184
2 -27.27238899
In [47]:
from sympy import *

x = Symbol('x')
NewtonsMethod(0, pow(x, 3) - cos(x), output="table", stopping_criterion="iteration", iteration=2)
---------------------------------------------------------------------------
ZeroDivisionError                         Traceback (most recent call last)
<ipython-input-47-12f7e0c09919> in <module>()
      2 
      3 x = Symbol('x')
----> 4 NewtonsMethod(0, pow(x, 3) - cos(x), output="table", stopping_criterion="iteration", iteration=2)

<ipython-input-42-4b17c606d1d1> in NewtonsMethod(p0, equation, tolerance, max_iterations, output, stopping_criterion, iteration)
     35 
     36     for i in range(1, max_iterations):
---> 37         p = p0 - (f(p0)/df(p0))
     38 
     39         if (output == "table"):

ZeroDivisionError: float division by zero
In [52]:
from sympy import *

x = Symbol('x')
equation = pow(x, 3) - cos(x)
df = lambdify(x, equation.diff(x))

df(0)
Out[52]:
0.0

No, we can't use p₀ = 0, because if we do, then derivate of f returns 0 and then we try to divide by that 0 in the formula to obtain p₁ , which causes a division by zero error.

5. Use Newton’s method to find solutions accurate to within 10⁻⁴ for the following problems.

a. x³ − 2x² − 5 = 0, [1,4]

In [6]:
from sympy import *

x = Symbol('x')
NewtonsMethod(p0=2, equation=pow(x, 3) - 2*pow(x, 2) - 5, tolerance=pow(10, -4), output="table", plot=True)
Out[6]:
p
0 2.0
1 3.25
2 2.811036789
3 2.697989502
4 2.690677153
5 2.690647449

b. x³ + 3x² − 1 = 0, [−3,−2]

In [8]:
from sympy import *

x = Symbol('x')
NewtonsMethod(p0=-3, equation=pow(x, 3) + 3*pow(x, 2) - 1, tolerance=pow(10, -4), output="table", plot=True)
Out[8]:
p
0 -3.0
1 -2.888888889
2 -2.879451567
3 -2.879385245

c. x − cosx = 0, [0,π/2]

In [9]:
from sympy import *

x = Symbol('x')
NewtonsMethod(p0=0, equation=x - cos(x), tolerance=pow(10, -4), output="table", plot=True)
Out[9]:
p
0 0.0
1 1.0
2 0.7503638678
3 0.7391128909
4 0.7390851334

d. x − 0.8 − 0.2sinx = 0, [0,π/2]

In [10]:
from sympy import *

x = Symbol('x')
NewtonsMethod(p0=0, equation=x - 0.8 - 0.2*sin(x), tolerance=pow(10, -4), output="table", plot=True)
Out[10]:
p
0 0.0
1 1.0
2 0.9644529683
3 0.964333889
4 0.9643338877