Reputation: 167
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:
Upvotes: 0
Views: 109
Reputation: 3001
Newton
function shouldn't yield and return at the same timeThis 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()
Upvotes: 1