Over the past few days I was debating what I should write about next. I really enjoy sharing my thoughts and have been getting great responses from my readers. I thanks you for your comments. At your request I have decided to cover another math topic, solving equations using python. Equation solving is something fundamental in computing, mathematics and other scientific applications. Using computers enables users to solve a wider range of equations than we could by non-computer methods. However, there is a cost at hand, accuracy. Although computers can solve equations using a few methods, we will witness that accuracy will change between methods. In some cases you might not get to a solution at all.

Before we dive into code behind solving equations, let us step back. A good programmer should always consult other resources to see what else has been done already. The Internet (and Google) have given us great tools to research any topic our heart desires. There are many tools available to solve equations and do almost any other mathematical calculations. Matlab and Matamatica are both great mathematical softwares that can do amazing things. However, they are paid distributions. I do not have anything against paid software, but for the sake of educational useage I would like to keep everything in the realm of free (open source preferred). I like to post code that anyone can download, run, modify, distribute or do anything else they would like to do with it. Using paid languages will be counter productive to such goals and limit the readers of this blog. That being said, if there is a request I will put up Matlab code for the algorithms discussed in this post.

On another note, we will be using Python2.6 for all the examples that are to follow. We will not be using NumPy or SciPy. Both are great scientific math tools but an over kill for what we are trying to do. Our mission is simple, solving an equations. More precisely, we want to figure out where a given function intersect the x-axis without using “the big guns”, including Wolfram-Alpha. So let us get to it.

The Bisection Method:

At out first stop, the most fundamental approach to solving equations is the Bisection Method. Here is the theory in a nutshell:

“If f is a continues function on [a,b], where f(a)f(b) < 0, there must be a point r between a and b where f(r) = 0”

So what does this means? It means that if we have a function that cross the axises between point a and b, then the value of the function at f(a), multiplied by the value of the function at b, f(b), must be smaller that 0. Further more between a and b there is a point where the function f is 0. Here is a picture describing such event:

Here is a general solution in code that does exactly that:

def bisection(a,b,tol):
	c = (a+b)/2.0
	while (b-a)/2.0 > tol:
		if f(c) == 0:
			return c
		elif f(a)*f(c) < 0:
			b = c
		else :
			a = c
		c = (a+b)/2.0
		
	return c

Notice that we have a while loop that runs while we are within tolerance and that we assume function f(x). The latter part can be solved by writing a function f that outputs the desired function. Here is the function we are going to work with:

And here is the python equivalent:

def f(x):
	return x**3 + x -1

Once we add the code together we can run it and get a solution within tolerance. What is the tolerance? The Bisection method literally cuts in half the solution at each iteration, assuming the mid point between a and b is the solution. So in theory we can keep asking for a more accurate solution by setting a lower tolerance. The “real” solution to the function according to Wolfram Alpha is 0.6823278... We will use this solution as a guide to compare the accuracy of solutions we compute. Here are a few outputs for our code, for a solution between 0 and 1 using various tolerances:

(The full code for the Bisection Methon in Python)

As you can see, the first solution is only accurate to 2 decimal places. If we tried a lower tolerance we would have got maybe one or none of the digits correct. When we increased our tolerance we arrive to a more accurate solution, within tolerance range. Bisection is a good method to use when it applies. The Bisection method will not work on all functions, but is a quick solution to many equations.

Fixed Point Iteration (FPI)

Fixed point iteration is another method of solving equations. The definition of a fixed point r is when g(r)=r. Using this definition, we start with some initial guess, x and plug it into the function. However, unlike bisection, this method does require us to do some math work before code. We have to transform our function to be in the form x=g(x). This is where the heart of FPI is. Based on your setup and initial guess, you may reach a solution really fast, or not at all. Here is how the code for FPI looks like:

def fpi(x, k):
	for i in range(k):
		x = g(x)
	return x

You will note that now we are no longer using a while loop. FPI could be setup in a smilier way to Bisection, but for simplicity is setup in a fixed i-loop. Using the same function as before, let us say we transform the function to be like this:

And in python (using python’s math library):

def g(x):
	return math.pow(1-x,1/3.0)

You will notice that python only has a square root function. For any other roots we use python’s power function. Please remember that 1/3 in integer division is 0. Therefore, we have to use 1/3.0 to get the cubic root. Running the program for 10 iterations starting at 0.5 yields the following output:

(The full code for Fixed Point Iteration in Python)

This is not very close, we only have 1 digit and we started close to the solution. If we let the program run for 100 iterations we get really close to the solution:

How many iteration do you need? It depends on your function. A variation on this function might solve it in 7 iterations (can you find the function?). Again, in FPI it depends on the setup and the initial guess. As a rule of thumb, it will either conversation very fast, or not at all. But then again, sometimes it just needs one or more steps.

Newton’s Method

Our next stop, and the basis for our next method, is Newton’s method. In this case, like FPI we start with an initial guess and keep redefining x, however, we use this function to find the next x :

Newton’s method adds another math task we will have to do outside of our code, taking the derivative. We will see that the last method will take care of this, but at another cost. Here is how Newton’s method will look like in code:

def newt(x,n):
	for i in range(n):
		x = x - f(x)/f_prime(x)
	return x

This means we will have to define our original function as well as a derivative (prime) function. Not a problem, we can use the same function we have been solving all along. This is what it will look like:

And in Python:

def f(x):
	return x**3 + x - 1

def f_prime(x):
	return 3*x**2 + 1

Before we run our code, we need to take care of one special case. What if the derivative is 0 or end up computing to 0? Than we will have a division by zero error. This can be simply avoided by a couple of if statements. Like so:

def newt(x,n):
	for i in range(n):
		if f_prime(x) == 0:
			return x
		x = x - f(x)/f_prime(x)
	return x

This should take care of any division by zero errors. However, do not take my word for it, try it out for yourself!

And now running the code, using -0.7 as our initial x for 10 iterations will yield:

(The full code for Newton’s Method in Python)

If you add some debug statements and examine the output you will note that the solutions was found after only 7 iterations. That is really fast and a lot better than 100 iterations of FPI. You might have output like this:

Secant Method

Although Newton’s method works really well, taking derivatives is not always an option. For that, there is another method, mathematically based on Newton’s method and the definition of derivatives. This time we start off with 2 initial guesses and move forward by using one of the guess and a new point. The new point is calculated from the formula:

In Python the Secant Method might look like this:

def secant(x0,x1,n):
	for i in range(n):
		x_temp = x1 - (f(x1)*(x1-x0)*1.0)/(f(x1)-f(x0))
		x0 = x1
		x1 = x_temp
	return x1

Note that there is a chance that we will be dividing by 0, again we must add a check to ensure that this does not happen. Also note that we use a temp variable to swap values around. This is so we do not require much memory. A new, modified version of the Secant method might look like this:

def secant_2(x0,x1,n):
	for i in range(n):
		if f(x1)-f(x0) == 0:
			return x1
		x_temp = x1 - (f(x1)*(x1-x0)*1.0)/(f(x1)-f(x0))
		x0 = x1
		x1 = x_temp
	return x1

As before, we are using for loops, which might be swapped to while loops. I recommend you try it your self, but remember, while loops must terminate at some point. Running the Secant code starting with 0 and 1 for 10 iterations, on out “famous” function yields the following output:

(The full code for Secant Method in Python)

We have gone over 4 basic methods of solving equations. There are many more approaches to solve equations, this is by no means a complete survey. I trust that if you are in need for more information you can find it using a friendly search engine. Note that there is still much to say about the accuracy of the solutions computed by these methods. Also note that there are some programming “tricks” that can be done to shorten calculation. For more accurate mathematical definitions and information please seek out a mathematician.

This post was made possible by Numerical Analysis by Timothy Sauer