CoCalc Public Filesweek 5 / assignment / FindRoots.ipynbOpen with one click!
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
Answer that will be automatically graded below, ID: cell-73de33f8730ab565
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

Checking Your Understanding

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)=x312x2+45x52p(x)=x^3-12x^2+45x-52. How do you express p(x)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)p'(x) (as a list called polyDeriv) and the root(s) of p(x)p'(x) (as a list, in increasing order, called derivRoots)?
  3. What is the Lagrange bound as applied to the polynomial p(x)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)p(x), but you know the roots of p(x)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) (-109,3),(3,5),(5,109) You should do this without evaluating p(x)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)p(x)? Give your answer as a list in increasing order called polyRoots. Make sure that all the roots of p(x)p(x) appear in exactly one of the above ranges.
Marks: 5
Answer that will be automatically graded below, ID: cell-cdaeaff6f2fb9c05
In [ ]:
import math # YOUR CODE HERE 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) (-109,3),(3,5),(5,109) You should do this without evaluating p(x)p(x) at any point.

Marks: 2

Manually graded answer(2 points), ID: explain_root_intervals

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.

Readonly, ID: cell-0c08ceaeaa375aea
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[0]*i[1]<0 for i in zip(PosList[1:],PosList[:-1])]),sum([i[0]*i[1]<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
Answer that will be automatically graded below, ID: cell-cac064c0ffad704e
In [10]:
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 ''' # YOUR CODE HERE 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 emptyList.append(listAdd) 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 emptyList.append(listAdd) 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

Answer that will be automatically graded below, ID: cell-6260876b18d911a3
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, nn, what's the maximum number of ranges there could be? How many real roots could there be for a degree nn 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 ''' # YOUR CODE HERE 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 emptyList.append(listAdd) 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 emptyList.append(listAdd) 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[0]*i[1]<0 for i in zip(PosList[1:],PosList[:-1])]),sum([i[0]*i[1]<0 for i in zip(NegList[1:],NegList[:-1])])]) DescartesSigns([-52,45,-12,1])
Answer that will be automatically graded below, ID: cell-8808af59d41588bd
In [58]:
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''' # YOUR CODE HERE 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[0],x[1]) #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 [59]:
#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[0]-i[1])<=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]
Test your code from above here(1 point), ID: root_finder_comments_asserts
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
Manually graded answer(3 points), ID: root_finder_explanation

The function FindRealZeros has been created to find the zeros of a quadratic equation through repeated iterations of Rolle's theorem using the 'FindZeroInInterval' function. We can then locate the roots by narrowing the interval in which the root could possibly be in, using the root_bounds function. Also, by using Descartes' rule of signs and the 'RootLocationBound' we are able to find the upper bounds of both the positive and negative roots, which then provides an interval within which all of the roots of the given polynomial will exist. The actual FindRealZeros function is used to reorganise the roots found using multiple different functions and return them in increasing order.

Manually graded task(5 points), ID: code_style

General feedback on your code structure, style, commenting, tests, explanations etc. You do not have to do anything here.

Marks: 5
In [ ]: