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

forcal

原帖位置：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)
5. ans =
7. 1.0e+005 *
9. -2.825891403741351
10. -0.000425366260653 + 0.000158382367172i
11. -0.000425366260653 - 0.000158382367172i
12. -0.000087958432851 + 0.000371861708019i
13. -0.000087958432851 - 0.000371861708019i
14. 0.000165990822566 + 0.000401197797141i
15. 0.000165990822566 - 0.000401197797141i
16. 0.000288958204839
17. 0.000415706942420
18. -0.000000005695893
19. 0.000000001799567
20. -0.000000000311790 + 0.000000000248976i
21. -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,
9. -42.31708262 - 9.374009122 I,
10. -0.00003117898620 - 0.00002489764282 I,
11. -10.53709954 - 43.27123834 I,
12. 18.66814213 - 32.72325227 I,
13. 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 我查看了代码，依然没有找到错误，求解释？

风云剑

试了试Mathematica7。 Solve[-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 == 0, g] {{g -> -282587.}, {g -> -42.3171 - 9.37401 I}, {g -> -42.3171 + 9.37401 I}, {g -> -10.5371 - 43.2712 I}, {g -> -10.5371 + 43.2712 I}, {g -> -0.000569589}, {g -> -0.000031179 - 0.0000248976 I}, {g -> -0.000031179 + 0.0000248976 I}, {g -> 0.000179957}, {g -> 18.6681 - 32.7233 I}, {g -> 18.6681 + 32.7233 I}, {g -> 33.6859 - 12.3843 I}, {g -> 33.6859 + 12.3843 I}} 也是一样。

wayne

看了你在11楼的发帖,

9＃的方程有2处符号不明确

我认为,这不是符号不明， a-+b , a+-b 都表示 a-b 这说明只有你的forcal不支持+ , - 的单目运算功能 =========================== 从MATLAB输出的根只有13个，也看得出原方程是一元13次方程，即实质上是一 个方程

forcal

2# 风云剑 3# wayne 谢谢两位的Mathematica结果！ 同样是Mathematica结果，但结果竟然不一样，看来是所用数据精度的问题。 不过Forcal只使用double数进行计算，在isolve函数求解时，解-282587的误差为0，而我验证时误差又很大，虽仍然不合情理，但似乎仍然是精度问题了。

forcal

看了你在11楼的发帖, 我认为,这不是符号不明， a-+b , a+-b 都表示 a-b 这说明只有你的forcal不支持+ , - 的单目运算功能 =========================== 从MATLAB输出的根只有13个，也看得出原方程是一元13次 ... wayne 发表于 2011-1-11 09:47

楼主的问题图片如下： 楼主图片中只使用了+号，但楼主提供的公式中该位置为+-，看来公式中楼主不小心写错了。 由于matlab和maple等对单目运算符支持地很好，故不会检查到这个问题。那么a+-b 将解释为a-b了。这么说支持单目运算符连用也是有缺陷的。 Forcal中支持单目运算符+、-、++、--等，但不允许运算符连用。呵呵。

wayne

对于原方程:

-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=0

按无限精度的方式理解系数,即: 140.793343730037 等同于分数 -140793343730037/1000000000000 那么, 原方程换成整数系数的 形式 , 就是:

-17599167966254625 - 622443485455691250000 g - 8262533063207162500000000 g^2 + 48746612998153000000000000000 g^3 + 107847187500000000000000000000000 g^4 - 2472800000000000000000000000000 g^5 - 518954166666666250000000000 g^6 + 1191410714285713750000000000 g^7 - 29658883928571375000000000 g^8 + 390858134920635000000000 g^9 - 3165773809523812500000 g^10 + 15831818181818125000 g^11 + 15831818181818125000 g^12 + 56024354458042 g^13=0

哇,这个方程太壮观了, 下面我给出实数解的80位精度:

-282587.14037411985834653535951392755378605145645703118704732586315336234780773100, -0.00056958933184318274677167457899523885736948661880132982039073611898774134375555451, 0.00017995667287714467434990428451619304579009200924450499176060681086577897412057017

============================================ 我们对值为 -282587.14037 的根进行精度截取，然后回代到原整系数方程，测试如下: {50, -9.094245556*10^6}, {60, 0.001056582233}, {70, -4.409812818*10^-13}, {80, -8.257927053*10^-23}, {90, 1.380405199*10^-32}, {100, 4.350865223*10^-42} 即,当-282587.14037 这个根保留50位精度,误差高达10^6，保留60位精度,误差达到10^-3，保留70位精度,误差达到10^-13 这说明什么？ 说明 方程的解以低于60位的精度的形式给出，是没有意义的！！！

G-Spider

7# wayne "即,当-282587.14037 这个根保留50位精度,误差高达10^6" 比如： x+10x^2=0 其中一精确解为 x=-0.1 若得到x=-0.105，则误差为-0.00475 扩大10倍 10x+100x^2=0 x=-0.1 若得到x=-0.105，则误差为-0.0475

wayne

8# G-Spider 你放心,我这里的误差性质跟你的不一样. 我首先是求出方程的准确解(代数表达),然后分别以保留50,60,70,80位的精度 再换成代数表达的形式，回 代入原方程，所得的结果应该越接近0越好， 我上面给的误差就是这个值

G-Spider

9# wayne 我上面的，x=-0.1是精确解，就代入即为0， 好比你的代入原方程的结果应该越接近0越好， 你将原式两边放大125*10^12， 你截断的结果代入此式误差为误差高达10^6，就可能是我上面的情况。所以保留前50位精度未必不可用。

wayne

8# G-Spider 我又检查了一遍,我的计算应该不会有问题. ==================================== 总结一下, 关于求解高阶多项式方程的根,我们一定要分析根的精度的敏感性. 准确的说,要尤其关注那些模大于1的根,越大于1,越值得关注. 比如,本题,所有的13个根的精度与回代方程后的误差如下: 图中13个根保留10位精度,依次如下: {-282587.1404, -0.0005695893318, 0.0001799566729, -42.31708262 - 9.37400912 I, -42.31708262 + 9.37400912 I, -10.53709954 - 43.27123834 I, -10.53709954 + 43.27123834 I, -0.00003117898620 - 0.00002489764282 I, -0.00003117898620 + 0.00002489764282 I, 18.66814213 - 32.72325227 I, 18.66814213 + 32.72325227 I, 33.68591029 - 12.38430963 I, 33.68591029 + 12.38430963 I}