找回密码
 欢迎注册
查看: 19688|回复: 15

[讨论] 方程求解中的bug,令人费解,求解释

[复制链接]
发表于 2011-1-11 07:04:57 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?欢迎注册

×
原帖位置: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
我查看了代码,依然没有找到错误,求解释?
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 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}}

也是一样。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2011-1-11 09:47:43 | 显示全部楼层
看了你在11楼的发帖,
9#的方程有2处符号不明确

我认为,这不是符号不明, a-+b  , a+-b  都表示 a-b
这说明只有你的forcal不支持+ ,  - 的单目运算功能
===========================
从MATLAB输出的根只有13个,也看得出原方程是一元13次方程,即实质上是一 个方程
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2011-1-11 10:43:02 | 显示全部楼层
2# 风云剑
3# wayne
谢谢两位的Mathematica结果!
同样是Mathematica结果,但结果竟然不一样,看来是所用数据精度的问题。

不过Forcal只使用double数进行计算,在isolve函数求解时,解-282587的误差为0,而我验证时误差又很大,虽仍然不合情理,但似乎仍然是精度问题了。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2011-1-11 10:56:05 | 显示全部楼层
看了你在11楼的发帖,

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

楼主的问题图片如下:
多项式公式.PNG

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

Forcal中支持单目运算符+、-、++、--等,但不允许运算符连用。呵呵。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 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位的精度的形式给出,是没有意义的!!!
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 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
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2011-1-11 14:06:44 | 显示全部楼层
8# G-Spider
你放心,我这里的误差性质跟你的不一样.

我首先是求出方程的准确解(代数表达),然后分别以保留50,60,70,80位的精度 再换成代数表达的形式,回 代入原方程,所得的结果应该越接近0越好, 我上面给的误差就是这个值
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2011-1-11 14:37:25 | 显示全部楼层
9# wayne
我上面的,x=-0.1是精确解,就代入即为0, 好比你的代入原方程的结果应该越接近0越好,
你将原式两边放大125*10^12,   你截断的结果代入此式误差为误差高达10^6,就可能是我上面的情况。所以保留前50位精度未必不可用。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2011-1-11 14:39:00 | 显示全部楼层
8# G-Spider
我又检查了一遍,我的计算应该不会有问题.
====================================

总结一下, 关于求解高阶多项式方程的根,我们一定要分析根的精度的敏感性.
准确的说,要尤其关注那些模大于1的根,越大于1,越值得关注.

比如,本题,所有的13个根的精度与回代方程后的误差如下:
截图01.png

图中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}
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

小黑屋|手机版|数学研发网 ( 苏ICP备07505100号 )

GMT+8, 2024-3-28 21:04 , Processed in 0.049313 second(s), 19 queries .

Powered by Discuz! X3.5

© 2001-2024 Discuz! Team.

快速回复 返回顶部 返回列表