Finding a Functions Roots with Python

Python is a fantastic programming language designed to be both intuitive to new programmers and flexible in terms of capability. But when you want to do something that on the surface seems simple like finding the roots of a function with python it’s much more complex than you may expect.

In this tutorial we are going to be writing a Python program that will attempt to find the roots of any function without importing in any modules. If you are looking for the easy way out SciPy has a great root finding library you can read the documentation to here.

Although you may think that python would have some way of isolating x in math functions and solving just like they taught you to in middle school. You’d unfortunately be wrong. Isolating variables is actually a fairly hard thing to do with a computer because there are too many factoring strategies Such as trig substitutions like tan(x)=sin(x)/cos(x) or the double angle identities or other factoring strategies like adding and subtracting something that it would take forever for us to program in.. That’s why there are so few competitors in the online equation solving space, wolfram alpha and symbolab being the two most recognizable.

Our solution will be much simpler than trying to isolate x and set it equal to zero. Instead we will be using the Newton-Raphson method to find roots. If you’ve never heard of this method before that is completely ok, I will be summarizing it below. I’d like to start off by pointing out that this method is not an exact science, and can both miss roots and create false positives. we will also only be checking for roots over a finitely large domain.(so we could miss infinitely many haha…)

The Newton Raphson method as pictured above is just a loop that follows this general pattern:

· We should start with a guess of where we think our root is going to be.

· We check if our guess is accurate by checking if it is very close to zero

· If it is, then we have found a solution and return that

· If not, then we will calculate the tangent line of our function at that point (the slope of one the dotted lines in the image above)

· We then calculate where the tangent line intersects with our x-axis. Then we use this as the x-value for our “guess” and repeat the process from the first bullet point.

After a few iterations we should see a solution if not we just give up.

This is called an iterative root finding algorithm. Although you will be able to choose how close you want your “guess” to be from the root. you generally won’t get an exact solution.

To make this algorithm more useful to us, we will be putting it in a loop to iterate over thousands of starting “guess” values. To hopefully find all of our function’s roots.

We will begin by writing a python function to perform the Newton-Raphson algorithm. Before we start, we want to decide what parameters we want our function to take in.

1. f: the math function
2. x0: our guess for the root’s location
3. epsilon: The distance from the x axis that we are willing to call our margin of error
4. Max_iterations: the maximum number of iterations that we want to attempt before aborting hopeless attempts such as if we hit a vertical asymptote.

We can now start writing our function

we are going to start off by setting xn=x0 as x0 represents our first guess, we will then write a for loop that will iterate max_iter many times:

We are now going to solve for our y value, if its absolute distance is within epsilon of zero, then we know we have found our solution and will return the solution.

if it is not our solution, we want to calculate the equation of the dotted line. We will start with the slope of the dotted line, it will be the derivative of the function at that point which we can approximate with (f(x+h)-f(x))/h using a very small h.

This function takes in 2 parameters, f: the math function, x: the x coordinate

Now we want to calculate the equation of our line which is given by:

You can play around with this equation here by dragging the x_n slider or changing the function to whatever you want.

we want to isolate x to find our point of intersection with the x axis.

Our code should now look as follows:

If the slope is zero, we want to return none as we have hit a peak/valley (local min/max) that is not a solution, so we abort this attempt.

We just have one last thing to do, outside the for loop we are going to return None. This is just in case a solution is never found it will return None.

We can now test our code out:

Note: By varying our x0 value that represents our “guess”, we are able to find different solutions to our equation. Also note that the smaller you make epsilon the more accurate your guess becomes as follows:

Our solution has now become so accurate that python rounds it to the correct solution. Also note that max_iter variable serves essentially no purpose here as our function is very well behaved very few iterations are needed to find a very accurate solution.

I’d also like to take a moment to point out that the algorithm we have just written is almost identical to how many scientific calculators approach the problem. Our algorithm just like my Casio fx-991EX requires an initial guess after you input the equation, with epsilon and max_iter simply being set to default values by the calculator for simplicity.

Isn’t it crazy how we can emulate one of Casios best features on their scientific calculators with just 14 lines of python?

To expand this solution to one that will find all roots and doesn’t require guesses we are simply going to write a loop that will repeatedly call our solver function.

I’m going to call this new function loop, it will take in 3 parameters, the left bound, right bound, increment. The left, right bound representing the interval we want to check for solutions, this is one of the greatest limitations of any root finding algorithm that don’t feature a CAS (computer algebra system) like fancy graphing calculators, wolfram alpha and symbolab do. Because we aren’t using a CAS we have no idea what the function looks or behaves like, we are only able to check for roots in a finite domain.

Now I will call our solver function

Note how I’ve hardcoded in the epsilon value and max_iter. These are values I choose randomly feel free to make epsilon larger and/or max_iter smaller if you need to increase performance (e.g if you had a very large domain).

If our solution is not None we are going to round it: we do this to reduce the number of “duplicate” solutions we are going to see. For example, if we had a single solution of zero. We wouldn’t want to see a dozen results incredibly close to zero e.g 0.00001, 0.00002, 0.0003,-0.00001.

Now we are going to check if our solution is in an array, we are going to call solutions if it isn’t we are going to add it in.

below the loop we are going to add in some print statements to visualize the fruit of our labour:

We can now test:

I’d like to point out that many of the solutions found are outside of our provided domain, this is because the domain we provide are used as guesses for the program but does not determine what parts of the function it will look for a solution.

We still have one function left to write; we are going to write f(x) to prevent zero division:

All this function does is it takes an x value attempts to plug it into the equation, it returns the result if it exists, if it has a zero-division error then it returns a number very close to zero.

Now we can test it out!

Our code will throw false positives if the function gets within epsilon of the x axis but never “touch’s”. e.g

This actually has no solutions as the function is undefined at x=0 as it becomes division by zero.

Our code also will generally run slowly for equations with no solutions.

I’d also like to briefly touch on the accuracy of our function. I can only really guarantee good accuracy when the function is continuous. Although in my personal testing it works fine with discontinuous functions like tan(x). if you put in a function with discontinuities on either side of a solution you may miss solutions. You may also miss solutions if your solutions are “denser” than your “increment” I use 0.5 in my examples. but if your function has more than 2 roots within an interval 0.5 wide you will likely miss roots and should decrease incrementation size.

I would also like to mention that this can also be used to solve equations, eg the intersection of two quadratics.

For example to find the points of intersection between x²-4x-2 and 3x²+3x+1 you could move all like terms to one side.

You can verify this with Desmos

Here’s a copy of the final code:

import math
def derivative(f, x):
 h=1e-8
 return (f(x+h)-f(x))/h
def solver(f, x0, epsilon, max_iter):
 xn=x0
 for n in range(0,max_iter):
 y=f(xn)
 if abs(y)<epsilon:
 return xn
 slope=derivative(f,xn)
 if(slope==0):
 return None
 xn=xn-y/slope
 return None
def loop(f, L_bound, R_bound, increment):
 solutions=[]
 while L_bound<=R_bound:
 solution=solver(f, L_bound, 1e-10, 1000)
 if solution is not None:
 solution=round(solution,4)
 if solution not in solutions:
 solutions.append(solution)
 L_bound+=increment
 print(sorted(solutions))
 print(“we found “+str(len(solutions))+” solutions!”)
equation=””
def f(x):
 try:
 y=eval(equation)
 except ZeroDivisionError:
 y= 1e-10
 return y
equation=”x**2–9"
loop(f,-100,100,0.5)