相同代码在各个版本matlab中执行结果差异问题

ChaChing

原帖转载自中文 : Matlab相同的代码在各个版本matlab中执行结果问题
http://www.ilovematlab.cn/viewthread.php?tid=61099&page=1#pid578522

为什么相同的代码, 在不用版本的matlab中执行结果不同? 就是下面的代码：
syms D u
a=int(exp(-u^2),u,1/(2*sqrt(D*0.1)),inf);
b=int(exp(-u^2),u,1/(2*sqrt(D*0.3)),inf);
double(solve((a/b),D))

各版本的结果如下
  在matlab7.0中执行的结果：.1010
  在r2009中执行的结果：0.0259
  在r2009a中执行的结果：0.0259
  在v5.3中执行的结果：0.0099
  r2008a中的执行结果： 0.0101000854526561
  matlab2008A ：   0.0101
  maple：     0.01877898738

现在问题的关键是不同版本之间的执行结果为什么不相同? 这样的执行结果的可信度还高吗?

个人使用过ezplot画下a/b vs D的图, 发现D<0.07几乎已为零! 猜测是否各版本间有关solve的收敛准则有所更动了!?
个人水平有限, 想听听诸位大牛的看法!

lengyunfeng

个人认为你这个等式本身就有问题。据Matlab7的帮助信息“solve(eq,var) solves the equation eq (or eq=0 in the two cases cited above) for the variable var.”，solve((a/b),D),是指计算a/b=0时的变量D。不管等式是否有解，首先要保证a/b有意义，也就是说b不能等于0，那也就是说a=0｜b~=0这两个等式必须同时成立（即一个方程组）。根据exp(-u^2)图形，a和b都是确定有解的（我相信每个版本的结果也是一样的），但同时a和b之间相差很小，因为它们只是积分下限的细微差异，当积分上限为inf时，这种细微差异几乎可以忽略不计的，也就是说你现在相当于要求b≈a=0|b~=0这个方程组成立。它们本身就是矛盾的，对于不同程序、同一程序不同版本来说出现不同结果是肯定的，程序没进入死循环就不错了。

feynmand

>> a

a =

-(pi^(1/2)*(erf(1/(2*(D/10)^(1/2))) - 1))/2

>> b

b =

-(pi^(1/2)*(erf(1/(2*((3*D)/10)^(1/2))) - 1))/2

以上是2009b中的a和b的结果。注意到这里面有erf函数，也就是一个积分。
其中erf函数注释有以下内容
%   $Revision: 5.13.4.5 $  $Date: 2009/03/16 22:18:30 $
solve函数注释有以下内容
%   $Revision: 1.1.6.3 $  $Date: 2009/05/07 19:13:57 $
另外涉及到的积分函数最近版本应该也有更新，具体的改动我不清楚还是等别人来补充吧。
这些改动综合到一起会影响到最终的计算结果。我相信这些算法内部的修改在多数情况下会改善计算正确性和计算速度。多数情况下我认为还是更新版本的计算结果更加可信。不过我觉得或许某些特定情况下，算法的改动正好增大了误差，因为不了解具体的内部算法，所以这一点并不能肯定。

针对于本题目来讲，double(solve((a/b),D))完全可以由double(solve(a,D))代替，因为方程等于0必然是分子等于0。而代替以后得到的结果并不一致：
>> double(solve((a),D))

ans =

   0.025633039359632
这是因为计算误差的原因。
实际上exp(-u^2)在u比较大的情况下趋近于0，比如D小于0.07时，exp(-(1/(2*sqrt(D*0.1)))^2)=3.086617322462758e-016，已经小于了eps，所以电脑在计算的时候此时积分式已经等于0了，因此求等式为零成立的条件时难免会产生误差。solve在解方程时，在众多的解中选取了0.0259为最终解，我想这应该是计算时的步长的选取和以前版本有所不同，这个步长的不同或许是由于积分的计算方法不同而造成的，而收敛准则应该没有大的变化，因为是否为最终解只要判断是否是等于0就可以了，不用很复杂的其他方法的。


关于可信度问题
对于本题目来说，实际上是在一个有很大相对误差范围的解区域内求解，得到的结果各版本不同是难免的，如果你想要得到更接近0.07方向的解，那么就需要自己来编制相应的程序，减小步长，从指定方向逼近，这样可能会比较理想。这一过程中需要仔细研究误差传递的理论，尽量减少计算过程中的舍入误差等。


--------
个人使用过ezplot画下a/b vs D的图, 发现D<0.07几乎已为零!
-------
对于计算机来说，D<0.07的情况下已经为0了。

messenger

.
我觉得是误差问题。

这个问题其实就是一个类似於2x=x这样的同解方程，其解为零。

但是 lz 将问题复杂化，将同解方程做了误差函数、半无限区间上的积分运算，最重要的是又做了除法运算。

D本来的值就在零附近，做了除法运算使其值很大，导致微小的变化就会产生很大的误差。

所以这个问题，运算精度越高就会越接近真实解零。

试了一下将精度分别设为30、90、150，精度越高越接近零解，但越往后接近的越慢。

>> syms D u
>> aa=vpa(int(exp(-u^2),u,1/(2*sqrt(D*0.1)),inf),150);
>> bb=vpa(int(exp(-u^2),u,1/(2*sqrt(D*0.3)),inf),150);
>> solve(aa/bb)

至於为什么不同的版本结果会不一样，我觉得可能是因为 inf 取值不一样，因为有近似除零运算，使 inf 的微小区别，都会带来很大误差。
                       .

feynmand

个人认为你这个等式本身就有问题。据Matlab7的帮助信息“solve(eq,var) solves the equation eq (or eq=0 in the two cases cited above) for the variable var.”，solve((a/b),D),是指计算a/b=0时的变量D。不管等式 ...
lengyunfeng 发表于 2009-12-28 17:01

被2楼提醒了，这个题目中b是一个很小的数，作为分母的话会将误差成倍的放大，在计算中应该避免这种表达方式。程序没有陷入死循环还是很正常的，只是过程中将相对误差放大了。

lengyunfeng

5# feynmand
谢谢指教。我晚上想了一会儿，觉得这种题如果以纸面的形式出给我做的话，我肯定要首先考虑这个式子本身是否有意义。Matlab作为数学科学运用软件的巨头之一，难道他们在编写源代码的时候会忽略这个问题？有点不可思议～～～不过既然Matlab软件对它有确定结果显示，既不报错也不警告，的确有点费解。。。。

feynmand

5# feynmand  
谢谢指教。我晚上想了一会儿，觉得这种题如果以纸面的形式出给我做的话，我肯定要首先考虑这个式子本身是否有意义。Matlab作为数学科学运用软件的巨头之一，难道他们在编写源代码的时候会忽略这个问 ...
lengyunfeng 发表于 2009-12-28 19:31

作为一个通用软件，编写算法程序的时候考虑的肯定是这一算法是否能够覆盖比较多的情况。那么在方程求解或者积分的时候肯定不能够针对某一特定情况给出最优的算法。从竖直解法角度来说本例中matlab给出了一个可行的解，并没有错误，也用不着发出警告。

一般而言数值解给出的都是特解，针对多解的情况软件只给出一个那么就是一个成功的软件了。

sogeno

请教版主，你说的求解精度设置是指设置小数点后的位数吗？本人疑惑中