在“问题五十九”中解一元六次方程时，就有提到，对应的方法可以推广到更高次多项式方程。
现在，由于要用ray tracing画sphere sweeping图形，若对应的曲线是3次b-spline曲线，则需要解一元十次方程（10=2*(2*3-1)）。
我们接下来就看看一元十次方程怎么解？
62.1 理论分析
参考“问题五十九：怎么求一元六次方程在区间内的所有不相等的实根（1）”
http://blog.csdn.net/libing_zeng/article/details/54605140
参考“问题五十九：怎么求一元六次方程在区间内的所有不相等的实根（2）”
http://blog.csdn.net/libing_zeng/article/details/54606500
参考“问题五十九：怎么求一元六次方程在区间内的所有不相等的实根（3）——修正一个问题”
http://blog.csdn.net/libing_zeng/article/details/54619148
一元六次方程和一元十次方程的求解过程有什么不同的地方呢？简单来说，只要将其中出现“六”的地方改成“十”就可以了。主要要改的地方在判断方程在区间内不等实根的总个数的过程中；牛顿迭代在子区间找根的过程几乎完全一样。
62.2 看C++代码实现
代码也几乎和求解一元六次方程时的一样，所以接下来不会有特别注释。
----------------------------------------------vec3.h ------------------------------------------
vec3.h
bool roots_num_equation_10th(float ee10[11], double a, double b, int &num);
bool roots_equation_10th(float ee10[11], float a, float b, float tol, float (&roots)[11]);
----------------------------------------------vec3.cpp ------------------------------------------
vec3.cpp
bool roots_num_equation_10th(float ee10[11], double a, double b, int &num) {
double temp1, temp2;
double ee1111[11][11] = {0.0};
double fun_a[11] = {0.0};
double fun_b[11] = {0.0};
int change_num_a = 0;
int change_num_b = 0;
for (int i=0; i<10; i++) {
//determine f0, f1
ee1111[0][i] = ee10[i];
ee1111[1][i] = (10-i)*ee10[i];
}
ee1111[0][10] = ee10[10];
for (int i=2; i<10; i++) {
//determine f2, f3, f4, f5, f6, f7, f8, f9
temp1 = ee1111[i-2][0] / ee1111[i-1][0];
temp2 = (ee1111[i-2][1] - temp1*ee1111[i-1][1]) / ee1111[i-1][0];
for (int j=0; j<(10-i); j++) {
ee1111[i][j] = temp1*ee1111[i-1][j+2] + temp2*ee1111[i-1][j+1] - ee1111[i-2][j+2];
if (fabs(ee1111[i][j])<1e-5) {
ee1111[i][j] = 0.0;
}
}
ee1111[i][10-i] = temp2*ee1111[i-1][10-i+1] - ee1111[i-2][10-i+2];
if (fabs(ee1111[i][10-i])<1e-5) {
ee1111[i][10-i] = 0.0;
}
if ((ee1111[i][0] == 0.0) && (ee1111[i][1] == 0.0) && (ee1111[i][2] == 0.0) &&
(ee1111[i][3] == 0.0) && (ee1111[i][4] == 0.0) && (ee1111[i][5] == 0.0) &&
(ee1111[i][6] == 0.0) && (ee1111[i][7] == 0.0) && (ee1111[i][8] == 0.0) &&
(ee1111[i][9] == 0.0) && (ee1111[i][10] == 0.0)) {
break;
}
}
if (fabs(ee1111[9][0])>1e-5) {
//determine f10
temp1 = ee1111[8][0] / ee1111[9][0];
temp2 = (ee1111[8][1] - temp1*ee1111[9][1]) / ee1111[9][0];
ee1111[10][0] = temp2*ee1111[9][1] - ee1111[8][2];
if (fabs(ee1111[10][0])<1e-5) {
ee1111[10][0] = 0.0;
}
}
else {
ee1111[10][0] = 0.0;
}
for (int j=0; j<11; j++) {
fun_a[0] = fun_a[0] + ee1111[0][j]*pow(a, (6-j));
fun_b[0] = fun_b[0] + ee1111[0][j]*pow(b, (6-j));
}
for (int i=1; i<11; i++) {
for (int j=0; j<(11-i); j++) {
fun_a[i] = fun_a[i] + ee1111[i][j]*pow(a, ((6-i)-j));
fun_b[i] = fun_b[i] + ee1111[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;
}
bool get_equation_10th_function_and_derivative(float ee10[11], 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<11; i++) {
fnc = fnc + ee10[i]*pow(xxx, (10-i));
if (i<10) {
drv = drv + (10-i)*ee10[i]*pow(xxx, (9-i));
}
}
return true;
}
bool get_root_by_newton_iteration_for_equation_10th(float ee10[11], 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;
for (int i=0; i<2; i++) {
t_k = x0[i];
for (int k=0; k<20; k++) {
if (!(isnan(t_k))) {
get_equation_10th_function_and_derivative(ee10, 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_10th(float ee10[11], float a, float b, float tol, float (&roots)[11]) {
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_10th(ee10, 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_10th(ee10, 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_10th(ee10, 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_10th(ee10, left, right, interval_num);
}
}
return true;
}
----------------------------------------------main.cpp ------------------------------------------
main.cpp
int main(){
//equation of 10th degree
float ee10[11] = {1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};//x^10------------1
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^10=0" << 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(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;
}
}
测试结果如下：
接下来再测试一下如下方程：
// float ee10[11] = {1.0, -1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};//x^9*(x-1)------------2
// float ee10[11] = {1.0, 0.0, -1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};//x^8*(x+1)*(x-1)------------3
// float ee10[11] = {1.0, 0.0, -5.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};//x^6*(x+1)*(x-1)*(x+2)*(x-2)------------5
// float ee10[11] = {1.0, 0.0, -14.0, 0.0, 49.0, 0.0, -36.0, 0.0, 0.0, 0.0, 0.0};//x^4*(x+1)*(x-1)*(x+2)*(x-2)*(x+3)*(x-3)------------7
// float ee10[11] = {1.0, 0.0, -30.0, 0.0, 273.0, 0.0, -820.0, 0.0, 576.0, 0.0, 0.0};//x^2*(x+1)*(x-1)*(x+2)*(x-2)*(x+3)*(x-3)*(x+4)*(x-4)------------9
// 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
测试结果依次如下：
截至当前只测试了这几个方程，后续会将这个程序用于sphere sphere的ray tracing图形。如果遇到问题需要修正该程序的话，后续再来补充。
扫一扫
444
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
