数值分析之牛顿法多项式求根

很朴素的方法，如果在区间[a,b]内有根，那么f(a)*f(b)<0

不难得到以下代码：

 

#include <iostream>
#include <memory>

using namespace std;

double f(int m, double c [], double x)
{
	int i;
	double p = c[m];

	for (i = m; i > 0; i--)
		p = p*x + c[i - 1];
	return p;
}

int newton(double x0, double *r,
	double c [], double cp [], int n,
	double a, double b, double eps)
{
	int MAX_ITERATION = 1000;
	int i = 1;
	double x1, x2, fp, eps2 = eps / 10.0;

	x1 = x0;
	while (i < MAX_ITERATION) {
		x2 = f(n, c, x1);
		fp = f(n - 1, cp, x1);
		if ((fabs(fp)<0.000000001) && (fabs(x2)>1.0))
			return 0;
		x2 = x1 - x2 / fp;
		if (fabs(x1 - x2) < eps2) {
			if (x2<a || x2>b)
				return 0;
			*r = x2;
			return 1;
		}
		x1 = x2;
		i++;
	}
	return 0;
}

double Polynomial_Root(double c [], int n, double a, double b, double eps)
{
	double *cp;
	int i;
	double root;

	cp = (double *) calloc(n, sizeof(double) );
	for (i = n - 1; i >= 0; i--) {
		cp[i] = (i + 1)*c[i + 1];
	}
	if (a > b) {
		root = a; a = b; b = root;
	}
	if ((!newton(a, &root, c, cp, n, a, b, eps)) &&
		(!newton(b, &root, c, cp, n, a, b, eps)))
		newton((a + b)*0.5, &root, c, cp, n, a, b, eps);
	free(cp);
	if (fabs(root) < eps)
		return fabs(root);
	else
		return root;
}

int main()
{
	double c[2] = { 1,2 };
	int n = 1, a = -5, b = 5;
	double eps = 1e-8;
	double root = Polynomial_Root(c, n, a, b, eps);
	cout << root << endl;
}