前续:问题五十九:怎么求一元六次方程在区间内的所有不相等的实根(2)
我们在画“问题六十”的各种回旋体时,遇到这样的问题:
当“基本曲线”的控制点为:
//8--meet some problems
vec3 ctrl_points[6] = {vec3( 4.0, 4.0, 0.0), vec3( 2.0, 4.0, 0.0),
vec3( 2.0, 1.0, 0.0), vec3(-0.5, 1.0, 0.0),
vec3( 1.5, -3.0, 0.0), vec3( 3.0, 0.0, 0.0)};
会出现:在解某个一元六次方程时卡在那个“不断细分区间直到区间内只有一个实根”的while循环里。
原因是:
出现了两种情况:
情况一(s1):当left和right值非常接近,以至于,每次细分之后还是原来的区间;而且该区间上的实根个数大于1,所以,不断细分。
情况二(s2):在细分区间的过程中,竟然出现细分区间的left=right=b(right value of the whole interval),而且,该区间上没有实根,所以,不断重复。
修改方式:
情况一(s1):直接break跳出while循环。该区间的实根个数大于1,但是Newton迭代可能又找不到实根,所以在迭代完之后,再将right值作为实根。
情况二(s2):直接break跳出while循环。但是该区间上没有实根,所以,跳出循环后就没有必要再往后运行了,所以再break跳出for循环。
修改代码如下方红色字体标注:
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;
int offset = 0;
roots[0] = 0.0;
roots_num_equation_6th(ee6, double(a), double(b), total_num);
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) {
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;
if (left == right) {/*s2: handle the situation that both left and right value reach the right end of the whole interval*/
break;
}
}
else if (fabs(left-right)<tol) {
/*s1: left and right value are are so close that their values even don't change under sub-dividing*/
break;
}
}
if (interval_num == 0) {
/*s2: handle the situation that both left and right value reach the right end of the whole interval*/
break;
}
xxx0[0] = left;
xxx0[1] = right;
get_root_by_newton_iteration_for_equation_6th(ee6, xxx0, tol, root);
if (root[0] != 0.0) {
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)) {
/*s1: left and right value are are so close that their values even don't change under sub-dividing,
the interval_num may be more than one, so we modify "interval_num==1" to "interval_num>=1"*/
roots[0] = roots[0] + 1.0;
roots[int(roots[0])] = left;
left = right;
right = b;
}
else {
offset ++;
right = (left + right)/2;
}
roots_num_equation_6th(ee6, left, right, interval_num);
}
}
return true;
}