Python - Implementing a numerical equation solver (Newton-Raphson)

Question

I am warning you, this could be confusing, and the code i have written is more of a mindmap than finished code..

I am trying to implement the Newton-Raphson method to solve equations. What I can't figure out is how to write this

equation in Python, to calculate the next approximation (xn+1) from the last approximation (xn). I have to use a loop, to get closer and closer to the real answer, and the loop should terminate when the change between approximations is less than the variable h.

1. How do I write the code for the equation?

2. How do I terminate the loop when the approximations are not changing anymore?

   Calculates the derivative for equation f, in point x, with the accuracy h (this is used in the equation for solve())

   def derivative(f, x, h):
       deriv = (1.0/(2*h))*(f(x+h)-f(x-h))
       return deriv
   

   The numerical equation solver

   Supposed to loop until the difference between approximations is less than h

   def solve(f, x0, h):
       xn = x0
       prev = 0
   
       while ( approx - prev > h):
            xn = xn - (f(xn))/derivative(f, xn, h)
   
       return xn
   

Answer 1

Here is the implementation of a N-R solver expanding what you wrote above (complete, working). I added a few extra lines to show what is happening...

def derivative(f, x, h):
      return (f(x+h) - f(x-h)) / (2.0*h)  # might want to return a small non-zero if ==0

def quadratic(x):
    return 2*x*x-5*x+1     # just a function to show it works

def solve(f, x0, h):
    lastX = x0
    nextX = lastX + 10* h  # "different than lastX so loop starts OK
    while (abs(lastX - nextX) > h):  # this is how you terminate the loop - note use of abs()
        newY = f(nextX)                     # just for debug... see what happens
        print "f(", nextX, ") = ", newY     # print out progress... again just debug
        lastX = nextX
        nextX = lastX - newY / derivative(f, lastX, h)  # update estimate using N-R
    return nextX

xFound = solve(quadratic, 5, 0.01)    # call the solver
print "solution: x = ", xFound        # print the result

output:

f( 5.1 ) =  27.52
f( 3.31298701299 ) =  6.38683083151
f( 2.53900845771 ) =  1.19808560807
f( 2.30664271935 ) =  0.107987672721
f( 2.28109300639 ) =  0.00130557566462
solution: x =  2.28077645501

Edit - you could also check the value of newY and stop when it is "close enough to zero" - but usually you keep this going until the change in x is <=h (you can argue about the value of the = sign in a numerical method - I prefer the more emphatic < myself, figuring that one more iteration won't hurt.).

Answer 2

If code under 'try:' cannot be implemented and the compiler is given the error 'ZeroDivisionError' then it will execute the code under 'except ZeroDivisionError:'. Although, if you'd like to account for another compiler exception 'XYZ' with a specific code implementation then add an additional 'except XYZ:"

 try:
    nextX = lastX - newY / derivative(function, lastX, h)
 except ZeroDivisionError:
    return newY