thesecond.cpp

这个程序分别用二分法

/* 04计B  林传敏  2004131301 */

#include <iostream.h>
#include <iomanip.h>
#include <math.h>
const int N=45;

double x[10];

int checkRoot(double);								// 检查某个根是否为重根			
int binary(double, double);							// 二分法
int Newton(double, double);							// Newton法
int xianjie(double, double);						// 弦截法
int NewtonDownMou(double, double);					// Newton下山法

int main()
{
	cout << "请输入运行次数 t:";
	int t;
	cin >> t;
	for(int kk=0; kk < t; kk++)
	{
		cout << "Input the range of the root" << endl;
		int		a, b,count;
		cin >> a >> b;
		cout << "the range is : [" << a << "," << b << "]" << endl;

		cout << "m != 0,说明这个数是方程的m重根, or, 这个数不是方程的根." << endl;
	
		binary(a, b);
		cout.setf(ios::left, ios::adjustfield);
		cout << "二分法" << endl;
		cout.width(10);
		cout << "roots" << "\t" << "m" << endl << "-----------------" << endl;
		for(int i=0; i < 5; i++)
		{
			cout.width(10);
			cout << (x[2*i] + x[2*i+1]) / 2 << "\t";
			count = checkRoot((x[2*i] + x[2*i+1]) / 2);
			cout << count << endl;
		}  // for
		cout << "\n";

		Newton(a, b);
		cout << "Newton法" << endl;
		cout.width(10);
		cout << "roots" << "\t" << "m" << endl << "-----------------" << endl;
		for(i=0; i < 5; i++)
		{
			count = checkRoot(x[i]);		
			cout.width(10);
			cout << x[i] << "\t" << count << endl;
		} // for

		cout << "\n";

		xianjie(a, b);
		cout << "XianJie法" << endl;
		cout.width(10);
		cout << "roots" << "\t" << "m" << endl << "-----------------" << endl;
		for(i=0; i < 5; i++)
		{
			count = checkRoot(x[i]);		
			cout.width(10);
			cout << x[i] << "\t" << count << endl;
		} // for
		cout << "\n";

		NewtonDownMou(a, b);
		cout << "Newton下山法" << endl;
		cout.width(10);
		cout << "roots" << "\t" << "m" << endl << "-----------------" << endl;
		for(i=0; i < 5; i++)
		{
			count = checkRoot(x[i]);		
			cout.width(10);
			cout << x[i] << "\t" << count << endl;
		} // for
	} // for
	
	cout << "Bye-Bye!!!" << endl;
	return 0;
} // main


int checkRoot(double x0)							//检查某个根是否为重根
{
	double	x, y, y0;
	x = x0 + 0.001;
	y = pow(x, 5) - 3*pow(x, 3) + x - 1;
	y0 = pow(x0, 5) - 3*pow(x0, 3) + x0 - 1;

	int	i=0;
	double k;
	if(fabs(y0) <= 0.0001)
	{
		i++;	
		k = (y0-y) / (x0-x);						// 一阶导数
		if(fabs(k) <= 0.00001)	
		{
			i++;
			y = 5*pow(x, 4) - 9*pow(x, 2) + 1;
			y0 = 5*pow(x0, 4) - 9*pow(x0, 2) + 1;
	
			k = (y0-y) / (x0-x);					// 二阶导数
			if(fabs(k) <= 0.00001)
			{
				i++;
				y = 20*pow(x, 3) - 18*x;
				y0 = 20*pow(x0, 3) - 18*x0;
	
				k = (y0-y) / (x0-x);				// 三阶导数
				if(fabs(k) <= 0.00001)
				{
					i++;
					y = 60*pow(x, 2) - 18;
					y0 = 60*pow(x0, 2) - 18;

					k=(y0-y) / (x0-x);				// 四阶导数
					if(fabs(k) <= 0.00001)
						i++;
				}
			}
		}
	}

	return i;
}

int binary(double a, double b)						// 二分法
{
	int i,j;
	double n,x1,x2,y1,y2;

	for(i=0; i < 10; i++)
	{
		x[i]=0;
	}

	n = fabs(a-b) / N;
	for(i=0,j=0; (i<N) && (j<10); i++)
	{
		x1 = a + i*n;
		x2 = a + (i+1)*n;
		y1 = pow(x1, 5) - 3*pow(x1, 3) + x1 - 1;
		y2 = pow(x2, 5) - 3*pow(x2, 3) + x2 - 1;

		if((y1<0) && (y2>0))
		{
			x[j++] = x1;
			x[j++] = x2;
		}
		else if((y1>0) && (y2<0))
		{
			x[j++] = x2;
			x[j++] = x1;
		}
	} // for

	for(i=0; i < 5; i++)
	{
		while(fabs(x[2*i]-x[2*i+1]) > 0.00001)
		{

			x1 = (x[2*i] + x[2*i+1]) / 2;
			y1 = pow(x1, 5)-3*pow(x1, 3) + x1 - 1;

			if(y1 < 0)
			{
				x[2*i] = x1;
			}
			else if(y1 > 0)
			{
				x[2*i+1] = x1;
			}
		}
	}

	return 0;
}

int Newton(double a, double b)						// Newton法
{
	int			i,j;
	double		n,x1,x2,y1,y2,k;

	for(i=0; i < 10; i++)
	{
		x[i] = 0;
	}

	n = fabs(a-b) / N;
	for(i=0,j=0; (i < N) && (j < 5); i++)
	{
		x1 = a + i*n;
		x2 = a + (i+1)*n;
		y1=pow(x1, 5) - 3*pow(x1, 3) + x1 - 1;
		y2=pow(x2, 5) - 3*pow(x2, 3) + x2 - 1;

		if((y1<0) && (y2>0))
		{
			x[2*j] = x1;
			x[2*j+1] = x2;
			j++;

		} // if
		else if((y1>0) && (y2<0))
		{
			x[2*j] = x2;
			x[2*j+1] = x1;
			j++;
		} // else
	} // for

	for(i=0; i < j; i++)
	{
		x1 = x[2*i];		
		y1 = pow(x1, 5) - 3*pow(x1, 3) + x1 - 1;

		x2 = x[2*i] + 0.00001;
		y2 = pow(x2, 5) - 3*pow(x2, 3) + x2 - 1;

		k = (y2-y1) / (x2-x1);						//一次导数
		x2 = x1 - y2/k;

		while(fabs(x1-x2) > 0.00001)
		{
			y1 = pow(x2, 5) - 3*pow(x2, 3) + x2 - 1;

			x1 = x2;
			x2 = x2 + 0.00001;
			y2 = pow(x2, 5) - 3*pow(x2, 3) + x2 - 1;
			k = (y2-y1) / (x2-x1);

			x2 = x1 - y1/k;
		} // while 

		x[i] = x2;
		x[i+1] = 0;
	} // for
	x[i+1] = 0;

	return 0;
}

int xianjie(double a, double b)						// 弦截法
{
	double		n,x1,x2,y1,y2,temp;

	for(int i=0; i < 10; i++)
	{
		x[i]=0;
	} // for

	n = fabs(a-b)/N;
	int j;
	for(i=0,j=0; (i < N) && (j < 5); i++)
	{
		x1 = a + i*n;
		x2 = a + (i+1)*n;
		y1 = pow(x1, 5) - 3*pow(x1, 3) + x1 - 1;
		y2=pow(x2, 5) - 3*pow(x2, 3) + x2 - 1;

		if((y1 < 0) && (y2 > 0))
		{
			x[2*j] = x1;
			x[2*j+1] = x2;
			j++;

		} // if
		else if((y1 > 0) && (y2 < 0))
		{
			x[2*j] = x2;
			x[2*j+1] = x1;
			j++;
		} // else
	} // for

	for(i=0; i < j; i++)
	{
		x1 = x[2*i];
		x2 = x[2*i] + 0.1;

		y1 = pow(x1, 5) - 3*pow(x1, 3) + x1 - 1;
		y2 = pow(x2, 5) - 3*pow(x2, 3) + x2 - 1;
		
		temp = x2;
		x2 = x2 - y2*(x2-x1)/(y2-y1);
		x1 = temp;

		while(fabs(x2-x1) > 0.00001)
		{
			y1 = pow(x1, 5) - 3*pow(x1, 3) + x1 - 1;
			y2 = pow(x2, 5) - 3*pow(x2, 3) + x2 - 1;

			temp = x2;
			x2 = x2 - y2*(x2-x1)/(y2-y1);
			x1 = temp;
		}  // while

		x[i] = x2;
		x[i+1] = 0;
	} // for
	
	x[i+1] = 0;

	return 0;
} // xianjie



int NewtonDownMou(double a, double b)				// Newton下山法
{
	int			i, j, flag=0;
	double		n, x1, x2, y1, y2, temp, k, r;

	for(i=0; i < 10; i++)
	{
		x[i]=0;
	} // for 

	n = fabs(a-b) / N;
	i = 0;
label:
	for(r=1,j=0; i < N; i++)
	{
		x1 = a + i*n;
		y1 = pow(x1, 5) - 3*pow(x1, 3) + x1 - 1;

		x2 = x1 + 0.000001;
		y2 = pow(x2, 5) - 3*pow(x2, 3) + x2 - 1;

		k = (y2-y1) / (x2-x1);
		x2 = temp = x1 - r*y1/k;
		y2 = pow(x2, 5) - 3*pow(x2, 3) + x2 - 1;

		while(1)
		{
			if(fabs(y2) < fabs(y1))
			{
				x2 = temp;
				break;
			}
			else
			{
				if(r < 0.00001)
				{
					cout << "Fail to down moutain "
						 << "choose the oher first_value"
						 << endl;
					
					goto label;
				}
				r = r*0.5;
				temp = x1 - r*y1/k;
				y2 = pow(temp, 5) - 3*pow(temp, 3) + temp - 1;
			}
		} // while

		for(r=1; fabs(x2-x1) > 0.00001; r=1 )
		{
			x1 = x2;
			y1 = pow(x1, 5) - 3*pow(x1, 3) + x1 - 1;

			x2 = x2 + 0.000001;
			y2 = pow(x2, 5) - 3*pow(x2, 3) + x2 - 1;

			k = (y2-y1) / (x2-x1);
			x2 = temp = x1 - r*y1/k;
			y2 = pow(x2, 5) - 3*pow(x2, 3) + x2 - 1;

			while(1)
			{
				if(fabs(y2) < fabs(y1))
				{
					x2 = temp;
					break;
				}  // if
				else
				{
					if(r < 0.00001)
					{
						cout << "Fail to down moutain"
							 << "choose the oher first_value"
							 << endl;
						
						goto label;
					} // if
					r = r*0.5;
					temp = x1 - r*y1/k;
					y2 = pow(temp, 5) - 3*temp*temp*temp+temp-1;
					
				} // else
			} // while
		} // for

		j = 0;
		while(j < 10)
		{
			if(fabs(x2-x[j]) > 0.1)
			{
				j++;
			} // if
			else
				break;
		} // while
		
		if(j == 10)
		{
			x[flag++] = x2;
		} // if
	}
	return 0;
} // NewtonDownMoutain