问题五十九：怎么求一元六次方程在区间内的所有不相等的实根（2）

59.2 求一元六次方程在区间内的所有不相等的实根

59.2.1 理论分析

在“59.1”中，我们已经求得一元六次方程在区间内不相等的实根的总个数。接下来，我们将具体求出方程在区间内的所有不相等的实根。

分两步进行：

第一步：将区间分成若干个只包含一个不相等实根的子区间。

1.1，计算整个区间[l,r]内实根的个数N（59.1中已经完成）；

1.2，若该区间内实根的个数>1，用二分法将该区间对分，计算[l,(l+r)/2]内实根的个数；

1.3，若新区间内实根的个数>1，继续二分区间，计算新二分区间内实根的个数；

若新区间内实根的个数=0，计算[(l+r)/2,r]内实根的个数；

第二步：若新区间内实根的个数=1，用牛顿迭代法求得该实根。

在如下章节有介绍过牛顿迭代法：

“问题四十六：怎么用ray tracing画superellipsoid”

http://blog.csdn.net/libing_zeng/article/details/54564227

注意：

牛顿迭代法，求得的实根不一定是区间内的实根。

比如：

某一元六次方程在区间[-5,5]有-4、-3、-2、2、3、4六个实根，在区间[-2.5, -1.5]进行迭代时，我们是想求”-2”这个实根，但是迭代的结果可能是”2”。

这时，我们的做法是：舍弃结果，返回迭代失败，求根总次数+1（因为失败了一次，需要再来一次），区间二分回到第一步重新开始（理论上，只要区间足够小，一定可以迭代到满足精度要求的实根。也有例外情况，如果区间小到误差范围，又知道该区间内有一个实根，直接将区间端值作为要求的实根）

另外，数据类型最好用double

59.2.2 看C++代码实现

----------------------------------------------vec3.h ------------------------------------------

vec3.h

bool roots_equation_6th(float ee6[7], float a, float b, float tol, float (&roots)[7]);

----------------------------------------------vec3.cpp ------------------------------------------

vec3.cpp

bool get_equation_6th_function_and_derivative(float ee6[7], double xxx, double& fnc, double& drv) {/*determine function value and derivative value at xxx*/fnc = 0.0;drv = 0.0;for (int i=0; i<7; i++) {fnc = fnc + ee6[i]*pow(xxx, (6-i));if (i<6) {drv = drv + (6-i)*ee6[i]*pow(xxx, (5-i));}}return true;
}bool get_root_by_newton_iteration_for_equation_6th(float ee6[7], double x0[2], float tol, double (&root)[2]) {/*find the single root in the interval [x0[0], x0[1]] by newton iteration */double t_k, t_k1, ft_k, ft_d_k;int get = 0;double i0 = x0[0];double i1 = x0[1];for (int i=0; i<2; i++) {t_k = x0[i];for (int k=0; k<20; k++) {if (!(isnan(t_k))) {get_equation_6th_function_and_derivative(ee6, t_k, ft_k, ft_d_k);if ((ft_d_k != 0) && !(isnan(ft_k)) && !(isnan(ft_d_k))) {t_k1 = t_k - ft_k/ft_d_k;if (((fabs(t_k1)>=1e-2)&&(fabs((t_k1 - t_k)/t_k1) < tol)) ||((fabs(t_k1)<1e-2)&&(fabs((t_k1 - t_k)) < tol))) {if ((t_k1 > x0[0]) && (t_k1<=x0[1])) {root[1] = (fabs(t_k1)<(tol*10))? 0.0:t_k1;root[0] = 1.0;get ++;}else {x0[1] = (x0[0] + x0[1])/2;}break;}else {t_k = t_k1;}}else {break;}}else {break;}}if (get == 1) {break;}}root[0] = double(get);return true;
}bool roots_equation_6th(float ee6[7], float a, float b, float tol, float (&roots)[7]) {double root[2], xxx0[2], temp;double left = a;double right = b;int total_num, interval_num;
/*”total_num” for the whole interval; “interval_num” for every sub-interval*/int offset = 0;roots[0] = 0.0;roots_num_equation_6th(ee6, double(a), double(b), total_num);
/*determine the number of distinct real roots in the interval*/if (total_num == 0) {return true;}else {interval_num = total_num;for (int i=1; i<(total_num+1+offset); i++) {while(interval_num != 1) {
/*to find a sub-interval with one real root in*/if (left > right) {temp = left;left = right;right = temp;}right = (left + right)/2;roots_num_equation_6th(ee6, left, right, interval_num);if (interval_num == 0) {left = right;right = b;}}xxx0[0] = left;xxx0[1] = right;get_root_by_newton_iteration_for_equation_6th(ee6, xxx0, tol, root);if (root[0] != 0.0) {
/*if newton finds the root*/roots[0] = roots[0] + 1.0;roots[int(roots[0])] = root[1];left = right;right = b;}else if ((root[0] == 0.0)&&(fabs(left-right)<tol)&&(interval_num==1)) {
/*if newton doesn’t find the root, but the interval is small enough*/roots[0] = roots[0] + 1.0;roots[int(roots[0])] = left;left = right;right = b;}else {
/*newton failed to find the root in regular interval, we need another time to find the root with smaller interval*/offset ++;right = (left + right)/2;}roots_num_equation_6th(ee6, left, right, interval_num);}}return true;
}

----------------------------------------------main.cpp ------------------------------------------

main.cpp

    int main(){
//        float ee6[7] = {1.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0};//x^6------------1
//        float ee6[7] = {1.0, -1.0,  0.0,  0.0,  0.0,  0.0,  0.0};//x^5*(x-1)------------2
//        float ee6[7] = {1.0, 0.0,  -1.0,  0.0,  0.0,  0.0,  0.0};//x^4*(x-1)*(x+1)------------3
//        float ee6[7] = {1.0,  0.0,  -5.0,  0.0,  4.0,  0.0,  0.0};//x^2*(x-1)*(x+1)*(x-2)*(x+2)------------5
//        float ee6[7] = {1.0,  3.0,  -5.0, -15.0,  4.0,  12.0,  0.0};//x*(x-1)*(x+1)*(x-2)*(x+2)*(x+3)------------6float ee6[7] = {1.0,  0.0, -0.4236111, 0.0,  0.050347222,  0.0, -0.0017361111};//(x-1/2)*(x+1/2)*(x-1/3)*(x+1/3)*(x-1/4)*(x+1/4)------------6float a = -2.5;float b =  2.5;int num;float tol = 1e-6;float roots[7];roots_num_equation_6th(ee6, a, b, num);std::cout << "the coefficients of the equation of 6th degree are:" << endl;for (int i=0; i<7; i++) {std::cout << ee6[i] << "  ";}std::cout << endl;std::cout << "the number of the roots of this equation is: " << num << endl;if (num > 0) {roots_equation_6th(ee6, a, b, tol, roots);std::cout << "the roots are: ";for (int j=1; j<(num+1); j++){std::cout << roots[j] << " ";}std::cout << endl;}
}

输出结果如下：

float ee6[7] ={1.0, 0.0, -0.4236111, 0.0, 0.050347222, 0.0, -0.0017361111};

//(x-1/2)*(x+1/2)*(x-1/3)*(x+1/3)*(x-1/4)*(x+1/4)------------6

float ee6[7] ={1.0, 3.0, -5.0, -15.0, 4.0, 12.0, 0.0};

//x*(x-1)*(x+1)*(x-2)*(x+2)*(x+3)------------6

float ee6[7] ={1.0, 0.0, -5.0, 0.0, 4.0, 0.0, 0.0};

//x^2*(x-1)*(x+1)*(x-2)*(x+2)------------5

float ee6[7] ={1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};//x^6------------1

OK, 在这里就只测试这么几个方程。

后续在“用ray tracing画回旋体”时，生成每个图形可能会有成千上万次求解一元六次方程，区间在[0,1]。

by the way, 此刻发文时，已经完成“回旋体”图形的生成。所以，如上求解一元六次方程在区间内的所有不相等实根的程序在一定程度上是可靠的。