方程求解中的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
我查看了代码,依然没有找到错误,求解释? 试了试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}}
也是一样。 看了你在11楼的发帖,
9#的方程有2处符号不明确
我认为,这不是符号不明, a-+b, a+-b都表示 a-b
这说明只有你的forcal不支持+ ,- 的单目运算功能
===========================
从MATLAB输出的根只有13个,也看得出原方程是一元13次方程,即实质上是一 个方程 2# 风云剑
3# wayne
谢谢两位的Mathematica结果!
同样是Mathematica结果,但结果竟然不一样,看来是所用数据精度的问题。
不过Forcal只使用double数进行计算,在isolve函数求解时,解-282587的误差为0,而我验证时误差又很大,虽仍然不合情理,但似乎仍然是精度问题了。 看了你在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 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位的精度的形式给出,是没有意义的!!! 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 8# G-Spider
你放心,我这里的误差性质跟你的不一样.
我首先是求出方程的准确解(代数表达),然后分别以保留50,60,70,80位的精度 再换成代数表达的形式,回 代入原方程,所得的结果应该越接近0越好, 我上面给的误差就是这个值 9# wayne
我上面的,x=-0.1是精确解,就代入即为0, 好比你的代入原方程的结果应该越接近0越好,
你将原式两边放大125*10^12, 你截断的结果代入此式误差为误差高达10^6,就可能是我上面的情况。所以保留前50位精度未必不可用。 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}
页:
[1]
2