方程求解中的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等对非线性方程全部解求解能力如何？看其隐函数绘图能力很强大，里面必然包含了对非线性方程全部解的求解功能。