Bisection Method

[Moni]

Look at the code .............


[cpp]
// bisect.cpp
// Program to find a root of an equation by the bisection method
// A J Norman, Leicester University Engineering Department, Dec 1993
// Pascal -> C++ Jul 1995

#include <iostream.h>
#include <iomanip.h>
#include <math.h>
#include <string.h>
#include <time.h>

#define MIN_TOLERANCE 1.0e-10

int bisection (double, double, double, int, int &,
double &, double &, char []);
void get_conditions (double &, double &, double &, int &);
double user_function (double x);

int main (void)
{
clock_t start, finish;
double x1, x2, tolerance, root, value, time_taken;
int max_iterations, iterations;
char error, err_string[80];

get_conditions (x1, x2, tolerance, max_iterations);

start = clock();

// Do the calculation 50000 times so that the time taken is a few
// seconds, rather than a few milliseconds.
for (unsigned long i = 0 ; i < 50000; i++)
{
error = bisection(x1, x2, tolerance, max_iterations, iterations,
root, value, err_string);
}

finish = clock();
time_taken = (finish - start) / double(CLOCKS_PER_SEC);

// Print out the number of iterations, or an error message

if (error)
{
cout << "ERROR : " << err_string << endl;
}
else
{
cout << "Successful calculation, " << iterations
<< " iterations." << endl
<< "The root is " << root << endl
<< "The value of the function at this point is "
<< setiosflags(ios::fixed) << setprecision(2)
<< value << endl
<< "Time for 1000 bisection calculations was "
<< time_taken << " seconds." << endl;
}

return 0;
}

//======================================================================

void get_conditions (double & x1, double & x2, double & tolerance,
int & max_iterations)
{
cout << "Bisection method to find the root of an equation." << endl
<< "Input two values of x." << endl
<< "f(x1) must be positive and f(x2) negative, or vice versa."
<< endl << "X1 = ";
cin >> x1;
cout << "X2 = ";
cin >> x2;

cout << "Calculation will terminate when f(x) is less than the"
<< endl << "tolerance level input by the user." << endl;

do
{
cout << "Tolerance : ";
cin >> tolerance;
if (tolerance <= MIN_TOLERANCE)
{
cout << "ERROR : tolerance must be greater than "
<< MIN_TOLERANCE << endl;
}
}
while (tolerance <= MIN_TOLERANCE);

do
{
cout << "Maximum number of iterations (1 or more) : ";
cin >> max_iterations;
}
while (max_iterations < 1);
}

//======================================================================
// Function to find a root of an equation by the bisection method.
// x1, x2 : Range of x within which the root lies. The user
// must ensure that the values of the function at x1
// and x2 have opposing signs.
// tolerance : tolerance level of the answer
// max_iterations : maximum number of iterations to be performed
// root : final value of the root
// value : value of user_funct at the root
// iterations : number of iterations performed
// err_string : string containing an error message, if appropriate
//
// Function returns 1 (true) if an error has been encountered.
//======================================================================

int bisection (double x1, double x2, double tolerance,
int max_iterations, int & iterations,
double & root, double & value, char err_string[])
{
double y1, y2, y3, a1, a2, a3;

// Error checking, to make sure that the range x1 -> x2
// really does contain a root.

y1 = user_function(x1);
y2 = user_function(x2);

if (y1 < 0.0)
{
if (y2 < 0.0)
{
strcpy(err_string,
"The function is negative at both ends of the range");
return 1;
}
else
{
a1 = x1;
a2 = x2;
}
}
else
{
if (y2 > 0.0)
{
strcpy(err_string,
"The function is positive at both ends of the range");
return 1;
}
else
{
a1 = x2;
a2 = x1;
}
}

// Begin the iterations.

for (iterations = 1; iterations <= max_iterations; iterations++)
{
a3 = (a2 - a1) / 2.0 + a1;
y3 = user_function(a3);

if (y3 < 0.0) { a1 = a3; } else { a2 = a3; }

if (fabs(y3) <= tolerance)
{
root = a3;
value = y3;
return 0;
}
}

strcpy(err_string, "Maximum number of iterations was exceeded");
return 1;
}

//======================================================================
// This function describes the mathematical function whose root is
// to be found.
//======================================================================

double user_function (double x)
{
return x*x - 5.0*x + 6.0;
}


[/cpp]


I think is it clear to you now

[Larry]

Can someone explain either or both methods? Thanks! =)

[Moni]

Can someone explain either or both methods? Thanks! =)

Would you please explain me what you really looking for.......?

When f(x) is continuous in the environment of a single root r, then f(x) changes sign through r. There is a small interval [a,b] including r such that f(a).f(b) < 0. Taking the midpoint m of [a,b], there are three possibilities.
f(m) = 0 ; then m is the root r.
f(m).f(a) < 0 ; then the root r is in [a,m]
f(m).f(b) < 0 ; then the root r is in [m,b]
Now we can restart the procedure with the smaller interval [a,m] or [m,b]. The interval becomes smaller and smaller, so we can find an approximation of the root r.

In the above the code is given anything else ???

[anupam]

well moni bhai, do u have the class of nr method. or if any1 got this please post it here.
share it with all of us.

[Moni]

Before posting the C++ code I would like to make you read it ..............

http://www.eit.ihk-edu.dk/subjects/cpp/newton.html

And I have an idea...........

One can post his code and we all can develop and find mistake or do improvement on that like Linux kernel programming runs............sounds how ???

[Observer]

So, using NR method, we should pay extra attention to the values of x where f'(x) = 0. Is that correct??

Should I do the checking this way?
1. Calculate f'(x).
2. If f'(x) = 0, then
3. - if f(x) = 0, then it's a root
4. Go to the next value x <---- How? By substraction/addition by an arbitarily small value???

I'll try to use this method to solve Q10428(The roots)......
By the way, do I need sorting in this qq?? And do I have to solve f'(x) = 0 (and thus f''(x)=0 etc.)???

[IIUC GOLD]

A polynomial function of degree n can be expressed as p (x) = (x - xr) q (x) where xr is the root of the polynomial p(x) and q(x) is the quotient polynomial of degree (n-1). By synthetic division we can obtain q(x) without performing the actual division. Synthetic division is performed as follows:

p (x) = (x - xr) q (x) can be written as
an xn + an-1 xn-1 + ....+ a1x + a0 = (x

[Lars Hellsten]

I used the bisection method (binary search) to solve 10428. You just need an extra trick to ensure you only consider monotonic intervals: notice that every root of the polynomial f(x) must lie between some local minimum and local maximum. These local optima are precisely the points where the function f'(x) = 0. If f(x) has n roots then f'(x) has n-1 roots. So I just take the derivative of f(x), and recursively find the roots for it. The recursion will stop after the (n-1)th derivative which is just linear. Anyway, say the roots of f'(x) are r1, r2, r3, ... r_{n-1} in increasing order. Then I do a binary search in each the intervals (-inf,r1), (r1,r2), ..., (r_{n-1},+inf) to find the roots for f(x).

Anyway, I thought this method was kind of interesting since it's easy to come up with from scratch if you don't know NR. I am embarassed to admit that I did know about NR and only did it this way because I wasn't sure how to proceed after finding the first root -- I forgot I could just divide out (x-r) once I found a root r.

The NR way or secant methods are more efficient at converging, and handle roots with multiplicity > 1 or polynomials with fewer than n real roots more elegantly (not an issue for 10428). The bisection approach is easier to understand and implement.

[Larry]

Thanks Moni, I actually didn't look at the code, because I want to learn concepts, but it seems like to me it's kind of like binary search..

How do you find the intervals, Lars? Ya, I notice that the abs of roots are within 25, but how do you choose the intervals so that you are certain there's not 2 or more values between that interval?

[Moni]

The code I posted is really very educational & helpful for me.......I have learnt from it

Lars I really want to know the trick used by great SnapDragon(Kisman).

Can you explain it and share with us ???

[Lars Hellsten]

How do you find the intervals, Lars? Ya, I notice that the abs of roots are within 25, but how do you choose the intervals so that you are certain there's not 2 or more values between that interval?

There can't be 2 or more roots in an interval, because the intervals are chosen so that the function is monotonic on each of them, by a recursive call. Imagine a polynomial. It has a bunch of peaks and valleys. It shouldn't be hard to see that the interval between a peak and a neighbouring valley can only have one root. The roots of the first derivative correspond to the peaks and valleys. It's possible for an artibrary polynomial that there are no roots in an interval, but there can't be two or more. In this question, it is given that there are n distinct real roots, which means that each of the intervals does have one.

Moni: the "trick" isn't really a trick, just something I hadn't seen much before, a ternary search. It's is similar to binary search, but for optimizing a continuous function that doesn't have local optima aside from the one global optima. e.g.:

[c]
while( hi-lo > epsilon ) {
a = (2*lo+hi)/3;
b = (lo+2*hi)/3;
if( f(a) < f(b) ) lo = a;
else hi = b;
}
[/c]

Not needed here though since binary search is fine.

[anupam]

(to iiuc gold) we know the method Synthetic division as nested multiplication
to find all roots of a given eqn by nr method.
Yes, it's very easy because after finding out a root we can easily
obtain another polynomial of degree - 1 by replacing P(x)[degree n]
by Q(x)[degree n-1].
And i haven't try the problem The root yet.I can tell about the
problem more precisely after solving it.

(forgive my poor english)
--
forgive me for my preveous fault
Anupam

[Observer]

After so long, I finally manage to solve Q10428 The Roots with bisection method. Those learning bisection should really try this problem out!

[junbin]

In bisection method for a & b f(a) & f(b) should be of opposite sign.But how to determine a & b while implementing it in C???????

Easy!

For instance, in Q358, first you should derive a formula f(x) = 0.
Then, you can clearly see that the function f(x) is strictly increasing.
Now choose the upper and lower limit. You may set 0 as the lower limit, and 2R as the upper limit.
Next, do the bisection!!
Finally, you get the answer easily!!!!!

By the way, what's "N-R Method"??

Railway is quite interesting! I didn't use bisection, though...

Am I right to say that bisection method means binary search?

ie: You "guess" an answer, determine if the correct answer should be higher or lower, then guess a number that is closer and you keep guessing until you guessed the correct answer?

[Observer]

You may call it "binary search" if you wish.

However, we usually use the term "binary search" when we're doing searching in a sorted array, and "bisection method" when we're doing the "guessing" as you've mentioned.