newtonmodify.cpp

1.大型稀疏线性方程组的求解 A*X=b 。 2. 一维数组冒泡法排序算法 4.矩阵求逆 5. 改进的牛顿算法——弦割法

//=============================================================
//                                                            =
//              改进的牛顿算法——弦割法                      =
//                                                            =
//                                                            =
//                                                            =
//                                                            =
//                                                            =
//	与牛顿法相比，最大优点是避免了求导数。但收敛速度不如  =
//  牛顿法,且需提供两个较好的初始近似根x0、x1.                =
//                                                            =
//=============================================================

//=============================================================
//
//	参数说明：
//		xl	区间左端点
//		xh	区间右端点
//		eps	用户给定的精度要求
//		*ErrorFlag	指针，错误标志。为 -1 时，表示根不在区间内，
//                  此时，根被赋于一个很大的值 
//
//=============================================================


#include <math.h>
#include <stdio.h>

double RegFalsi(double xl,double xh,double eps,int *ErrorFlag);

void main()
{
	double xl,xh,root,eps;
	int ErrorFlag;

	xl=0;
	xh=2.0;
	eps=1.0e-8;

	root=RegFalsi(xl,xh,eps,&ErrorFlag);

	printf("ErrorFlag = %d\n",ErrorFlag);
	printf(" Root = %f\n",root);
}

/*
double AFunction(double x)	//自编函数
{
	return(2*x*exp(-5)+1-2*exp(-5*x));
}
*/


double AFunction(double x)	//自编函数
{
	return(x*x*x-x-1);
}


double RegFalsi(double xl,double xh,double eps,int *ErrorFlag)
{
	double x0,x1,x2,f0,f1,f2;

	*ErrorFlag=0;
	x0=xl;
	f0=AFunction(x0);

	if(fabs(f0)<eps) return(x0);

	x1=xh;
	f1=AFunction(x1);

	if(fabs(f1)<eps) return(x1);

	if(f0*f1>0)
	{
		*ErrorFlag=-1;
		return(1.0e+15);
	}

	while(1)
	{
		x2=x1-(x0-x1)/(f0-f1)*f1;
		f2=AFunction(x2);

		if(fabs(f2)<eps) return (x2);

		if(f1*f2<0) 
		{
			x0=x1;
			f0=f0/(f0+f2)*f1;
		}
		else 
			f0=f1/(f1+f2)*f0;

		while(1)
		{
			x1=x2;
			f1=f2;

			if(fabs(x0-x1)<eps)
			{
				if(fabs(f0)<fabs(f1))
					return(x0);
				else
					return(x1);
			}

			x2=x1-(x0-x1)/(f0-f1)*f1;
			f2=AFunction(x2);

			if(fabs(f2)<eps) return(x2);

			if(f1*f2>0)
			{
				f0=f1/(f1+f2)*f0;
				continue;
			}
			break;
		}

		x0=x1;
		f0=f1;
		x1=x2;
		f1=f2;

		if(fabs(x0-x1)<eps)
			break;
	}

	if(fabs(f0)<fabs(f1))
		return(x0);
	else
		return(x1);
}