I'm trying to implement Newton's method for fun in Python but I'm having a problem conceptually understanding the placement of the check.
A refresher on Newton's method as a way to approximate roots via repeated linear approximation through differentiation:
I have the following code:
# x_1 = x_0 - (f(x_0)/f'(x_0))
# x_n+1 - x_n = precision
def newton_method(f, f_p, prec=0.01):
x=1
x_p=1
tmp=0
while(True):
tmp = x
x = x_p - (f(x_p)/float(f_p(x_p)))
if (abs(x-x_p) < prec):
break;
x_p = tmp
return x
This works, however if I move the if
statement in the loop to after the x_p = tmp
line, the function ceases to work as expected. Like so:
# x_1 = x_0 - (f(x_0)/f'(x_0))
# x_n+1 - x_n = precision
def newton_method(f, f_p, prec=0.01):
x=1
x_p=1
tmp=0
while(True):
tmp = x
x = x_p - (f(x_p)/float(f_p(x_p)))
x_p = tmp
if (abs(x-x_p) < prec):
break;
return x
To clarify, function v1 (the first piece of code) works as expected, function v2 (the second) does not.
Why is this the case?
Isn't the original version essentially checking the current x
versus the x
from 2 assignments back, rather than the immediately previous x
?
Here is the test code I am using:
def f(x):
return x*x - 5
def f_p(x):
return 2*x
newton_method(f,f_p)
EDIT
I ended up using this version of the code, which forgoes the tmp
variable and is much clearer for me, conceptually:
# x_1 = x_0 - (f(x_0)/f'(x_0))
# x_n+1 - x_n = precision
def newton_method(f, f_p, prec=0.01):
x=1
x_p=1
tmp=0
while(True):
x = x_p - (f(x_p)/float(f_p(x_p)))
if (abs(x-x_p) < prec):
break;
x_p = x
return x
Let x[i]
be the new value to be computed in an iteration.
The statement x = x_p - (f(x_p)/float(f_p(x_p)))
translates to:
x[i] = x[i-2] - f(x[i-2])/f'(x[i-2])
- 1
But according to the actual mathematical formula, it should have been this:
x[i] = x[i-1] - f(x[i-1])/f'(x[i-1])
Similarly, x[i-1] = x[i-2] - f(x[i-2])/f'(x[i-2])
- 2
Comparing 1 and 2, we can see that the x[i]
in 1 is actually x[i-1]
according to the math formula.
The main point to note here is that x
and x_p
are always one iteration apart. That is, x
is the actual successor to x_p
, unlike what it might seem by just looking at the code.
Hence, it is working correctly as expected.
Just like the above case, the same thing happens at the statement x = x_p - (f(x_p)/float(f_p(x_p)))
.
But by the time we reach if (abs(x-x_p) < prec)
, x_p
had changed its value to temp
= x
= x[i-1]
.
But as deduced in the case of version 1, x
too is x[i-1]
rather than x[i]
.
So, abs(x - x_p)
translates to abs(x[i-1] - x[i-1])
, which will turn out to be 0, hence terminating the iteration.
The main point to note here is that x
and x_p
are actually the same values numerically, which always results in the algorithm terminating after just 1 iteration itself.