仅供参考:
// 牛顿法解多项式的根
// 输入:多项式系数c[],多项式度数n,求在[a,b]间的根
// 输出:根
// 要求保证[a,b]间有根
#include <stdio.h>
#include <stdlib.h>
#define MAXN1 10
double C[MAXN1];
double ROOT;
double fabs( double x ) {
return (x<0)? -x : x;
}
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;
}
void main() {
C[0]=-0.5;
C[1]=0.0;
C[2]=1;
ROOT=Polynomial_Root(C, 2, -1.0, 0.0, 0.0001);
printf("%g\n",ROOT);//x*x-0.5=0在[-1,0]之间的根精确到0.0001约是-0.707107
}