Did you ever want to find the roots of some polynomial ?
Here is 'Newton's Method' ... to the rescue.
Roots by Newton's (iterative) method of successive (improved) approximations
Xn = Xp - f(Xp)/f'(Xp) i.e. Xnext = Xprev - f'(Xprev)/f(Xprev)
.
./ |f(Xp)
./ |
./ |
. /__|
. Xn Xp Method to find Xn is derived from:
f(Xp)/(Xp-Xn) = rise/run = slope = f'(Xp)
# newtonsRoots.py # this ver. 2013-02-01 #
header = \
'''Demo program linked from 'Beginning Calculus' by dwzavitz
Finds roots and max's and min's of polynomial using Newton's method
Example used here is:
f(x) = x**5 - 4*x**4 - 5*x**3 + 20*x**2 + 4*x - 16
testData = [-10, -1.1, -0.1, 0.1, 1.1, 10]
'''
testData = [-10, -1.1, -0.1, 0.1, 1.1, 10] # x[i] to use in testData
print( header )
print( '\nloops \t\t prev x \t new x')
def f( x ):
return x**5 - 4*x**4 - 5*x**3 + 20*x**2 + 4*x - 16
def df( x ):
return 5*x**4 - 16*x**3 - 15*x**2 + 40*x + 4
def getx1( x ):
return x - f(x)/df(x)
import math
err = 10**-15
def isEqual( a, b ):
if 1 - err <= (a/b) <= 1 + err:
return True
return False
def inList( v, lst ):
for item in lst:
if isEqual( v, item ): return True
return False
roots = []
for x1 in testData:
if x1 == 0: continue
x0 = 0
count = 0
while not isEqual(x0,x1):
count += 1
x0 = x1
x1 = getx1(x1)
s = '{:<3} {:<20} {:<20}'
print( s.format(count, x0, x1) )
print( "Check if f({}) = {} = 0 (or nearly so?)".format( x1, f(x1) ) )
if not inList( x1, roots ):
roots.append( x1 )
input("Press 'Enter' to continue ... ")
print()
print("\n\nNow ... finding roots of df = f2 (to find max and min) ...")
extremes = []
def inList( v, lst ):
for item in lst:
if isEqual( v, item[1] ): return True
return False
def f2( x ):
return 5*x**4 - 16*x**3 - 15*x**2 + 40*x + 4
def df2( x ):
return 20*x**3 - 48*x**2 - 30*x + 40
def getx2( x ):
return x - f2(x)/df2(x)
def maxMin( x ):
if x < 0: return "=> a Max at"
elif x > 0: return "=> a Min at"
else: return "=> an inflection at"
for x1 in testData:
if x1 == 0: continue
x0 = 0
count = 0
while not isEqual(x0,x1): #math.fabs(x1-x0) > 10**-15: # x0/x1 != 1:
count += 1
x0 = x1
x1 = getx2(x1)
s = '{:<3} {:<20} {:<20}'
print( s.format(count, x0, x1) )
print( "Check if f2({}) = {} = 0 (or nearly so?)".format( x1, f2(x1) ) )
tmp = df2(x1)
s2 = '{:.2f}'
print( "Value of df2({}) = {}".format(x1, tmp), \
maxMin(tmp), s2.format(x1) )
if not inList( x1, extremes ):
extremes.append( [maxMin(tmp), x1, f(x1)] )
input("Press 'Enter' to continue ... ")
print()
numPoints = 50
intervalWidth = 9
end = intervalWidth/2
dw = intervalWidth/numPoints
s = "points (x,y), on curve y = f(x), with x in [{:.2f},{:.2f}] ...\n"
print( "Table of", numPoints, s.format(-end, end-dw) )
points = []
# create table #
for p in range( numPoints ):
v = -end + p*dw
points.append( [v, f(v)] )
def getSign(x):
if x<0: return 0
return 1
# show table #
pOld = points[0][1]
n = 0
for p in points:
print( '({:5.2f},{:7.1f})'.format(p[0],p[1]), end = '' )
n += 1
if n == 5:
print()
n = 0
else:
print( ' ', end = '' )
if getSign(p[1]) != getSign(pOld): print( "<- sign changed\n" )
pOld = p[1]
print( "\nRoots:" )
print( roots )
print( "MaxMin:" )
for mm in extremes:
print( mm )
input("Press 'Enter' to continue ... ")
The above, a little more polished ...
# newtonsRoots5.py # this ver. 2013-02-01 #
header = \
'''Demo program linked from 'Beginning Calculus' by dwzavitz
Finds roots and max's and min's of polynomial using Newton's method
Example used here is:
f(x) = x**5 - 4*x**4 - 5*x**3 + 20*x**2 + 4*x - 16
testData = [x/10 for x in range(-101, 102)]
'''
testData = [x/10 for x in range(-101, 102)] # x[i] to use in testData #
print( header )
def f( x ):
return x**5 - 4*x**4 - 5*x**3 + 20*x**2 + 4*x - 16
def df( x ):
return 5*x**4 - 16*x**3 - 15*x**2 + 40*x + 4
def getx1( x ):
return x - f(x)/df(x)
import math
err = 10**-15
def isEqual( a, b ):
if 1 - err <= (a/b) <= 1 + err:
return True
return False
def inList( v, lst ):
for item in lst:
if isEqual( v, item ): return True
return False
roots = []
for x1 in testData:
if x1 == 0: continue
x0 = 0
count = 0
while not isEqual(x0,x1):
count += 1
x0 = x1
x1 = getx1(x1)
if not inList( x1, roots ):
roots.append( x1 )
extremes = []
def inList( v, lst ):
for item in lst:
if isEqual( v, item[1] ): return True
return False
def f2( x ):
return 5*x**4 - 16*x**3 - 15*x**2 + 40*x + 4
def df2( x ):
return 20*x**3 - 48*x**2 - 30*x + 40
def getx2( x ):
return x - f2(x)/df2(x)
def maxMin( x ):
if x < 0: return "=> a Max at"
elif x > 0: return "=> a Min at"
else: return "=> an inflection at"
for x1 in testData:
if x1 == 0: continue
x0 = 0
count = 0
while not isEqual(x0,x1):
count += 1
x0 = x1
x1 = getx2(x1)
tmp = df2(x1)
if not inList( x1, extremes ):
extremes.append( [maxMin(tmp), x1, f(x1)] )
input("Press 'Enter' to continue ... ")
print()
roots.sort()
numPoints = 50
intervalWidth = roots[-1]+1 - roots[0]
end = roots[0] - 1
dw = intervalWidth/numPoints
s = "points (x,y), on curve y = f(x), with x in [{:.2f},{:.2f}] ...\n"
print( "Table of", numPoints+5, s.format(end, end+(numPoints+4)*dw) )
points = []
# create table #
for p in range( numPoints+5 ):
v = end + p*dw
points.append( [v, f(v)] )
def getSign(x):
if x<0: return 0
return 1
# show table #
pOld = points[0][1]
n = 0
for p in points:
print( '({:5.2f},{:7.1f})'.format(p[0],p[1]), end = '' )
n += 1
if n == 5:
print()
n = 0
else:
print( ' ', end = '' )
if getSign(p[1]) != getSign(pOld): print( "<- sign changed\n" )
pOld = p[1]
print( "\nRoots:" )
print('[', end = '' )
for i, r in enumerate(roots):
print( '{:.3f}'.format(r), end = '' )
if i < len(roots)-1: print(', ', end = '' )
else: print(']')
print( "\nMaxMin:" )
from operator import itemgetter
extremes.sort(key=itemgetter(1))
for mm in extremes:
print( mm )
input("Press 'Enter' to continue ... ")
The above, but a more general solution ...
# newtonsRoots6.py # this ver. 2013-02-01 #
header = \
'''Demo program linked from 'Beginning Calculus' by dwzavitz
Finds roots and max's and min's of polynomial using Newton's method
Example used here is:
f(x) = x**6 - 4*x**5 - 5*x**4 + 20*x**3 + 4*x**2 - 16*x
testData = [x/10 for x in range(-1001, 1002)]
'''
testData = [x/10 for x in range(-1001, 1002)] # x[i] to use in testData #
print( header )
def f( x ):
return x**6 - 4*x**5 - 5*x**4 + 20*x**3 + 4*x**2 - 16*x
def df( x ):
return 6*x**5 - 20*x**4 - 20*x**3 + 60*x**2 + 8*x -16
def getx1( x ):
return x - f(x)/df(x)
import math
err = 10**-15
def isEqual( a, b ):
if b != 0:
if 1 - err <= (a/b) <= 1 + err:
return True
return False
else:
if abs(b-a) < err: return True
else: return False
def inList( v, lst ):
for item in lst:
if isEqual(v, item): return True
return False
roots = []
for x1 in testData:
if x1 == 0: continue
x0 = 0
count = 0
while not isEqual(x0,x1):
count += 1
x0 = x1
x1 = getx1(x1)
if not inList( x1, roots ):
roots.append( x1 )
extremes = []
def inList( v, lst ):
for item in lst:
if isEqual(v, item[1]): return True
return False
def f2( x ):
return 6*x**5 - 20*x**4 - 20*x**3 + 60*x**2 + 8*x -16
def df2( x ):
return 30*x**4 - 80*x**3 - 60*x**2 + 120*x + 8
def getx2( x ):
return x - f2(x)/df2(x)
def maxMin( x ):
if x < 0: return "=> a Max at"
elif x > 0: return "=> a Min at"
else: return "=> an inflection at"
for x1 in testData:
if x1 == 0: continue
x0 = 0
count = 0
while not isEqual(x0,x1):
count += 1
x0 = x1
x1 = getx2(x1)
tmp = df2(x1)
if not inList( x1, extremes ):
extremes.append( [maxMin(tmp), x1, f(x1)] )
input("Press 'Enter' to continue ... ")
print()
roots.sort()
numPoints = 50
intervalWidth = roots[-1]+1 - roots[0]
end = roots[0] - 1
dw = intervalWidth/numPoints
s = "points (x,y), on curve y = f(x), with x in [{:.2f},{:.2f}] ...\n"
print( "Table of", numPoints+5, s.format(end, end+(numPoints+4)*dw) )
points = []
# create table #
for p in range( numPoints+5 ):
v = end + p*dw
points.append( [v, f(v)] )
def getSign(x):
if x<0: return 0
return 1
# show table #
pOld = points[0][1]
n = 0
for p in points:
print( '({:5.2f},{:7.1f})'.format(p[0],p[1]), end = '' )
n += 1
if n == 5:
print()
n = 0
else:
print( ' ', end = '' )
if getSign(p[1]) != getSign(pOld): print( "<- sign changed\n" )
pOld = p[1]
print( "\nRoots:" )
print('[', end = '' )
for i, r in enumerate(roots):
print( '{:.3f}'.format(r), end = '' )
if i < len(roots)-1: print(', ', end = '' )
else: print(']')
print( "\nMaxMin:" )
from operator import itemgetter
extremes.sort(key=itemgetter(1))
for mm in extremes:
print( mm )
input("Press 'Enter' to continue ... ")
Just find nth root of some value (using Newton's Method) ...
# newtonsRootsNth.py # this ver. 2013-02-01 #
# nth root of value demo #
err = 10**-15
header = \
'''Demo program linked from 'Beginning Calculus' by dwzavitz
1. Input positive value to find root ...
2. Input integer n > 1 to find nth root (of value entered above)
Roots by Newton's (iterative) method of successive (improved) approximations
Xn = Xp - f(Xp)/f'(Xp) i.e. Xnext = Xprev - f'(Xprev)/f(Xprev)
.
./|f(Xp)
./ |
./ |
. /___|
. Xn Xp Method to find Xn is derived from:
f(Xp)/(Xp-Xn) = rise/run = slope = f'(Xp)
So to find 4th root of 5,
Use/solve f(x) = x**4 - 5 = 0
'''
print( header )
print()
#import math
def f( x, n, v ):
try:
return x**n - v ## example 4th root of 5 ==> solve: x**4 - 5 = 0 ##
except:
print( "Error in f..." )
return 0
def df( x, n ):
try:
return n*x**(n-1)
except:
print( "Error in df..." )
return 1
def getx1( x, n, v ):
return x - f(x,n,v)/df(x,n)
def more():
return input("More (y/n) ? " ).lower() != 'n'
def takeIn():
done = False
while not done:
floatErr = True
overflowErr = True
try:
v = eval(input( "Enter number to find nth root : " ) )
if v < 0:
raise
floatErr = False
n = int(input( "Enter n for nth root to find : " ) )
if n < 2:
overflowErr = False
raise
case = 0
if v < 1:
if n*v >= 1:
test = v - 1/n - (v - 1/n)**n - v
case = 1
else:
test = v*n - (v*n)**n - v
case = 2
elif v > 1:
test = (1+v/n) - (1+v/n)**n - v
case = 3
overflowErr = False
done = True
except:
if floatErr:
print( "Entry error! Positive number expected ... Try again ... " )
elif overflowErr:
if case == 1: x = v - 1/n
elif case == 2: x = n*v
else: x = 1 + v/n
print( "Overflow error! {}**{} causes overflow ... Try again ... ".format(x,n) )
else:
print( "Entry error! Positive integer > 1 expected ... Try again ... " )
return v, n
while 1:
v, n = takeIn()
if v < 1:
if n*v >= 1:
x1 = v - 1/n # (n*v-1)/n
else:
x1 = n*v
else:
x1 = 1 + v/n
x0 = 0
count = 0
ok = True
while ok:
count += 1
x0 = x1
try:
x1 = getx1(x1,n,v)
#print( count, x0, x1 )
except:
print( "Some error occured ... perhaps overflow ...")
ok = False
if (1 - err) <= x0/x1 <= (1 + err):
break
if count > 10000:
print( "Loops > 10,000 ERROR!" )
break
if n > 3: th = 'th'
elif n == 3: th = 'rd'
else: th = 'nd' # since n == 2 here #
print( '{}{}_root({}) = {} found in {} loops.'.\
format(n, th, v, x1, count) )
print( "Check if f({}) = {} = 0.0 (or nearly so?)".\
format( x1, f(x1,n,v) ) )
if not more():
break
print()
input("Press 'Enter' to continue/exit ... ")
