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

forcal

对Forcal的结果，我写了个验证程序：先将isolve的求解结果存放到数组a中，输出结果（17位有效数字），然后将求得的解带入函数f求值，再输出结果（17位有效数字）。代码：

1. !using("fcopt","math");
2. 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;
3. main(:a)=
4. {
5. oo{a=array[6]},
6. isolve[HFor("f"),optarray,a],
7. a.outm(17),
8. a[1]=f[a(0)], a[3]=f[a(2)], a[5]=f[a(4)],
9. a.outm(17)
10. };

复制代码

结果： -282587.14037411986 0. -5.6958933184322706e-004 8.3159070654306564e-009 1.7995667287714099e-004 -1.0782341761389774e-010 -282587.14037411986 0. -5.6958933184322706e-004 8.3159070654306564e-009 1.7995667287714099e-004 -1.0782341761389774e-010 可见所得结果相同，那么问题出在哪儿呢？ 我将本次输出的-282587.14037411986和上次验证所用的-282587.1403741199 再次带入函数f验证，代码：

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. f[-282587.14037411986];
3. f[-282587.1403741199];

复制代码

结果： 0. -1.225996432692711e+055

mathe

x=[-282587.14037411986 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; 代入，那么最高项的值大概为0.45*282587^13=3.3*10^78 所以1*10^55相对它来说是0，正常误差范围

wayne

9# wayne 我上面的，x=-0.1是精确解，就代入即为0， 好比你的代入原方程的结果应该越接近0越好， 你将原式两边放大125*10^12， 你截断的结果代入此式误差为误差高达10^6，就可能是我上面的情况。所以保留前50位 ... G-Spider 发表于 2011-1-11 14:37

我算出来的那个误差就是楼主forcal给的函数 f(g) 在 g等于算出来的各种精度的根处的取值，并没有扩大，

wayne

我还是把我处理的过程详细的再过一遍吧， 我假设大家熟悉Mathematica语言： 一、把楼主给的方程用代数形式表达出来： 直接把方程复制到Mathematica是可以的，只是Mathematica默认会把系数理解成浮点数的形式，因此我们得先在每一个系数后面加一个NumberMark，告诉Mathematica这些是高精度的数，然后，才能 使用Rationalize函数进行有理化，代码如下：

1. f=Rationalize[-140.793343730037\`100 - 4979547.88364553\`100*g -
2. 66100264505.6573\`100*g^2 + 389972903985224\`100*g^3 +
3. 862777500000000000\`100*g^4 - 19782400000000000\`100*g^5 -
4. 4151633333333.33\`100*g^6 + 9531285714285.71\`100*g^7 -
5. 237271071428.571\`100*g^8 + 3126865079.36508\`100*g^9 -
6. 25326190.4761905\`100*g^10 + 126654.545454545\`100*g^11 +
7. 126654.545454545\`100*g^12 + 0.448194835664336\`100*g^13]

复制代码

二、解方程 对于Mathematica来说，这是再平凡不过的方程了，系数再怎么的大，只要你的机器硬件可以忍，它都可以忍， 代码：

1. ans = Solve[f == 0, g];

复制代码

我在此加了分号，目的是不显示结果，因为这个结果是代数的形式表达的，太繁杂，没啥好看的。如果你好奇，可以去掉这个分号。如果你不放心这结果的正确性，可以用求值函数显示，比如我要看看这13个根的20位精度是个什么样子，代码：g /. N[ans, 20] 三、把上面的13个根按各种精度截取，回代入原方程，计算出误差

1. Table[{ii, SetPrecision[f /. Rationalize[SetPrecision[N[ans, ii], 1000]], 10]}, {ii, 10, 100, 10}] // Grid

复制代码

注意,按精度截取后,必须转换成代数的形式才回代进去,这样做能确保误差计算的高度精确. 输出会有点乱，末尾加的grid是为了格式化输出，如果你对那些误差很小，基本上都是10的-7 以下的那种数可以忽视的话，可以使用Chop函数，Chop[%]，Chop函数只对很小很小的数感兴趣，它的作用是将很小很小的数置换成0 附件是nb文件： 5.nb (36.92 KB, 下载次数: 0)

2011-1-12 13:07 上传

点击文件名下载附件

forcal

感谢大家的讨论，特别wayne 的总结非常好！ 还想了解一下，Mathematica、maple、matlab等对非线性方程全部解求解能力如何？看其隐函数绘图能力很强大，里面必然包含了对非线性方程全部解的求解功能。