bisection.c

一些迭代求根法的C语言程序（简单迭代法、割线法、二分法等）些迭代求根法的C语言程序（简单迭代法、割线法、二分法等）

/*   This program is about the root-finding algorithm "bisection method", 
 * which is described in detail below. This program depends on the file 
 * "funct.c" in which the function f(x) is defined, you can alter the function
 * f(x) as you like, but be sure to put these two files together.
 *   This program is compiled using gcc 4.1.3 on ubuntu Linux operating system, 
 * you can also compile and run it with other standard C/C++ compiler and 
 * on other operating systems. The output of the debug information is intended
 * for visualization using gnuplot. The following steps tell how to do it.
 * In Linux shell(bash):
 * gcc bisection.c -o bisection -lm
 * ./bisection | tee dat
 * gnuplot
 * In gnuplot:
 * plot "dat" using points
 * Then the convergency process with respect to the iteration step will 
 * be displayed. :) Have fun!
 *
 *   This program was written by Xiaoke Yang @ School of Automation Science and
 * Electrical Engineering, Beihang University. You can copy, modify or 
 * redistribute it freely with this header kept intact. Welcome reports of 
 * bugs, suggestions, etc. [e-mail] 
 * xiaoke
 * 
 * Created by Xiaoke on Nov.1, 2007
 * Last modified by Xiaoke Yang on Nov.5, 2007
 *
 * */
#include <math.h>
#include <stdio.h>
#include "funct.c"		/*funct.c includes f(x)*/
#define eps		1e-20	/*accuracy of f(x), the absolute error tolerance*/
#define delta	1e-20	/*accuracy of x, the absolute error tolerance*/
#define IMAX	1000	/*maximum iteration step*/
#define DEBUG	1		/*debug flag, set 0 to turn off the debug info*/


/* FUNCTION: double bisection(double (*f)(double), double a, double b)
 * DESCRIPTION: Bisection method is a root-finding algorithm which works 
 * repeately dividing an interval in half and selecting the subinterval 
 * in which the root exitst. In this procedure, the equation is defined as 
 * f(x)=0 which is passed by a function pointer. Two points a and b should 
 * also be provided such that f(a) and f(b) have opposite signs, then 
 * this procedure can always converge to a root lying in the 
 * interval [a,b] with the predefined error tolerance.
 * The bisection method is less efficient than Newton's method but it is 
 * much less prone to odd behavior.
 * INPUT:	f- a function pointer to f(x) 
 *			a- left bound of closed interval [a,b]
 *			b- right bound of closed interval [a,b]
 * OUTPUT:	one of the roots lying in [a,b], PAY ATTENTION TO THE ERROR MESSAGE
 *			if there exists on the standard error which might be indicating the
 *			errors of computation process.
 * CALLING SYNTAX: val=bisection(f,a,b);
 * CREATED DATE: Nov.1, 2007
 * CREATED BY: Xiaoke Yang
 * */
double bisection(double (*f)(double), double a, double b){
	int i;					/*the iterator*/
	double an,bn;			/*the boundary of the interval in each iteration*/
	double m,fm;			/*m:=the mid-point of an and bn. fm:=f(m)*/
	an=a;					/*the initial boundary, passed as arguments*/
	bn=b;
	/*check the initial interval to make sure a<b*/
	if(a>=b){
		fprintf(stderr,"#INPUT ERROR!\n");
		fprintf(stderr,"#a=%e,b=%e\n",a,b);
		fprintf(stderr,"#Make sure a<b for the initial interval [a,b]!\n");
		return 0;
	}
	/*the roots' number in the interval should be odd*/
	if(f(a)*f(b)>0){
		fprintf(stderr,"#f(%f)f(%f)>0, computation failed!",a,b);
		return 0;			/*in this case, 0 is NOT the root*/
	}
	/*check whether the left bound "a" is the root*/
	if(fabs(f(a))<eps){		/*f(a) is less than the error tolerance*/
		if(DEBUG==1){
			printf("#The root of f(x) is x=\n");
			printf("%e\n",a);
		}
		return a;			/*return a as the root*/
	}
	/*check whether the right bound b is the root*/
	if(fabs(f(b))<eps){		/*f(b) is less than the absolute error tolerance*/
		if(DEBUG==1){
			printf("#The root of f(x) is x=\n");
			printf("%e\n",b);
		}
		return b;			/*return b as the root and terminate*/
	}
	/*start the iteration loop*/
	for(i=1;i<=IMAX;i++){
		m=(an+bn)/2;		/*get the mid-point*/
		fm=f(m);			/*evaluate f(m)*/
		if(fabs(fm)<eps){	/*whether the mid-point is the root*/
			if(DEBUG==1){
				printf("#The current iteration count is:\n");
				printf("%d\n",i);
				printf("#The error tolerance of f(x) achieved!\n");
				printf("#The root of f(x) is x=\n");
				printf("#%0.15e\n",m);
				printf("#The corresponding f(x)=\n");
				printf("#%0.15e\n",fm);
			}
		return m;			/*if so, return it*/
		}
		if(f(an)*f(m)<0)	/*else select the new sub-interval*/
			bn=m;
		else
			an=m;
		if(bn-an<delta){	/*check the length of the interval*/
			if(DEBUG==1){
				printf("#The lower limit of the interval achieved!\n");
				printf("#The root of f(x) is x=\n");
				printf("#%e\n",m);
				printf("#The corresponding function value is f(x)=\n");
				printf("#%e\n",fm);
			}
			return m;/*if small enough, return the mid-point as the root*/
		/*else print debug information of the current step if needed*/
		}else{	
			if(DEBUG==1){
				printf("#The current iteration count is:\n");
				printf("%d\n",i);
				printf("#The value finded in this step is x=\n");
				printf("%0.15e\n",m);
				printf("#The corresponding function value is f(x)=\n");
				printf("#%0.15e\n\n",fm);
			}
		}
	}
	/*iteration finishes, indicating the specified accuracy is not reached.*/
	fprintf(stderr,"#Maximum steps of iteration reached! The result might not be accurate enough!\n");
	fprintf(stderr,"#Consider increase the maximum iteration step.\n");
	if(DEBUG==1){
		printf("#The value returned at last is x=\n");
		printf("#%e\n",m);
		printf("#The corresponding function value is f(x)=\n");
		printf("#%e\n",fm);
	}
	return m;		/*return the midpoint, although NOT accurate enough*/
}
/*the main function, for the purpose of testing this algorithm*/
int main(void){
	double x;
	x=bisection(f,0,3.5);
	printf("#One of the roots/or the root is %f\n",x);
	return 0;
}