newton.cpp

Newton法求解方程的非重根

#include <stdio.h>
#include <math.h>
#define ZERO 0.000000001 /* X is considered to be 0 if |X|<ZERO */
#define MAXN 11   /* Max Polynomial Degree + 1 */

double Polynomial_Root(int n, double c[], double a, double b, double EPS);
int main()
{	
	freopen("input.txt","r",stdin);
	freopen("output.txt","w",stdout);
	int n;
	double c[MAXN], a, b;
	double EPS = 0.00005;
	int i;
	
	while (scanf("%d", &n)!= EOF){
		for (i=n; i>=0; i--) 
			scanf("%lf", &c[i]);
		scanf("%lf %lf", &a, &b);
		printf("%.4lf\n", Polynomial_Root(n, c, a, b, EPS));
	}
 return 0;
}

/* Your function will be put here */

double newton(int n, double c[], double d[], double a, double b, double EPS, double p0,int * flag);
double func(int n, double c[], double x);


double Polynomial_Root(int n, double c[], double a, double b, double EPS)
{
	double d[MAXN];

	int flag=0;
	if(a>b)
	{
		double temp;
		temp=a;
		a=b;
		b=temp;
	}
	int i;
	for (i=n;i>=1;i--)
			d[i-1]=i*c[i];              //diff
	
	double p0= (a+b)/2;
	
	double ans=newton( n,c,d,a,b,EPS/10, a,&flag);
	
	if(!flag)
			 ans=newton( n,c,d,a,b,EPS/10, b,&flag);
	if(!flag)
			 ans=newton( n,c,d,a,b,EPS/10, p0,&flag);
	flag=0;
	if (fabs(ans)<EPS)
		return(fabs(ans));
	else
	return ans;
}
double newton(int n, double c[], double d[], double a, double b, double EPS, double p0,int *flag)
{
	double p,t;
	do 
	{
		double x=func(n,c,p0);
		double y=func(n-1,d,p0);

		if (fabs(y)<ZERO && fabs(x)>1)
			return 0;
		p=p0-x/y;
		
		t=p0-p;
		
		if( fabs(t)< EPS )
			{
			*flag=1;
			return(p);
			}
		else
			p0=p;
		
	}while (p>a&&p<b);
	return 0;
}

double func(int n, double c[], double x)
{
		int i;
		double sum=c[n];
		for(i=n;i>0;i--)
		{
			sum=sum*x+c[i-1];
		}
		return sum;
}