 CoCalc Public Filesweek 5 / assignment / FindRoots.ipynb
Authors: Boyuan Qi, Rajvir Sidhu
Views : 31
Compute Environment: Ubuntu 20.04 (Default)

# Assignment 2: Polynomial Root Finding

## Markdown

Give the address (url) for a website that has documentation for markdown, including tables (i.e. a webpage, not a pdf or other document, that describes how you use markdown to create a table). Enter your answer as a variable named website='http://url.goes.here'

Marks: 1
In [ ]:
# YOUR CODE HERE

markdownWebsite='https://www.markdownguide.org/cheat-sheet'

Test your code from above here(1 point), ID: validate_website
In [ ]:
#test that it's a suitable website link you've given


Work out the solutions to the following questions long-hand (i.e. you should not be coding anything at this point). In the first box, set the appropriate variable equal to your answer (there's no need to define any functions). For example, to answer the last question, you will simply type polyRoots=[4-math.sqrt(3),4,4+math.sqrt(3)]

1. Define the polynomial $p(x)=x^3-12x^2+45x-52$. How do you express $p(x)$ as a list, such as will be input to our algorithms? Set the variable polyX equal to this.
2. What is the derivative, $p'(x)$ (as a list called polyDeriv) and the root(s) of $p'(x)$ (as a list, in increasing order, called derivRoots)?
3. What is the Lagrange bound as applied to the polynomial $p(x)$? Your answer should be a real, positive number, called lagrange. (In case you missed it, this was an exercise in the workshops)
4. Pretend that you don't know the roots of $p(x)$, but you know the roots of $p'(x)$. Use this information to explain why we know that there is no more than one root in each of the ranges: $(-109,3),(3,5),(5,109)$ You should do this without evaluating $p(x)$ at any point. Write your solution in the free text box below (this will be manually graded).
5. What are the roots of $p(x)$? Give your answer as a list in increasing order called polyRoots. Make sure that all the roots of $p(x)$ appear in exactly one of the above ranges.
Marks: 5
In [ ]:
import math

polyX=[-52,45,-12,1]
polyDeriv=[45,-24,3]
derivRoots=[3,5]
lagrange=109
polyRoots=[4-math.sqrt(3),4,4+math.sqrt(3)]

Test your code from above here(1 point), ID: polyRoots_correct
In [ ]:
#the question told you how to answer this, so no point in hiding the test
assert polyRoots==[4-math.sqrt(3),4,4+math.sqrt(3)]

Test your code from above here(1 point), ID: polyX_correct
In [ ]:
#check that polyX is correct

Test your code from above here(1 point), ID: polyDeriv_correct
In [ ]:
#check that polyDeriv is correct

Test your code from above here(1 point), ID: derivRoots_correct
In [ ]:
#check that derivRoots is correct

Test your code from above here(1 point), ID: lagrange_bound_correct
In [ ]:
#check that lagrange is correct


This is where you give your answer to (4) above: Explain why we know that there is no more than one root in each of the ranges: $(-109,3),(3,5),(5,109)$ You should do this without evaluating $p(x)$ at any point.

Marks: 2

Given any two roots b>a of p(x), f is continuous and differentiable on [a,b]. Using Rolle's theorem, there exists exactly one root, 'k', such that f'(k)=0 and a<c<b. Due to p'(x) having 2 roots, p(x) therefore has at most 3 roots. The roots of p'(x) are 3 and 5 which means p(x) must have a root within the interval [3,5]. Also, the limits of the function are 109 and -109, so if there are 3 roots, then exactly one root must exist in the each of the intervals [-109,3] and [5,109]. This means there that no more than one root can exist in each of the intervals [-109,3],[3,5] and [5,109]

## Collecting our work so far

In working through the lab sessions, we have built up a library of subroutines that will be useful for finding the roots of a polynomial. Those that we have already produced are given first.

In [ ]:
from math import sqrt

def EvaluatePoly(polynomial,x):
'''evaluate a polynomial at a value x
input: polynomial as a list, value x to evaluate at
output: value of polynomial'''
return(sum(val*x**n for n,val in enumerate(polynomial)))
def DifferentiatePoly(polynomial):
'''return the derivative of polynomial'''
return [(i+1)*val for i,val in enumerate(polynomial[1:])]
def EnsureStandardForm(polynomial):
'''make sure that the highest order of the polynomial has non-zero coefficient by removing all higher order zeros'''
while not polynomial[-1]:
del polynomial[-1]
return polynomial
def DescartesSigns(polynomial):
'''Descartes rules of signs returns an upper bound on the number of positive roots and the number of negative roots'''
PosList=list(polynomial)
NegList=[((-1)**i)*val for i,val in enumerate(polynomial)]
while 0 in PosList:
PosList.remove(0)
NegList.remove(0)
return([sum([i*i<0 for i in zip(PosList[1:],PosList[:-1])]),sum([i*i<0 for i in zip(NegList[1:],NegList[:-1])])])
def FindZeroInInterval(polynomial,range_lower,range_upper):
'''perform an interval bisection on polynomial between the values range_lower<range_upper
return x such that polynomial(x)=0 (or a good enough approximation to it)'''
#this is slightly modified compared to what we presented.
assert range_lower<=range_upper,"second argument should be smaller than first"
accuracy=10**(-15)
value_lower=EvaluatePoly(polynomial,range_lower)
value_upper=EvaluatePoly(polynomial,range_upper)
if value_lower==0: # an exact root on the lower boundary
return(range_lower)
elif abs(value_upper/value_lower)<accuracy and abs(value_upper)<accuracy: #define this as being close enough to being a root on the upper boundary
return(range_upper)
elif abs(value_lower/value_upper)<accuracy and abs(value_lower)<accuracy: #close enough to root on lower boundary
return(range_lower)
if (range_lower==range_upper):#we already know it's not a root.
return('')
extent=2*(range_upper-range_lower)     #initialise with a dummy variable larger than anything that will ver be calculated

while range_upper-range_lower>max(abs(range_upper),abs(range_lower),1)*accuracy and extent>(range_upper-range_lower):#stopping criteria based on numerical accuracy
if value_upper*value_lower>0:      #we don't seem to have a root in this range
return('')
range_mid=(range_upper+range_lower)/2    #bisection
value_mid=EvaluatePoly(polynomial,range_mid)
extent=range_upper-range_lower

if value_mid==0:              #we have found the root
return range_mid
elif value_mid*value_lower>0: #the root is between value_mid and value_upper. redefine range and repeat
range_lower=range_mid
value_lower=value_mid
else:
range_upper=range_mid
value_upper=value_mid

return range_mid

return range_mid
def Smooth(mylist):
'''remove all the empty entries in a list'''
assert len(mylist)>0,'Should not smooth an empty list'
while '' in mylist:
mylist.remove('')
return mylist
def ZeroTest(mylist):
"fudge the fact that it's hard to tell if a very small number is actually 0"
return([(not abs(a)<10**(-14))*a for a in mylist])
def RootLocationBound(polynomial):
'use a Lagrange bound to crudely limit where the roots can be found'
return max(1,sum(abs(a/polynomial[-1]) for a in polynomial[:-1]))

def QuickTests():
#tests for the above code
assert EnsureStandardForm([1,1,1])==[1,1,1]
assert EnsureStandardForm([1,1,1,0,0,0])==[1,1,1]
assert [EvaluatePoly([24,-50,35,-10,1],i)-(24-50*i+35*i*i-10*i**3+i**4) for i in range(-5,5,1)]==[0,0,0,0,0,0,0,0,0,0]
assert DifferentiatePoly([24,-50,35,-10,1])==[-50,70,-30,4]
assert ZeroTest([-1,0.1,10**(-15),2])==[-1,0.1,0,2]
assert DescartesSigns([-1,-1,1,1])==[1,2]
assert DescartesSigns([-1,0,0,1])==[1,0]
return True

assert QuickTests()


Let's say you were given a polynomial as a list polynomial and had found the roots of its derivative, as an ordered list: deriv_roots. Write a function root_bounds that returns a list of pairs of numbers between which there is no more than one root of the polynomial. The list should be ordered, going from smallest to largest, and the pairs of values should also be ordered smallest to largest.

Hint: If you did the free text box above (question 4) clearly and correctly, this should be straightforward.

Hint: Look at the first test to understand the input and output format.

Marks: 4
In :
def root_bounds(deriv_roots,upper):
'''return the pairs of numbers which bound individual roots (or none at all)
in: deriv_roots: a list of roots of the deriviative of a polynomial, in increasing order, repeated according to their multiplicity
in: upper. if x is a root of polynomial, |x|<=upper
'''

assert upper>=0
assert all([val1<=val2 for val1,val2 in zip (deriv_roots[:-1],deriv_roots[1:])])
assert isinstance (deriv_roots,list)
assert len(deriv_roots)>0  # This is needed to actually apply Rolle's Theorem
emptyList=[]   # A list which is used later on in the interation
if deriv_roots[-1]>upper:
for x in range(len(deriv_roots)-1):
listAdd=[deriv_roots[x],deriv_roots[x+1]]  # A nested listed will be added to the interval that contains the root
if deriv_roots[-1]<=upper:
lower=[-upper]
high=[upper]
z=list(deriv_roots)
bigList=lower+z+high  # The completed list which has the upper and lower bounds of the function
for x in range(len(bigList)-1):
listAdd=[bigList[x],bigList[x+1]]   # A nested listed will be added to the interval that contains the root
return emptyList

print(root_bounds([-1,0,0,1],5))

[[-5, -1], [-1, 0], [0, 0], [0, 1], [1, 5]]
Test your code from above here(2 points), ID: root_bounds_code_visible
In [ ]:
#test the root_bounds function
assert list(root_bounds([-1,0,0,1],5))==[[-5,-1],[-1,0],[0,0],[0,1],[1,5]] or list(root_bounds([-1,0,0,1],5))==[(-5,-1),(-1,0),(0,0),(0,1),(1,5)]

Test your code from above here(1 point), ID: root_bounds_with_float
In [ ]:
#a basic test

Test your code from above here(1 point), ID: root_bounds_preconditions
In [ ]:
#check some preconditions


The function FindZeroInInterval (above) is a little different to the code that we presented in the notes (given as OriginalFindZeroInInterval in cell below). Aside from the assert statement, how is it different? (explanation not required.)

To demonstrate the difference, give a set of parameters as input that gives a different outcome between the two functions. Give your answer as 3 variables poly, below and above.

Marks: 1

In [ ]:
# YOUR CODE HERE

poly=[3,1,0,0]
above=0
below=-3

Test your code from above here(2 points), ID: difference_FindZero
In [ ]:
def OriginalFindZeroInInterval ( polynomial , range_lower , range_upper ):
''' perform an interval bisection on polynomial between the values
range_lower < range_upper
return x such that polynomial (x)=0 (or a good enough approximation to it)
'''
value_lower = EvaluatePoly ( polynomial , range_lower )
value_upper = EvaluatePoly ( polynomial , range_upper )

while range_upper - range_lower >10**( -15) :
if value_upper * value_lower >0: #we don 't seem to have a root in this range
return('')
range_mid =( range_upper + range_lower )/2 # bisection
value_mid = EvaluatePoly ( polynomial , range_mid )
if value_mid ==0: #we have found the 0, so return itsposition
return range_mid
elif value_mid * value_lower >0: # crossing between range_mid andrange_upper
range_lower = range_mid
value_lower = value_mid
else : # crossing between range_lower and
range_mid
range_upper = range_mid
value_upper = value_mid
return range_mid

assert above>below
assert FindZeroInInterval(poly,below,above)!=OriginalFindZeroInInterval(poly,below,above)


## The Final Root Finder

Using the previous functions (and defining any additional functions that wish in the cell below), complete the function FindRealZeros which accepts a polynomial as a list and returns a list of the roots of the polynomial, in increasing order, appearing as many times as they are repeated.

You might structure your algorithm by considering:

• How did you work out the solution by hand?
• What bits of code do you have to hand that could replace some of the by-hand calculation.
• Are there any simple types of polynomial for which you can immediately give the root?
• Is there a minimum size of polynomial that has a solution?
• Can you give a set of ranges between which no more than 1 root of the polynomial lies, given the roots of the derivative?
• For a polynomial of a given degree, $n$, what's the maximum number of ranges there could be? How many real roots could there be for a degree $n$ polynomial?

You need to make use of the functions we've already defined, particularly FindZeroInInterval, DifferentiatePoly and RootLocationBound

In [ ]:
def EvaluatePoly(polynomial,x):
'''evaluate a polynomial at a value x
input: polynomial as a list, value x to evaluate at
output: value of polynomial'''
return(sum(val*x**n for n,val in enumerate(polynomial)))

def RootLocationBound(polynomial):
'use a Lagrange bound to crudely limit where the roots can be found'
return max(1,sum(abs(a/polynomial[-1]) for a in polynomial[:-1]))

def root_bounds(deriv_roots,upper):
'''return the pairs of numbers which bound individual roots (or none at all)
in: deriv_roots: a list of roots of the deriviative of a polynomial, in increasing order, repeated according to their multiplicity
in: upper. if x is a root of polynomial, |x|<=upper
'''

assert upper>=0
assert all([val1<=val2 for val1,val2 in zip (deriv_roots[:-1],deriv_roots[1:])])
assert isinstance (deriv_roots,list)
assert len(deriv_roots)>0  # This is needed to actually apply Rolle's Theorem
emptyList=[]   # A list which is used later on in the interation
if deriv_roots[-1]>upper:
for x in range(len(deriv_roots)-1):
listAdd=[deriv_roots[x],deriv_roots[x+1]]  # A nested listed will be added to the interval that contains the root
if deriv_roots[-1]<=upper:
lower=[-upper]
high=[upper]
z=list(deriv_roots)
bigList=lower+z+high  # The completed list which has the upper and lower bounds of the function
for x in range(len(bigList)-1):
listAdd=[bigList[x],bigList[x+1]]   # A nested listed will be added to the interval that contains the root
return emptyList

print(root_bounds([-1,0,0,1],5))

def EvaluatePoly(polynomial,x):
'''evaluate a polynomial at a value x
input: polynomial as a list, value x to evaluate at
output: value of polynomial'''
return(sum(val*x**n for n,val in enumerate(polynomial)))

def DifferentiatePoly(polynomial,order=1):
'''return the derivative of polynomial of order n'''
derivative=[(i+1)*val for i,val in enumerate(polynomial[1:])]
if order==1:
return derivative
else:
return DifferentiatePoly(derivative, order-1)

def DescartesSigns(polynomial):
'''Descartes rules of signs returns an upper bound on the number of positive roots and the number of negative roots'''
PosList=list(polynomial)
NegList=[((-1)**i)*val for i,val in enumerate(polynomial)]
while 0 in PosList:
PosList.remove(0)
NegList.remove(0)
return([sum([i*i<0 for i in zip(PosList[1:],PosList[:-1])]),sum([i*i<0 for i in zip(NegList[1:],NegList[:-1])])])

DescartesSigns([-52,45,-12,1])


In :
import numpy
from math import sqrt
import math

def EvaluatePoly(polynomial,x):
'''evaluate a polynomial at a value x
input: polynomial as a list, value x to evaluate at
output: value of polynomial'''
return(sum(val*x**n for n,val in enumerate(polynomial)))

def FindZeroInInterval(polynomial,range_lower,range_upper):
'''perform an interval bisection on polynomial between the values range_lower<range_upper
return x such that polynomial(x)=0 (or a good enough approximation to it)'''
#this is slightly modified compared to what we presented.
assert range_lower<=range_upper,"second argument should be smaller than first"
accuracy=10**(-15)
value_lower=EvaluatePoly(polynomial,range_lower)
value_upper=EvaluatePoly(polynomial,range_upper)
if value_lower==0: # an exact root on the lower boundary
return(range_lower)
elif abs(value_upper/value_lower)<accuracy and abs(value_upper)<accuracy: #define this as being close enough to being a root on the upper boundary
return(range_upper)
elif abs(value_lower/value_upper)<accuracy and abs(value_lower)<accuracy: #close enough to root on lower boundary
return(range_lower)
if (range_lower==range_upper):
return('')
extent=2*(range_upper-range_lower)     #initialise with a dummy variable larger than anything that will ver be calculated

while range_upper-range_lower>max(abs(range_upper),abs(range_lower),1)*accuracy and extent>(range_upper-range_lower):#stopping criteria based on numerical accuracy
if value_upper*value_lower>0:
return('no zero')
range_mid=(range_upper+range_lower)/2    #bisection
value_mid=EvaluatePoly(polynomial,range_mid)
extent=range_upper-range_lower

if value_mid==0:
return range_mid
elif value_mid*value_lower>0: #the root is between value_mid and value_upper. redefine range and repeat
range_lower=range_mid
value_lower=value_mid
else:
range_upper=range_mid
value_upper=value_mid

return range_mid

def RemoveTrailingZero(list):
"Remove the '.0' of a float to transform it in integer"
trail=[]
for x in list:
print(x)
if (x!=0 and x%1==0):
idk=int(x)
trail.append(idk)
elif x==0.0:
idk=0
trail.append(idk)
else:
trail.append(x)
return trail

def FindRealZeros(polynomial):
'''return the zeros of the polynomial in an ordered list'''

ReversePoly=list(reversed(polynomial)) #Reverse the order of the list 'polynomial' for numpy because numpy consider the first entry the lower degree
ListRoot=numpy.roots(ReversePoly) #list of root(s) of the polynomial found by numpy with a 'okay' precision
print(ListRoot)
newPrecision=[]
output=[]
Final=[]
nestedList=[]
withoutImaginary=[]

for x in range(len(ListRoot)):
newPrecision+=[numpy.around(ListRoot[x],decimals=7,out=None)]
for x in newPrecision: #get rid of complex numbers if present in the previous list
if x.imag==0: #keep the real numbers and remove the "0j" part
i=x.real
withoutImaginary.append(i)
for x in withoutImaginary: #Determine intervals around the roots with great precision
DetermineIntervals=[x-0.0000001,x+0.0000001]
output+=DetermineIntervals #Create a list whose entries are end points of intervals
i=0
while i<len(output): #Create a nested list
nestedList.append(output[i:i+2])
i+=2
print (nestedList)
for x in nestedList:
Zeropolynomial=FindZeroInInterval(polynomial,x,x) #The upper bound and lower bound correspond to the endpoints of each interval in the nested list
if Zeropolynomial =='no zero':
pass
else:
Final.append(Zeropolynomial)
IncreasingOrder=list(sorted(Final,key=float)) #Reorganized the roots found by increasing order
end=RemoveTrailingZero(IncreasingOrder) #In the case where one of the root/multiple roots are supposed integer(s), we need to remove the '.0' part

print(end)
return end



We will test if your solution works by comparing ideal solutions to your calculated solutions via the function EqualityWithinTolerance, which aims to tolerate numerical inaccuracy. Some examples are given below. There will be further tests that have to be passed.

Marks: 10
Test your code from above here(3 points), ID: root_finder
In :
#basic testing examples
def EqualityWithinTolerance(list1,list2,epsilon):
"""since we don't have perfect accuracy, cannot just test equality of two lists in our tests.
return true if all elements are close enough (within epsilon)"""
assert len(list1)==len(list2),'cannot compare lists of different lengths'
return sum([abs(i-i)<=epsilon for i in zip(list1,list2)])==len(list1)

#integer solutions
assert EqualityWithinTolerance(FindRealZeros([24,-50,35,-10,1]),[1,2,3,4],2*10**(-10))
assert EqualityWithinTolerance(FindRealZeros([-16,0,0,0,1]),[-2,2],10**(-10))
#irrational solutions
assert EqualityWithinTolerance(FindRealZeros([70,-60,12]),[(15-sqrt(15))/6,(15+sqrt(15))/6],10**(-10))
#no real solutions
assert EqualityWithinTolerance(FindRealZeros([16,0,0,0,1]),[],10**(-10))
assert EqualityWithinTolerance(FindRealZeros([16,0,1]),[],10**(-10)),"You should not be returning complex roots"
#roots very close to 0
assert EqualityWithinTolerance(FindRealZeros([-10**(-12),0,1]),[-10**(-6),10**(-6)],10**(-10))
#try a longer polynomial to make sure it's not just an explicit formula going up to degree 3 or 4
assert EqualityWithinTolerance(FindRealZeros([-362880, 1026576, -1172700, 723680, -269325, 63273, -9450, 870, -45,1]),[1,2,3,4,5,6,7,8,9],10**(-10))

[4. 3. 2. 1.] [[3.9999999, 4.0000001], [2.9999999, 3.0000001], [1.9999999, 2.0000001], [0.9999999, 1.0000001]] 1.0 2.0 3.0 4.0 [1, 2, 3, 4] [-2.00000000e+00+0.j 1.66533454e-16+2.j 1.66533454e-16-2.j 2.00000000e+00+0.j] [[-2.0000001, -1.9999999], [1.9999999, 2.0000001]] -2.0 2.0 [-2, 2] [3.14549722 1.85450278] [[3.1454971, 3.1454972999999997], [1.8545026999999998, 1.8545029]] 1.8545027756320966 3.145497224367904 [1.8545027756320966, 3.145497224367904] [-1.41421356+1.41421356j -1.41421356-1.41421356j 1.41421356+1.41421356j 1.41421356-1.41421356j] [] [] [-0.+4.j 0.-4.j] [] [] [-1.e-06 1.e-06] [[-1.1e-06, -9e-07], [9e-07, 1.1e-06]] -1e-06 1e-06 [-1e-06, 1e-06] [9. 8. 7. 6. 5. 4. 3. 2. 1.] [[8.9999999, 9.0000001], [7.9999999, 8.0000001], [6.9999999, 7.0000001], [5.9999999, 6.0000001], [4.9999999, 5.0000001], [3.9999999, 4.0000001], [2.9999999, 3.0000001], [1.9999999, 2.0000001], [0.9999999, 1.0000001]] 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 [1, 2, 3, 4, 5, 6, 7, 8, 9]
In [ ]:
#check that there's a reasonable number of pre- and post-conditions inside FindRealZeros and plenty of commenting

Test your code from above here(1 point), ID: find_large_roots
In [ ]:
#here follow some hidden tests.
#very large roots

Test your code from above here(1 point), ID: root_finder_hidden_tests1
In [ ]:
#a few more general tests, not checking anything in particular, just keeping you honest

Test your code from above here(1 point), ID: root_finder_hidden_tests2
In [ ]:
#a few more general tests, not checking anything in particular, just keeping you honest

Test your code from above here(1 point), ID: root_finder_repeated
In [ ]:
#what happens when there are repeated roots?

Test your code from above here(1 point), ID: root_finder_uses_helpers
In [ ]:
#Ensure that the correct helper functions are being used

Test your code from above here(1 point), ID: root_finder_no_external
In [ ]:
#check that FindRealZeros does not rely on external functions


Explain how your function FindRealZeros works. Each member of the pair programming pair should write their own explanation, independently, in a separate copy of the file (which will contain the same solutions to the programming questions).

Marks: 3