论坛 › 数学软件 › 方程求解中的bug，令人费解，求解释

forcal 发表于 2011-1-11 07:04:57

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

原帖位置：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对方程形式没有要求，也不认识是不是多项式，但目前只能求得全部实数解。代码：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;
fcopt::isolve;结果（每行前一个数是解，后一个是误差）：
-282587.1403741199      0.
-5.695893318432271e-004   8.315907065430656e-009
1.79956672877141e-004   -1.078234176138977e-010

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

matlab代码：>>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');
>>p=sym2poly(expr);
>>roots(p)

ans =

1.0e+005 *

-2.825891403741351                     
-0.000425366260653 + 0.000158382367172i
-0.000425366260653 - 0.000158382367172i
-0.000087958432851 + 0.000371861708019i
-0.000087958432851 - 0.000371861708019i
0.000165990822566 + 0.000401197797141i
0.000165990822566 - 0.000401197797141i
0.000288958204839                     
0.000415706942420                     
-0.000000005695893                     
0.000000001799567                     
-0.000000000311790 + 0.000000000248976i
-0.000000000311790 - 0.000000000248976imaple结果：0.0001799566729,
33.68591029 + 12.38430963 I,
18.66814213 + 32.72325227 I,
      -10.53709954 + 43.27123834 I,
      -0.00003117898620 + 0.00002489764282 I,
      -42.31708262 + 9.374009122 I,
    -0.0005695893318, -282587.1404,

      -42.31708262 - 9.374009122 I,
      -0.00003117898620 - 0.00002489764282 I,
      -10.53709954 - 43.27123834 I,
       18.66814213 - 32.72325227 I,
      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
我查看了代码，依然没有找到错误，求解释？

风云剑 发表于 2011-1-11 09:17:05

试了试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 发表于 2011-1-11 09:47:43

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

forcal 发表于 2011-1-11 10:43:02

2# 风云剑
3# wayne
谢谢两位的Mathematica结果！
同样是Mathematica结果，但结果竟然不一样，看来是所用数据精度的问题。

不过Forcal只使用double数进行计算，在isolve函数求解时，解-282587的误差为0，而我验证时误差又很大，虽仍然不合情理，但似乎仍然是精度问题了。

forcal 发表于 2011-1-11 10:56:05

看了你在11楼的发帖,

我认为,这不是符号不明， a-+b, a+-b都表示 a-b
这说明只有你的forcal不支持+ ,- 的单目运算功能
===========================
从MATLAB输出的根只有13个，也看得出原方程是一元13次 ...
wayne 发表于 2011-1-11 09:47 http://bbs.emath.ac.cn/images/common/back.gif
楼主的问题图片如下：


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

Forcal中支持单目运算符+、-、++、--等，但不允许运算符连用。呵呵。

wayne 发表于 2011-1-11 12:30:26

方程的根的精度的敏感性分析

本帖最后由 wayne 于 2011-1-11 13:12 编辑

对于原方程:
-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 发表于 2011-1-11 13:29:19

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 发表于 2011-1-11 14:06:44

8# G-Spider
你放心,我这里的误差性质跟你的不一样.

我首先是求出方程的准确解(代数表达),然后分别以保留50,60,70,80位的精度 再换成代数表达的形式，回 代入原方程，所得的结果应该越接近0越好， 我上面给的误差就是这个值

G-Spider 发表于 2011-1-11 14:37:25

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

wayne 发表于 2011-1-11 14:39:00

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}

查看完整版本: 方程求解中的bug，令人费解，求解释