1

In the following code I have implemented Newtons method in Python.

import math
def Newton(f, dfdx, x, eps):
    f_value = f(x)
    iteration_counter = 0
    while abs(f_value) > eps and iteration_counter < 100:
        try:
            x = x - float(f_value)/dfdx(x)
        except ZeroDivisionError:
            print ("Error! - derivative zero for x = ", x)
            sys.exit(1)     # Abort with error
        f_value = f(x)
        iteration_counter += 1

    # Here, either a solution is found, or too many iterations
    if abs(f_value) > eps:
        iteration_counter = -1
    return x, iteration_counter
def f(x):
    return (math.cos(x)-math.sin(x))
def dfdx(x):
    return (-math.sin(x)-math.cos(x))
solution, no_iterations = Newton(f, dfdx, x=1, eps=1.0e-14)
if no_iterations > 0:    # Solution found
    print ("Number of function calls: %d" % (1 + 2*no_iterations))
    print ("A solution is: %f" % (solution))
else:
    print ("Solution not found!")

However now I am looking to plot the convergence diagram on that same interval. This would be the absolute error as a function of the number of iterations.

I attempted to make an iterable that each time yields a 2-tuple with the absolute error, and the iteration. Here is my code below, with the error I am getting from it included also,

import math
def Newton(f, dfdx, x, eps):
    f_value = f(x)
    iteration_counter = 0
    while abs(f_value) > eps and iteration_counter < 100:
        try:
            x = x - float(f_value)/dfdx(x)
            yield interation_counter, abs(f(x))
        except ZeroDivisionError:
            print ("Error! - derivative zero for x = ", x)
            sys.exit(1)     # Abort with error
        f_value = f(x)
        iteration_counter += 1

    # Here, either a solution is found, or too many iterations
    if abs(f_value) > eps:
        iteration_counter = -1
def f(x):
    (math.cos(x)-math.sin(x))
def dfdx(x):
    (-math.sin(x)-math.cos(x))

With this I tried to put the result into an array so I could graph my results

import numpy as np
np.array(list(Newton(f,dfdx, 1,10e-4)))

However I am getting the following error:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-20-9378f4e2dbe3> in <module>()
      1 import numpy as np
----> 2 np.array(list(Newton(f,dfdx, 1,10e-4)))

<ipython-input-19-40b67c2c3121> in Newton(f, dfdx, x, eps)
      4     f_value = f(x)
      5     iteration_counter = 0
----> 6     while abs(f_value) > eps and iteration_counter < 100:
      7         try:
      8             x = x - float(f_value)/dfdx(x)

TypeError: bad operand type for abs(): 'NoneType'
fr14
  • 167
  • 11

1 Answers1

3

I suggest storing each value x and corresponding output f for each iteration in two respective arrays, and then returning each array. If your function is correct (you should know if it converges), then this should be easy to do. From there, just plot the array.

I did this a couple months ago. You can see how I worked through it here: How to tell if Newtons-Method Fails

So A couple things to note here.

1) I did not check to make sure that your function is converging correctly or not.

2) Newtons method commonly has a very large step size and you usually won't get any pretty visualizations for its convergences. It's not uncommon for there to be <10 iterations, and they usually don't follow a smooth path to the convergence point either (once again, large step size). Notice that for the code I implemented, there were only 3 iterations. But that could just be because your function is wrong. Frankly, I'm not sure

import numpy as np
import math
import matplotlib.pyplot as plt

def Newton(f, dfdx, x, eps):
    xstore=[]
    fstore=[]
    f_value = f(x)
    iteration_counter = 0
    while abs(f_value) > eps and iteration_counter < 100:
        try:
            x = x - float(f_value)/dfdx(x)
        except ZeroDivisionError:
            print ("Error! - derivative zero for x = ", x)
            sys.exit(1)     # Abort with error
        f_value = f(x)
        xstore.append(x)
        fstore.append(f_value)
        iteration_counter += 1

    # Here, either a solution is found, or too many iterations
    if abs(f_value) > eps:
        iteration_counter = -1

    return x, iteration_counter,xstore,fstore

def f(x):
    return (math.cos(x)-math.sin(x))
def dfdx(x):
    return (-math.sin(x)-math.cos(x))

solution, no_iterations,xvalues,fvalues = Newton(f, dfdx, x=1, eps=1.0e-14)

if no_iterations > 0:    # Solution found
    print ("Number of function calls: %d" % (1 + 2*no_iterations))
    print ("A solution is: %f" % (solution))
else:
    print ("Solution not found!")

x = np.array([i for i in xvalues])
f = np.array(fvalues)
fig = plt.figure()
plt.scatter(x,f,label='Newton')
plt.legend()

enter image description here

RocketSocks22
  • 391
  • 1
  • 4
  • 20
  • could you show me how this is done with my code? I can accept the answer then – fr14 Nov 04 '18 at 22:04
  • thanks for your edit. Is this a convergence plot that shows the absolute error as a function of the number of iterations? right now my code only needs 7 iterations to find my solution – fr14 Nov 04 '18 at 22:28
  • I believe this just shows the answer, not the absolute error as a function of the number of iterations – fr14 Nov 04 '18 at 22:41
  • No, this is only plotting the values your Newtons method iterates over (e.g. x_1, x_2, x_3, etc. w/ the corresponding objective function output for a given x) Also, your code only iterates 3 times. I do not know why you defined the number of function calls as `1 + 2*no_iterations`. I'm not following what you mean when you say 'absolute error as a function of the number of iterations.' Regardless, you should get the gist of how it's possible to visualize your output for whatever you are trying to output – RocketSocks22 Nov 04 '18 at 22:43
  • how would I fix the number of iterations then? – fr14 Nov 04 '18 at 22:52
  • i guess just using no_iterations would fix this – fr14 Nov 04 '18 at 22:53
  • I tried with the implementation in `scipy` and I get the same answer after three evaluations of `f`, so the OP's code seems to be right – Sam Mason Nov 04 '18 at 23:04
  • @Sam Mason so its 3 iterations? I had the output statement wrong – fr14 Nov 04 '18 at 23:18
  • @fr14 I suppose it depends on your definition of what an "iteration" is, i.e. whether you want to count invocations of `f` or something else. I also left the error tolerance as default in scipy, if you want to do more in depth testing that's up to you – Sam Mason Nov 05 '18 at 14:06