Solving equation using bisection method

Question

Is there a bisection method I can find online, specifically for python?

For example, given these equations how can I solve them using the bisection method?

x^3 = 9  
3 * x^3 + x^2 = x + 5  
cos^2x + 6 = x  

Strand

I wish my numerical methods course had used Python. :/ This is really instructive to implement yourself; just read Wikipedia's description for the algorithm. – Metrorrhagia

Best to use something that's already been in use by many people than to try to write it yourself. 75%-90% of binary search implementations are incorrect. – Audiology

Answer 1

Using scipy.optimize.bisect:

import scipy.optimize as optimize
import numpy as np

def func(x):
    return np.cos(x)**2 + 6 - x

# 0<=cos(x)**2<=1, so the root has to be between x=6 and x=7
print(optimize.bisect(func, 6, 7))
# 6.77609231632

optimize.bisect calls _zeros._bisect, which is implemented in C.

Shivery

how to get the a and b value? and also the number of loop? – Strand

like in the example math.fullerton.edu/mathews/n2003/BisectionMod.html it is trying to solve this equation (x^3+4x^2-10=0) and it uses function Bisection[1,2,30] how to get the number 1 2 and 30? is it from the equation? – Strand

@bn: To use bisect, you must supply a and b such that func(a) and func(b) have opposite signs, thus guaranteeing that there is a root in [a,b] since func is required to be continuous. You could try to guess the values for a and b, use a bit of analysis, or if you want to do it programmatically, you could devise some method of generating candidate a and b until you find two that have opposite signs. That's all beyond the scope of the simple bisect method however. – Shivery

@bn: Regarding the number of iterations, scipy.optimize.bisect has a maxiter argument. For example, so.bisect(func,6,8,maxiter=500). See docs.scipy.org/doc/scipy/reference/generated/…. – Shivery

@unutbu, thanks! for (equation x^3 = 9 ) i tried value of 0 for a and 1000 for b (2.10571289062) . and it gives me different answer with 0 for a and 5 for b (2.0802307129) according to wolfram, this 2.0802307129 is the answer wolframalpha.com/input/?i=x^3+%3D+9++ – Strand

@bn: For a odd-order polynomial like f(x) = x^3 + 4x^2 - 10, the x^3 term is going to dominate when x is large and negative, and when x is large and positive. When x is large and negative, f(x) is going to be negative, and similarly, when x is large and positive, f(x) is going to be positive. So it is easy to guess a and b in this case: Just choose a large negative value (like -100) for a and a large positive value (like 100) for b. As long as the number you choose is larger (in magnitude) than the coefficients, it should work. – Shivery

@bn: The bisect method quits when the root is known to lie within xtol of the returned value. You can try to make the answer more accurate by making the xtol parameter smaller. For example, bisect(func,0,5,xtol=1e-15). – Shivery

@bn: Correction: the number you choose should be significantly larger than the coefficients. Being just slightly larger than the largest coefficient might not cut it. – Shivery

@bn: If you're looking for the a and b values just plot the function over the range you are interested in. This is especially important if the function has a lot of zeros. – Baalman

Please remove the link to the execrable implementation in Python. Thanks. – Animism

Answer 2

This could help you!

import numpy as np

def fn(x):
    # This the equation to find the root
    return (x**3 - x - 1) #x**2 - x - 1

def find_root_interval():
    for x in range(0, 1000):
        if fn(x) < 0:
            lower_interval = x
            if fn(x+1) > 0:
                higher_interval = x + 1
                return lower_interval, higher_interval
    return False

def bisection():
    a,b = find_root_interval()
    print("Interval: [{},{}]".format(a,b))
    # Create a 1000 equally spaced values between interval
    mid = 0
    while True:
        prev_mid = mid
        mid = (a+b)/2
        print("Mid value: "+str(mid))
        # 0.0005 is set as the error range
        if abs(mid-prev_mid) < 0.0005:
            return mid
        elif fn(mid) > 0:
            b = mid
        else:
            a = mid

root = bisection()
print(root)

Marcelo