正数开方咋就出来复数了？Bug？很简单的一道计算

TBE_Legend

位于 c 盘 根目录下 的 data1.txt文件：

1. 200 279 105.60 519.43
   
2. 400 270 65.13 458.25
   
3. 1600 255 76.18 254.00
   
4. 3200 230 100.76 149.45
   
5. 6400 202 114.51 110.53
   
6. 12800 170 117.56 80.29
   
7. 25600 150 114.84 40.82
   
8. 51200 125 108.54 8.46
   
9. 102400 100 90.43 -15.41
   

复制代码

正数开方居然算出复数了？

1. Clear["Global`*"];
   
2. data = Import["C:\\data1.txt", "Data"] // ToExpression;
   
3. {rexdata, imxdata, reydata, imydata} = {data[[All, 1]],
   
4.    data[[All, 2]], data[[All, 3]], data[[All, 4]]};
   
5. f = ((a x^b)/c^x^d /. {x -> rexdata + imxdata I});
   
6. problem = f - (reydata + imydata I);
   
7. min = NMinimize[Total[Re[problem]^2 + Im[problem]^2]^(
   
8.    1/2), {a, b, c, d},
   
9.    Method -> {"RandomSearch", "SearchPoints" -> 10,
   
10.      "RandomSeed" -> 2}] // Timing
    

复制代码

waynebuaa

目标函数有Re，Im操作，于是Mathematica就将目标函数的参数a,b,c,d自动提升为复域，所以，我们有必要告诉Mathematica不要将a,b,c,d的“类型"自动转化为复数”类型“，有点像C语言。

1. min = NMinimize[{Sqrt[Total[Re[problem]^2 + Im[problem]^2]], a^2 > 0 && b^2 > 0 && c^2 > 0 && d^2 > 0}, {a, b, c, d}, Method -> {"RandomSearch", "SearchPoints" -> 10,
   
2.      "RandomSeed" -> 2}] // Timing

复制代码

lin2009

表达式用到了Re[problem]^2 + Im[problem]^2的对复数的操作，所以MM就把这段程序放在复数域进行的。
可以注意到，画圈的部分虚部为0，对结果应该没有影响。
正如0^0 问题一样(如下)，软件计算有时不如人工计算来的简洁。
http://forum.simwe.com/thread-937601-1-1.html

TBE_Legend

表达式用到了Re[problem]^2 + Im[problem]^2的对复数的操作，所以MM就把这段程序放在复数域进行的。
可以注意到，画圈的部分虚部为0，对结果应该没有影响。
正如0^0 问题一样(如下)，软件计算有时不如人工计算来的 ...
lin2009 发表于 2010-6-17 18:13

什么啊？ 我分别把实部虚部提取出来，哪还有什么复数的问题？再说了，我的第二问又该如何解释呢？

lin2009

Re和Im本身是对复数进行操作的，所以整个计算是放在复数域进行，计算结果可能就以实部＋0虚部的形式显示出来，这应该是MM的计算机制，Maple也有这个现象。

TBE_Legend

Re和Im本身是对复数进行操作的，所以整个计算是放在复数域进行，计算结果可能就以实部＋0虚部的形式显示出来，这应该是MM的计算机制，Maple也有这个现象。
lin2009 发表于 2010-6-17 19:06

这样模糊的、似是而非的结论毫无意义！

我都说了，实部虚部我都已经非别提出来了，根本没有什么复数的事情了，问题有可能出现在优化NMinimize上。我还说了，你看看我的第二个图片中的问题！我把出现复数时的a,b,c,d带入后却有没有复数了。

TBE_Legend

希望这个问题能有个清晰的结论或解释。 版主帮忙置顶下吧，等解决了再取消置顶。

lin2009

5# TBE_Legend
楼主没仔细看我的意思吧，可能是我写得还不是很清楚吧，不作解释了。
下面的贴子虽然是关于maple的，但对于MM可能也有所帮助。
关键词： zero imaginary part
出处：http://www.mapleprimes.com/posts/41752-Bug--Evalf-Hypergeometric-2F1
“
The initial problem you've described is a bug which we will fix. Thanks for pointing it out. For the other issue, note that as a general rule, in the floating point domain, Maple closes branch cuts from *both* sides, using the sign of the zero imaginary part to associate "sides". For your example, we thus have: > tstData:=[a=1.2+.4*I, b=-2.6+.85*I, c=-1.0+.7*I, z= 2.8 - 0.0*I]; tstData := [a = 1.2 + 0.4 I, b = -2.6 + 0.85 I, c = -1.0 + 0.7 I, z = 2.8 - 0. I] > TstData:= op(eval([a,b,c,z], (tstData))): > MPL_2F1(TstData); -720.102249494927486 + 440.743419599155560 I > tstData:=[a=1.2+.4*I, b=-2.6+.85*I, c=-1.0+.7*I, z= 2.8 + 0.0*I]; tstData := [a = 1.2 + 0.4 I, b = -2.6 + 0.85 I, c = -1.0 + 0.7 I, z = 2.8 + 0. I] > TstData:= op(eval([a,b,c,z], (tstData))): > MPL_2F1(TstData); 1.49139640953106626 - 3.74804421753558801 I showing the limiting values from below (the first) and above (the second). The question then becomes: If the user is not explicit in specifying the side of the branch cut they want, which side should be assumed? That is, do we associate z=2.8 with z=2.8+0.0*I or z=2.8-0.0*I? There is no "right" answer to this, and the best we can hope to achieve is consistency. To make this decision, Maple uses the convention of "counter-clockwise continuity" around a (finite) branch point, meaning that the closure rule is determined by taking a limit along a path around the branch point in a counter-clockwise direction. For branches extending to positive infinity, this implies the closure is from below, so for a 2F1 such as yours, the association is z=2.8 => z=2.8-0.0*I. Note that this convention determines other things, such as the logarithmic form of inverse trig functions and the identities relating them. As a final note, I am confused by your remarks concerning an apparent difference in this behaviour for 2F1 between Maple 10 and Maple 11. My tests indicate that the same results obtain in both versions. Can you provide an explicit example, with results, showing the discrepancy? Thx. Dave Hare Manager Mathematical Software Development Group Maplesoft
”
"Dave, Thank you for the detailled answer, it is very welcome. Certainly it would be worth to add such to the docu, especially for the numerics ... I was not aware of this rule in Maple. And also not of the subtility of the 'complex zero'. Please find attached the version for M10 and M11, printed as pdf. There seems to be a difference in handling 2.8 + 1e-300*I*0 giving 2.8 or 2.8 + 0.0*I (and after that it is covered by what you said)."

TBE_Legend

5# TBE_Legend
楼主没仔细看我的意思吧，可能是我写得还不是很清楚吧，不作解释了。
下面的贴子虽然是关于maple的，但对于MM可能也有所帮助。
关键词： zero imaginary part
出处：http://www.mapleprimes.com/ ...
lin2009 发表于 2010-6-18 10:08

多谢！看了下，还是想看看直接的mmtc用户或官方的解释。

waynebuaa

当Mathematica报错的时候，出现不可理解的结果是可以理解的，就怕Mathematica不报错却出现不可理解的结果了，当然了，这个时候也就肯定是bug了。

如果你对Mathematica报错的时候返回的结果还信任，只是不喜欢那个几乎为0的虚部的话，可以试试Chop命令
========================================================
另外，最好多看看出错原因吧：

Typically this occurs when the function is not smooth or the values of the function are much less accurate than the WorkingPrecision

waynebuaa

BTW，:victory:：

Total[Re[problem]^2 + Im[problem]^2]^ (1/2)

其实就是求problem的模，你可以直接用Norm函数

waynebuaa

如果你信任Mathematica报错的时候给的结果的话，可以这么来写：

1. 
   
2. min = Chop@Quiet@NMinimize[Norm[problem], {a, b, c, d}, Method -> {"RandomSearch", "SearchPoints" -> 10,  "RandomSeed" -> 2}] // Timing

复制代码

waynebuaa

1# TBE_Legend
另外，前面读取数据部分的代码可以这样来：

1. {rexdata, imxdata, reydata, imydata} =  Import["/data1.txt", "Data"] // Transpose

复制代码

TBE_Legend

复数的出现主要是因为用了开方。

这个问题是我做优化问题时偶然发现的，如果不能开方，可以找到最优解，但如果开方的话，则会出现上面的警告而且找不到最优值了。

TBE_Legend

而且我试着把原始数据都用准确表示（如，105.60用 10560/100），优化时我也把精度调大或调小，都无济于事。

所以我认为，问题不在于精度，而在于那个开方运算。

TBE_Legend

问题总算是得到了一个满意的解决，非常感谢 lin2009，wayne 很有建设性的讨论，受益良多，有机会来哈尔滨的话一定要来找我 :handshake

waynebuaa

这么快就入选为最佳答案了，俺还没回过神来呢，呵呵，其实lin2009已经说出问题的根源了。
俺只是推测，瞎试出来的