/* bisection.c: Implements the bisection method for finding roots.
*
* By Terry R. McConnell
*
/*
Solves f(x) = 0, given an initial interval [a,b] such that f has
differing signs at the endpoints. Quits when f((a+b)/2) is close
enough to zero. Function f must be continuous. It must be implemented
as double f(double) in this file or elsewhere. We include an example
implementation with f(x) = 3x^3 - x - 1. The endpoints may be
specified with -a and -b options. They default to a = 0, b = 1.0.
If the endpoints don't have the required property the program quits
with an error.
The algorithm bisects the initial interval successively, until
the value of f at the midpoint differs from zero by less than
TOLERANCE, which defaults to .00000001. (It can be set on the
command line with -t option.) The final result is the midpoint
of the last interval. Which of the two bisected halves is retained
at a given stage is determined by calculating f at the endpoints,
so that the condition of different signs at the endpoints is
preserved.
The program reports the current interval, and the values of f at the
endpoints at each stage.
The algorithm is also implemented in a user callable function.
*/
/* compile: cc -o bisection.c -lm
Use -D_SHORT_STRINGS if your compiler does not support multiline
string constants.
Use -DHAVEF and link object file containing implementation of
function declared as double f(double) if you want to solve a non
default equation. In this case, also -DEQUATION= so the printout
will give the correct heading.
Use -DNOMAIN and -DHAVEF if you only want to link in function implemented
below with declaration
double bisection(double a, double b, double(*f)(double), double tolerance);
It is up to the caller to make sure the endpoints and f are suitable.
You're toast if they're not.
*/
#include
#include
#include
#ifndef EQUATION
#define EQUATION "x^2 + 2x + 1 = 0"
#endif
#define VERSION "1.0"
#define USAGE "bisection [ -a float -b float -t float -h -v]"
#ifndef _SHORT_STRINGS
#define HELP "\nbisection [ -a float -b float -t float -h -v ]\n\n\
Find root of linked function f using bisection method.\n\n\
-a: Use next argument as left endpoint of initial interval.\n\
-b: Use next argument as right endpoint of initial interval.\n\
-t: Use next argument as tolerance. Quit when f(midpoint) dips below this.\n\
-v: Print version number and exit. \n\
-h: Print this helpful information. \n\n"
#else
#define HELP USAGE
#endif
/* Default values */
#define TOLERANCE .00000001
#define LEFT 0.0
#define RIGHT 1.0
#define MAXPASS 256
#ifndef HAVEF
/* Here is the default function f's implementation. */
double f(double x){
return x*x + 2*x + 1;
}
#else
extern double f(double);
#endif
#ifndef NOMAIN
int
main(int argc, char **argv)
{
double tolerance = TOLERANCE;
double a = LEFT;
double b = RIGHT;
double c,d,e,mid;
int j=0,n=1;
/* Process command line */
while(++j < argc){
if(argv[j][0] == '-')
switch(argv[j][1]){
case 'a':
case 'A':
if(j+1 >= argc){
fprintf(stderr,"%s\n",USAGE);
exit(1);
}
a = atof(argv[j+1]);
j++;
continue;
case 'b':
case 'B':
if(j+1 >= argc){
fprintf(stderr,"%s\n",USAGE);
exit(1);
}
b = atof(argv[j+1]);
j++;
continue;
case 't':
case 'T':
if(j+1 >= argc){
fprintf(stderr,"%s\n",USAGE);
exit(1);
}
tolerance = atof(argv[j+1]);
if(tolerance <= 0.0){
fprintf(stderr,"tolerance must be positive.\n");
exit(1);
}
j++;
continue;
case 'v':
case 'V':
printf("%s\n",VERSION);
exit(0);
case '?':
case 'h':
case 'H':
printf("%s\n",HELP);
exit(0);
default:
fprintf(stderr,"bisection: unkown option %s\n",
argv[j]);
exit(1);
}
}
/* Check necessary conditions */
if( b < a ){
fprintf(stderr,"Error: initial endpoints are swapped.\n");
exit(1);
}
c = f(a); d = f(b);
if(((c>0.0)&&(d>0.0))||((c<0.0)&&(d<0.0))){
fprintf(stderr,"f has same sign at enpoints. Cannot continue.\n");
exit(1);
}
/* loop until desired tolerance */
printf("\n\nBisection method solution of %s:\n\n",EQUATION);
while(1){
mid = (a+b)/2.0;
e = f(mid);
printf("%2d. f[%10.8f,%10.8f] = [%10.8f,%10.8f]\n",
n,a,b,c,d,e);
if(fabs(e) < tolerance){
printf("\nSolution is x = %10.8f at tolerance %9.8f.\n\n",mid,tolerance);
exit(0);
}
if( c > 0.0)
if(e < 0.0) b = mid;
else a = mid;
if( c < 0.0)
if(e > 0.0) b = mid;
else a = mid;
c = f(a); d = f(b);
n++;
if(n>MAXPASS) exit(1); /* obviously not converging */
}
return 0;
}
#endif
double bisection(double a, double b, double(*f)(double),double tolerance)
{
double mid,c,d,e;
c = f(a); d = f(b);
while(1){
mid = (a+b)/2.0;
e = f(mid);
if(fabs(e) < tolerance)break;
if( c > 0.0)
if(e < 0.0) b = mid;
else a = mid;
if( c < 0.0)
if(e > 0.0) b = mid;
else a = mid;
c = f(a); d = f(b);
}
return mid;
}