问题六十五：二叉查找树的一个应用实例——求解一元十次方程时单实根区间的划分

65.1 概述

回忆一下：

“问题五十九：怎么求一元六次方程在区间内的所有不相等的实根”和“问题六十二：怎么求一元十次方程在区间内的所有不相等的实根”中求一元六次方程和一元十次方程的实根时，我们需要将求根区间划分成若干个单根子区间，然后再用牛顿迭代法求出方程在子区间的该单实根。

关于“划分单根子区间”，我们之前的做法是：

1，若子区间内实根的个数大于1，则将区间二分（将区间的右端值移至区间中点）；

2，若子区间内实根的个数小于1（没有实根），则将区间的左端值移至当前区间的右端值，右端值移至初始区间的右端值；

 

该算法的一种可能的区间划分方式示意图如下：

“问题六十三”中最后一张图形生成的时间大概为680s（因计算机负载情况而异），贴图如下：

接下来，我们用二叉查找树的结点来表示区间的划分过程：

1，根结点为初始区间；

2，若区间内实根的个数大于1，则将区间二分为左右两个子区间，分别用二叉查找树的左右孩子结点表示；

3，最终得到的二叉查找树的叶子结点只有两种情况：单实根区间和没有实根的区间；

4，所以，只需要遍历二叉查找树的有实根的叶子结点即可得到最终的迭代区间；

65.2 看C++代码实现

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

vec3.h

bool roots_equation_10th_bt(float ee10[11], float a, float b, float tol, float (&roots)[11]);

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

vec3.cpp

/*表示子区间的结构体*/
struct IntervalStruct
{
    double key;//成员key，可有可无。程序中赋值子区间中点值。
    double leftValue;//子区间左端点值
    double rightValue;//子区间右端值
    int rootNum;//子区间内实根的个数
};
typedef IntervalStruct ItemType;//给子区间定义一个别名
struct IntervalTreeNode//定义二叉查找树的结点
{
    ItemType info;//data
    IntervalTreeNode* left;//左孩子结点的指针
    IntervalTreeNode* right; //右孩子结点的指针
};

/*划分单实根区间，创建二叉查找树*/
void CreateIntervalTree(float ee10[11], double a, double b, int num, IntervalTreeNode *rootNode) {
    double left = a;
    double right = b;
    int leftNum, rightNum;
    rootNode->info.key = (left+right)/2;
    rootNode->info.leftValue = left;
    rootNode->info.rightValue = right;
    rootNode->info.rootNum = num;
    rootNode->left = NULL;
    rootNode->right = NULL;
if ((num > 1) && (fabs(left-right)>1e-6)) {
/*该算法忽略存在实根的小区间，所有会有漏根情况出现。之前的算法是将小区间的一个端点值直接作为一个实根（貌似更为合理）*/
        IntervalTreeNode *leftNode = new IntervalTreeNode;
        rootNode->left = leftNode;
        left = rootNode->info.leftValue;
        right = (rootNode->info.leftValue + rootNode->info.rightValue) / 2;
        roots_num_equation_10th(ee10, left, right, leftNum);
        CreateIntervalTree(ee10, left, right, leftNum, leftNode);//递归左子树

        IntervalTreeNode *rightNode = new IntervalTreeNode;
        rootNode->right = rightNode;
        left = (rootNode->info.leftValue + rootNode->info.rightValue) / 2;
        right = rootNode->info.rightValue;
//        roots_num_equation_10th(ee10, left, right, rightNum);
        rightNum = num-leftNum;
        CreateIntervalTree(ee10, left, right, rightNum, rightNode); //递归右子树
    }
}

/*遍历有实根的叶子结点。函数返回的是所有“有实根的叶子结点”表示的区间*/
void TraverseLeavesForValidIntervals(IntervalTreeNode *rootNode, int &counter, ItemType (&leaves)[11]) {
    if(rootNode!=NULL) {
        TraverseLeavesForValidIntervals(rootNode->left, counter, leaves);
        if ((rootNode->left == NULL) && (rootNode->right == NULL) && (rootNode->info.rootNum != 0)) {
            if (counter == 0) {counter = 1;}
            leaves[counter].leftValue = rootNode->info.leftValue;
            leaves[counter].rightValue = rootNode->info.rightValue;
            leaves[counter].rootNum = rootNode->info.rootNum;
            leaves[0].rootNum = counter;
            std::cout << "[" << leaves[counter].leftValue << ", " << leaves[counter].rightValue << "]: rootNum " << leaves[counter].rootNum << endl;/*测试时用这句代码打印区间信息；实际画图时不需要这句代码，以免log太多*/
            counter ++;
        }
        TraverseLeavesForValidIntervals(rootNode->right, counter, leaves);
    }
}

/*as to single real root intervals, the interval ends may be not close enough to the root, 
so we use this function to narrow the intervals*/ 
void NarrowValidIntervals(float ee10[11], double a, double b, double &na, double &nb) {
    int num = 0;
    na = a;
    nb = (a+b)/2;
    roots_num_equation_10th(ee10, na, nb, num);
    if (num == 0) {
        na = (a+b)/2;
        nb = b;
    }
}

bool roots_equation_10th_bt(float ee10[11], float a, float b, float tol, float (&roots)[11]) {
    IntervalTreeNode *rootNode = new IntervalTreeNode;
    ItemType intervals[11];
    double root[2], xxx0[2];
    double left = a;
    double right = b;
    double left_v, right_v;
    int totalNum, num;
    bool found;
    num = 0;
    roots[0] = 0.0;
    roots_num_equation_10th(ee10, left, right, totalNum);//计算总的实根个数
    CreateIntervalTree(ee10, left, right, totalNum, rootNode);//创建二叉查找树
    TraverseLeavesForValidIntervals(rootNode, num, intervals);/*遍历有实根的叶子结点（由于叶子结点表示的区间内最多只有一个实根，所以该函数得到的为单实根区间）*/
    for (int i=1; i<num; i++) {
        left = intervals[i].leftValue;
        right = intervals[i].rightValue;
        found = false;
        while (!found) {/*对于每一个单实根区间，如果区间太宽，该算法会在缩小区间后反复迭代，直到找到对应的实根为止*/
            xxx0[0] = left;
            xxx0[1] = right;
            get_root_by_newton_iteration_for_equation_10th(ee10, xxx0, tol, root);
            if (root[0] != 0.0) {
                roots[0] = roots[0] + 1.0;
                roots[int(roots[0])] = root[1];
                found = true;
            }
            else if ((root[0] == 0.0)&&(fabs(left-right)<tol)) {
                roots[0] = roots[0] + 1.0;
                roots[int(roots[0])] = left;
                found = true;
            }
            else {
                NarrowValidIntervals(ee10, left, right, left_v, right_v);//缩小单实根区间
                left = left_v;
                right = right_v;
            }
        }
    }
    return true;
}

65.3 测试算法：求解一元十次方程

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

mian.cpp

    #include <iostream>
    #include <fstream>
    #include <limits>
    #include <stdlib.h>
    #include <vec3.h>

    using namespace std;

    int main(){

        float ee10[11] = {1.0, -5.0, -30.0,  150.0,  273.0, -1365.0,  -820.0, 4100.0,  576.0, -2880.0,  0.0};//x*(x+1)*(x-1)*(x+2)*(x-2)*(x+3)*(x-3)*(x+4)*(x-4)*(x-5)------------10
        float a = -10;
        float b =  10;
        int num;
        float tol = 1e-6;
        float roots[11];
        roots_num_equation_10th(ee10, a, b, num);
        std::cout << "the coefficients of the equation of 10th degree are:" << endl;
        for (int i=0; i<11; i++) {
            std::cout << ee10[i] << "  ";
        }
        std::cout << endl;
        std::cout << "the corresponding equation is:" << endl;
        std::cout << "x^2*(x+1)*(x-1)*(x+2)*(x-2)*(x+3)*(x-3)*(x+4)*(x-4)" << endl;
        std::cout << endl;
        std::cout << "the number of roots of this equation in the interval [" << a << ", " << b << "] is: " << num << endl;

        if (num > 0) {
            roots_equation_10th_bt(ee10, a, b, tol, roots);
//            roots_equation_10th(ee10, a, b, tol, roots);
            std::cout << "the roots are: ";
            for (int j=1; j<(num+1); j++){
                std::cout << roots[j] << " ";
            }
            std::cout << endl;
        }
    }

测试结果截图如下：

65.3 测试算法：图形生成

如本章节开始时提到：

“问题六十三”中最后一张图形生成的时间大概为680s（因计算机负载情况而异）

 

使用该算法，生成该图形所需的时间大概为520s（因计算机负载情况而异）

贴图如下：

65.4 之前算法的另一版本

1，若子区间内实根的个数大于1，则将区间二分（将区间的右端值移至区间中点）；

2，若子区间内实根的个数小于1（没有实根），则将区间的左端值移至当前区间的右端值，右端值移至父区间的右端值；

 

该算法的一种可能的区间划分方式示意图如下：

对应的修改：

roots_equation_10th()函数中left、right的赋值都做如下修改：

//                        left = right;
//                        right = b;
                    right = 2*right - left;
                    left = (left + right) / 2;

这个修改之后，可以减少区间的划分次数，但是最终生成的图形有缺失（先不细究原因），贴图如下：