Bisection Method

Let's talk about algorithms!

Moderator: Board moderators

Moni
Experienced poster
Posts: 202
Joined: Fri Mar 22, 2002 2:00 am
Location: Chittagong. CSE - CUET
Contact:

Post by 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 :)
ImageWe are all in a circular way, no advances, only moving and moving!

Larry
Guru
Posts: 647
Joined: Wed Jun 26, 2002 10:12 pm
Location: Hong Kong and New York City
Contact:

Post by Larry »

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

Moni
Experienced poster
Posts: 202
Joined: Fri Mar 22, 2002 2:00 am
Location: Chittagong. CSE - CUET
Contact:

Post by Moni »

Larry wrote: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 ???
ImageWe are all in a circular way, no advances, only moving and moving!

anupam
A great helper
Posts: 405
Joined: Wed Aug 28, 2002 6:45 pm
Contact:

Post by 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. :oops: :oops:
"Everything should be made simple, but not always simpler"

Moni
Experienced poster
Posts: 202
Joined: Fri Mar 22, 2002 2:00 am
Location: Chittagong. CSE - CUET
Contact:

Post by 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 ??? :)
ImageWe are all in a circular way, no advances, only moving and moving!

Observer
Guru
Posts: 570
Joined: Sat May 10, 2003 4:20 am
Location: Hong Kong

HI!

Post by 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.)???
7th Contest of Newbies
Date: December 31st, 2011 (Saturday)
Time: 12:00 - 16:00 (UTC)
URL: http://uva.onlinejudge.org

IIUC GOLD
New poster
Posts: 19
Joined: Tue Jun 11, 2002 4:27 pm
Location: Bangladesh
Contact:

Post by 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
New poster
Posts: 20
Joined: Thu Apr 10, 2003 10:53 pm

Post by 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. :oops:

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.
Last edited by Lars Hellsten on Tue Jul 01, 2003 4:09 am, edited 2 times in total.

Larry
Guru
Posts: 647
Joined: Wed Jun 26, 2002 10:12 pm
Location: Hong Kong and New York City
Contact:

Post by 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
Experienced poster
Posts: 202
Joined: Fri Mar 22, 2002 2:00 am
Location: Chittagong. CSE - CUET
Contact:

Post by 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 ??? :)


ImageWe are all in a circular way, no advances, only moving and moving!

Lars Hellsten
New poster
Posts: 20
Joined: Thu Apr 10, 2003 10:53 pm

Post by Lars Hellsten »

Larry wrote: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
A great helper
Posts: 405
Joined: Wed Aug 28, 2002 6:45 pm
Contact:

Post by 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 :oops: :oops:
Anupam
"Everything should be made simple, but not always simpler"

Observer
Guru
Posts: 570
Joined: Sat May 10, 2003 4:20 am
Location: Hong Kong

Post by Observer »

After so long, I finally manage to solve Q10428 The Roots with bisection method. Those learning bisection should really try this problem out! :wink:
7th Contest of Newbies
Date: December 31st, 2011 (Saturday)
Time: 12:00 - 16:00 (UTC)
URL: http://uva.onlinejudge.org

junbin
Experienced poster
Posts: 174
Joined: Mon Dec 08, 2003 10:41 am

Post by junbin »

Observer wrote:
Faizur wrote: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
Guru
Posts: 570
Joined: Sat May 10, 2003 4:20 am
Location: Hong Kong

Post by 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.
7th Contest of Newbies
Date: December 31st, 2011 (Saturday)
Time: 12:00 - 16:00 (UTC)
URL: http://uva.onlinejudge.org

Post Reply

Return to “Algorithms”