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

为什么要求一元六次方程在某区间的所有根？

原因是：

后面在用ray tracing画回旋体（rotational sweeping/ revolution）时，若侧面曲线是三次b样条曲线，求光线和回旋体的交点时会出现一元六次方程，而且我们要求的是离光线起点最近的交点，所有，我们需要先求出所有交点，然后选出最近的交点。

 

接下来，我们会分两步来做这个事情：

第一步：判断方程在区间内不相等的实根的总个数N；

第二步：求出所有N个不相等的实根

 

59.1判断一元六次方程在区间内不相等的实根的总数

59.1.1 理论分析

对于多项式方程在区间内不等实根总数的判断，我们可以参照Sturm's theorem。

维基百科：https://en.wikipedia.org/wiki/Sturm%27s_theorem

 

判断过程分两步:

（为什么要这么做？参考如上链接，研究一下“Sturm’s theorem”,我们这一章节的内容相当于“Sturm’s theorem”的一个实例）

 

首先，构造一个叫做“sturm序列”的东东。

 

 

关于“多项式长除法”，看这里：

https://en.wikipedia.org/wiki/Polynomial_long_division

 

然后，将区间两端点值l,r代入各个序列项，计算符号的变化次数。

59.1.2 看C++代码实现

 

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

vec3.h

bool roots_num_equation_6th(float ee6[7], double a, double b, int &num);

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

vec3.cpp

bool roots_num_equation_6th(float ee6[7], double a, double b, int &num) {
        double temp1, temp2;
        double ee77[7][7] = {0.0};
        double fun_a[7] = {0.0};
        double fun_b[7] = {0.0};
        int change_num_a = 0;
        int change_num_b = 0;

        for (int i=0; i<6; i++) {
        //determine f0, f1
            ee77[0][i] = ee6[i];
            ee77[1][i] = (6-i)*ee6[i];
        }
        ee77[0][6] = ee6[6];

        for (int i=2; i<6; i++) {
        //determine f2, f3, f4, f5
            temp1 = ee77[i-2][0] / ee77[i-1][0];
            temp2 = (ee77[i-2][1] - temp1*ee77[i-1][1]) / ee77[i-1][0];
            for (int j=0; j<(6-i); j++) {
                ee77[i][j] = temp1*ee77[i-1][j+2] + temp2*ee77[i-1][j+1] - ee77[i-2][j+2];
                if (fabs(ee77[i][j])<1e-5) {
                    ee77[i][j] = 0.0;
                }
            }
            ee77[i][6-i] = temp2*ee77[i-1][6-i+1] - ee77[i-2][6-i+2];
            if (fabs(ee77[i][6-i])<1e-5) {
                ee77[i][6-i] = 0.0;
            }

            if ((ee77[i][0] == 0.0) && (ee77[i][1] == 0.0) && (ee77[i][2] == 0.0) &&
                (ee77[i][3] == 0.0) && (ee77[i][4] == 0.0) && (ee77[i][5] == 0.0) &&
                (ee77[i][6] == 0.0)) {
                break;
            }
        }
        if (fabs(ee77[5][0])>1e-5) {
        //determine f6
            temp1 = ee77[4][0] / ee77[5][0];
            temp2 = (ee77[4][1] - temp1*ee77[5][1]) / ee77[5][0];
            ee77[6][0] = temp2*ee77[5][1] - ee77[4][2];
            if (fabs(ee77[6][0])<1e-5) {
                ee77[6][0] = 0.0;
            }
        }
        else {
            ee77[6][0] = 0.0;
        }

        for (int j=0; j<7; j++) {
            fun_a[0] = fun_a[0] + ee77[0][j]*pow(a, (6-j));
            fun_b[0] = fun_b[0] + ee77[0][j]*pow(b, (6-j));
        }
        for (int i=1; i<7; i++) {
            for (int j=0; j<(7-i); j++) {
                fun_a[i] = fun_a[i] + ee77[i][j]*pow(a, ((6-i)-j));
                fun_b[i] = fun_b[i] + ee77[i][j]*pow(b, ((6-i)-j));
            }
            if ((fun_a[i]*fun_a[i-1]) < 0) {
                change_num_a ++;
            }
            if ((fun_b[i]*fun_b[i-1]) < 0) {
                change_num_b ++;
            }
        }
        num = abs(change_num_a - change_num_b);
        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)------------6
        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 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;
}

输出结果如下：

        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

求得一元六次方程在区间内不相等的实根的总个数之后，怎么具体求出区间内的这些所有的不相等实根呢？

下一章节就干这事~~~~