Bisection Method
Moderator: Board moderators

 Experienced poster
 Posts: 202
 Joined: Fri Mar 22, 2002 2:00 am
 Location: Chittagong. CSE  CUET
 Contact:
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.0e10
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
[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.0e10
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
We are all in a circular way, no advances, only moving and moving!

 Experienced poster
 Posts: 202
 Joined: Fri Mar 22, 2002 2:00 am
 Location: Chittagong. CSE  CUET
 Contact:
Would you please explain me what you really looking for.......?Larry wrote:Can someone explain either or both methods? Thanks! =)
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 ???
We are all in a circular way, no advances, only moving and moving!

 Experienced poster
 Posts: 202
 Joined: Fri Mar 22, 2002 2:00 am
 Location: Chittagong. CSE  CUET
 Contact:
Before posting the C++ code I would like to make you read it ..............
http://www.eit.ihkedu.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 ???
http://www.eit.ihkedu.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 ???
We are all in a circular way, no advances, only moving and moving!
HI!
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.)???
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
Date: December 31st, 2011 (Saturday)
Time: 12:00  16:00 (UTC)
URL: http://uva.onlinejudge.org
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 (n1). 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 + an1 xn1 + ....+ a1x + a0 = (x
p (x) = (x  xr) q (x) can be written as
an xn + an1 xn1 + ....+ a1x + a0 = (x

 New poster
 Posts: 20
 Joined: Thu Apr 10, 2003 10:53 pm
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 n1 roots. So I just take the derivative of f(x), and recursively find the roots for it. The recursion will stop after the (n1)th derivative which is just linear. Anyway, say the roots of f'(x) are r1, r2, r3, ... r_{n1} in increasing order. Then I do a binary search in each the intervals (inf,r1), (r1,r2), ..., (r_{n1},+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 (xr) 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.
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 (xr) 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.
Last edited by Lars Hellsten on Tue Jul 01, 2003 4:09 am, edited 2 times in total.

 Guru
 Posts: 647
 Joined: Wed Jun 26, 2002 10:12 pm
 Location: Hong Kong and New York City
 Contact:
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?
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?

 New poster
 Posts: 20
 Joined: Thu Apr 10, 2003 10:53 pm
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.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?
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( hilo > 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.
(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 n1].
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
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 n1].
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
"Everything should be made simple, but not always simpler"
After so long, I finally manage to solve Q10428 The Roots with bisection method. Those learning bisection should really try this problem out!
7th Contest of Newbies
Date: December 31st, 2011 (Saturday)
Time: 12:00  16:00 (UTC)
URL: http://uva.onlinejudge.org
Date: December 31st, 2011 (Saturday)
Time: 12:00  16:00 (UTC)
URL: http://uva.onlinejudge.org
Observer wrote:Easy!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???????
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 "NR 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?
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.
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
Date: December 31st, 2011 (Saturday)
Time: 12:00  16:00 (UTC)
URL: http://uva.onlinejudge.org