方程求解中的bug，令人费解，求解释

wanglu

原帖位置：http://www.matlabsky.com/forum-viewthread-tid-3789-extra-%26page%3D1-page-1.html

问题描述：求以下方程的全部解：
f(g)=-140.793343730037-4979547.88364553*g-66100264505.6573*g^2+389972903985224*g^3+862777500000000000*g^4-19782400000000000*g^5-4151633333333.33*g^6+9531285714285.71*g^7-237271071428.571*g^8+3126865079.36508*g^9-25326190.4761905*g^10+126654.545454545*g^11+126654.545454545*g^12+0.448194835664336*g^13;

这是个多项式方程，Forcal对方程形式没有要求，也不认识是不是多项式，但目前只能求得全部实数解。代码：

1. f(g)=-140.793343730037-4979547.88364553*g-66100264505.6573*g^2+389972903985224*g^3+862777500000000000*g^4-19782400000000000*g^5-4151633333333.33*g^6+9531285714285.71*g^7-237271071428.571*g^8+3126865079.36508*g^9-25326190.4761905*g^10+126654.545454545*g^11+126654.545454545*g^12+0.448194835664336*g^13;
   
2. fcopt::isolve[HFor("f")];

复制代码

结果（每行前一个数是解，后一个是误差）：
-282587.1403741199        0.
-5.695893318432271e-004   8.315907065430656e-009
1.79956672877141e-004     -1.078234176138977e-010

matlab和maple可用多项式解法求得复数域内的全部解。

matlab代码：

1. >>expr=sym('-140.793343730037-4979547.88364553*g-66100264505.6573*g^2+389972903985224*g^3+862777500000000000*g^4-19782400000000000*g^5-4151633333333.33*g^6+-9531285714285.71*g^7-237271071428.571*g^8+3126865079.36508*g^9-25326190.4761905*g^10+-126654.545454545*g^11+126654.545454545*g^12+0.448194835664336*g^13');
   
2. >>p=sym2poly(expr);
   
3. >>roots(p)
   
4. ans =
   
5.   1.0e+005 *
   
6. -2.825891403741351                     
   
7. -0.000425366260653 + 0.000158382367172i
   
8. -0.000425366260653 - 0.000158382367172i
   
9. -0.000087958432851 + 0.000371861708019i
   
10. -0.000087958432851 - 0.000371861708019i
    
11.   0.000165990822566 + 0.000401197797141i
    
12.   0.000165990822566 - 0.000401197797141i
    
13.   0.000288958204839                     
    
14.   0.000415706942420                     
    
15. -0.000000005695893                     
    
16.   0.000000001799567                     
    
17. -0.000000000311790 + 0.000000000248976i
    
18. -0.000000000311790 - 0.000000000248976i

复制代码

maple结果：

1.   0.0001799566729,
   
2.   33.68591029 + 12.38430963 I,
   
3.   18.66814213 + 32.72325227 I,
   
4.       -10.53709954 + 43.27123834 I,
   
5.         -0.00003117898620 + 0.00002489764282 I,
   
6.         -42.31708262 + 9.374009122 I,
   
7.     -0.0005695893318, -282587.1404,
   
8.         -42.31708262 - 9.374009122 I,
   
9.         -0.00003117898620 - 0.00002489764282 I,
   
10.         -10.53709954 - 43.27123834 I,
    
11.        18.66814213 - 32.72325227 I,
    
12.         33.68591029 - 12.38430963 I

复制代码

验证了一下matlab和maple的实根，复根没有验证。
-2.825891403741351        误差大，不是正解，与Forcal和maple结果一致
  0.000288958204839       误差大，不是正解
  0.000415706942420       误差大，不是正解
-0.000000005695893        误差小，是正解，与Forcal和maple结果一致
  0.000000001799567       误差小，是正解，与Forcal和maple结果一致

令人费解的是，matlab、maple和Forcal用不同的算法都求得了一个错误的解：-2.825891403741351
我查看了代码，依然没有找到错误，求解释？

qibbxxt

是不是可以从帮助文件中

1. Algorithm
   
2. 
   
3. The algorithm simply involves computing the eigenvalues of the companion matrix:
   
4. 
   
5. A = diag(ones(n-1,1),-1);
   
6. A(1,:) = -c(2:n+1)./c(1);
   
7. eig(A)
   
8. 
   
9. It is possible to prove that the results produced are the exact eigenvalues of a matrix within roundoff error of the companion matrix A, but this does not mean that they are the exact roots of a polynomial with coefficients within roundoff error of those in c.

复制代码

去找找解释呢？

wanglu

谢谢qibbxxt ！
这个帖子中也有一些解释，可参考一下：http://bbs.emath.ac.cn/thread-2853-1-1.html

bainhome

所给链接说明为正解，但这种系数极度震荡的多项式一般十分少见，感觉如果不是专项搞数值分析算法研究的人，没必要深究，里面黑幕重重，实在是有害身心健康。
即使您是个不幸的人，无可避免遇到这种害虫式的方程，我想多数人会具体专业情况下做具体分析，比如给出近似方程，比如由所求内容具体专业给出系数所在范围用数值方法寻找部分合理的根，然后到此为止。
本题中最后四个根(含最后一组共轭根)感觉都还比较合适。
个人见解，未必准确。

feiyuzhen

在maple中用命令Student[Calculus1][Roots]试试

yuehua1985

学习了 多谢多谢额