fr14
fr14

Reputation: 167

Graphing diagram In Python

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 on the interval [0,1]. Meaning the number of iterations on the x axis with the corresponding absolute error on the y axis

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 output and graph. Is my output correct? should the graph look like this? All help is greatly appreciated! The number of iterations from my code is 3

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 iteration_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
    return x, iteration_counter
def f(x):
    return (math.cos(x)-math.sin(x))
def dfdx(x):
    return (-math.sin(x)-math.cos(x))

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

which produces the following output:

array([[0.00000000e+00, 4.74646213e-03],
   [1.00000000e+00, 1.78222779e-08]])

and finally:

import numpy as np
import matplotlib.pyplot as plt

data = np.array(list(Newton(f,dfdx, 1, 10e-14)))
plt.plot(data[:,0], data[:,1])
plt.yscale('log')
plt.show()

which produces the graph:

enter image description here

Upvotes: 0

Views: 109

Answers (1)

  1. Your Newton function shouldn't yield and return at the same time
  2. Use a more slowly converging function to test your results

This is what I would do:

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


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 iteration_counter, x, 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

def f(x):
    return x ** 2 - 1.34

def dfdx(x):
    return 2 * x

data = np.array(list(newton(f, dfdx, 10, 10e-14)))

# plt.plot(data[:, 0], data[:, 1]) # x-axis: iteration, y-axis: x values
plt.plot(data[:, 0], data[:, 2]) # x-axis: iteration, y-axis: f(x) values
plt.yscale('log')
plt.show()

enter image description here

Upvotes: 1

Related Questions