fixedpointiter.c

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

/*   This program is about the root-finding algorithm "Fixed Point Iteration", 
 * which is described in detail below. This program depends on the file 
 * "funct.c" in which the function g(x) is defined, you can alter the function
 * g(x) as you like, but be sure to put these two files together.
 *   This program is compiled using gcc 3.4.2 on Windows operating system 
 * with mingw32, 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 FixedPointIter.c -o fpiter -lm
 * ./fpiter | tee dat
 * gnuplot
 * In gnuplot:
 * plot "dat" using points
 * Then the specified information 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]
 *
 * Last modified by Xiaoke Yang, Nov.5,2007
 *
*/
#include <math.h>
#include <stdio.h>
#include "funct.c"		/*funct.c includes g(x)*/
#define eps		1e-20	/*accuracy of g(x) or f(x), the error tolerance*/
#define IMAX	1000	/*maximum iteration step*/
#define DEBUG	1		/*debug flag, set 0 to turn off the debug info*/


/* FUNCTION: double FixedPointIter(double (*g)(double),double p0)
 * DESCRIPTION: This procedure uses the "Fixed Point Iteration" algorithm to 
 * solve fixed point problems or root-finding problems. The equation f(x)=0 
 * and the corresponding iterated function g(x)=x are defined in a seperate
 * file "funct.c". The aim of this procedure is to find a real value p such 
 * that f(p)=0 or g(p)=p, with a initial guess p0. The following conditions
 * assure the convergence of this algorithm. 
 * For a given interval [a,b]:
 * (1) g belongs to C[a,b](the set of all continuous functions defined on [a,b])
 *  and g(x) belongs to [a,b] for any x in [a,b];
 * (2) g'(x) exists on (a,b) and a positive constant K<1 exists with
 *  |g'(x)|<=K<1;
 * then this procedure will generate a sequence pn which converges 
 * linearly(when g'(p)!=0) to the unique fixed point in [a,b], when
 * g'(p)=0, higher-order convergence can be gotten.
 * The bounds of the error given by this algorithm is 
 * |p-pn|<=(K^n)/(1-K)|p1-p0|, in which n is the iteration number.
 * INPUT: 	g- a function pointer to g(x)
 * 			p0- the initial guess point
 * OUTPUT:	the unique fixed point of x=g(x) or root of f(x) 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=FixedPointIter(g,p0);
 * CREATED DATE: Nov.5,2007
 * CREATED BY: Xiaoke Yang
 * */
double FixedPointIter(double (*g)(double),double p0){
	int i;						/*the iterator*/
	double x,gx,err;			/*gx=g(x), err is the relative error of g(x)*/
	x=p0;						/*initialize x*/
	for(i=1;i<=IMAX;i++){		/*iteration starts*/
		gx=g(x);				/*calculate the right hand side of x=g(x)*/
		err=fabs((x-gx)/gx);	/*calculate the error*/
		if(err<eps){			/*stopping cirteria*/
			if(DEBUG==1){		/*print the debug information if needed*/
				printf("#Program successfully exits with the given accuracy achieved.\n#relative error=%0.15e<=%0.15e.\n",err,eps);
				printf("#The root of f(x) is x=\n");
				printf("%0.15e\n",gx);
			}
			return gx;			/*return the root*/
		}	
		if(DEBUG==1){			/*print the debug information*/
			printf("#The current iteration count is:\n");
			printf("%d\n",i);
			printf("#The point in the previous step is x=\n");
			printf("%0.15e\n",x);
			printf("#g(x) in this step is g(x)=\n");
			printf("%0.15e\n",gx);
			printf("#The relative error of x in this step is log10(err)=\n");
			printf("%0.15e\n\n",log10(err));
		}
		x=gx;
	}
	/*iterations exceeded, print an ERROR MESSAGE on stand error*/
	fprintf(stderr,"#Maximum steps of iteration exceeded! The result may not be accurate enough!\n");
	fprintf(stderr,"#Consider increase the maximum iteration step.\n");
	if(DEBUG==1){		/*print the debug information if needed*/
		printf("#The root returned at last is x=\n");
		printf("#%0.15e\n",x);
		printf("#The relative error of x at last is log10(err)=\n");
		printf("%0.15e\n",log10(err));
	}
	return x;			/*return the root, although NOT accurate enough*/
}

int main(void){
	FixedPointIter(g,10000);
}