来自MATLAB的python新手。在
我用的是一个震级函数的双曲正切截断。
我在将0.5 * math.tanh(r/rE-r0) + 0.5
函数应用于范围值r = np.arange(0.1,100.01,0.01)
的数组时遇到了问题。我在接近0的一侧得到了函数的几个0.0
值,当我执行对数时,这些值会导致域问题:
P1 = [ (0.5*m.tanh(x / rE + r0 ) + 0.5) for x in r] # truncation function
我使用以下方法:
^{pr2}$这对我现在所做的已经足够了,但有点像创可贴。在
根据数学明确性的要求:
在天文学中,震级的工作原理大致如下:
mu = -2.5log(flux) + mzp # apparent magnitude
式中,mzp是每秒看到1个光子的大小。因此,较大的通量等于较小(或更负)的表观震级。我正在为使用多个组件函数的源创建模型。例如,两个具有不同sersic索引的sersic函数,其内部组件上有一个P1
外部截断,外部组件上有一个1-P1
内部截断。这样,当向每个分量添加截断函数时,由radius定义的大小将变得非常大,因为当P1
渐近接近零时,mu1-2.5*log(P1
)会变得非常小。在
TLDR:我想知道的是,是否有一种方法可以保存精度不足以与零区分的浮点(特别是在渐近接近零的函数的结果中)。这一点很重要,因为当取这些数的对数时,结果是一个域错误。在
非对数P1中的输出开始读取零之前的最后一个数字是5.551115123125783e-17
,这是一个常见的浮点算术舍入错误结果,其中所需的结果应为零。在
如有任何意见,我们将不胜感激。在
在@用户:丹 不把我的整个剧本放进去:
xc1,yc1 = 103.5150,102.5461;
Ee1 = 23.6781;
re1 = 10.0728*0.187;
n1 = 4.0234;
# radial brightness profile (magnitudes -- really surface brightness but fine in ex.)
mu1 = [ Ee1 + 2.5/m.log(10)*bn(n1)*((x/re1)**(1.0/n1) - 1) for x in r];
# outer truncation
rb1 = 8.0121
drs1 = 11.4792
P1 = [ (0.5*m.tanh( (2.0 - B(rb1,drs1) ) * x / rb1 + B(rb1,drs1) ) + 0.5) for x in r]
P1 = [ -2.5*m.log10(x) if x!=0.0 else np.inf for x in P1 ] # band-aid for problem
mu1t = [x+y for x,y in zip(P1,mu1)] # m1 truncated by P1
式中,bn(n1)=7.72,B(rb1,drs1)=2.65-4.98*(r_b1/(-drs1))
mu1是要截断的分量的幅度分布。P1是截断函数。P1的许多最终条目都是零,这是由于浮点精度的原因,浮点与零没有区别。在
了解问题的简单方法:
>>> r = np.arange(0,101,1)
>>> P1 = [0.5*m.tanh(-x)+0.5 for x in r]
>>> P1
[0.5, 0.11920292202211757, 0.01798620996209155, 0.002472623156634768, 0.000335350130466483, 4.539786870244589e-05, 6.144174602207286e-06, 8.315280276560699e-07, 1.1253516207787584e-07, 1.5229979499764568e-08, 2.0611536366565986e-09, 2.789468100949932e-10, 3.775135759553905e-11, 5.109079825871277e-12, 6.914468997365475e-13, 9.35918009759007e-14, 1.2656542480726785e-14, 1.7208456881689926e-15, 2.220446049250313e-16, 5.551115123125783e-17, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
还要注意零之前的浮动。在
有两件事你可以试试。
(1)蛮力法:找到一个可变精度浮点运算包,用它代替内置的固定精度。我在Maxima[1]中讨论你的问题,我发现为了避免下溢,我不得不提高浮点精度,但这是可能的。如果你愿意的话,我可以发布Maxima代码。我可以想象有一些适合Python的可变精度浮点库。
(2)用泰勒级数或其它近似方法近似对数((1/2)(1+tanh(-x)),以避免完全对数(tanh(…))。
[1]http://maxima.sourceforge.net
回想一下双曲正切可以表示为(1-e^{-2x})/(1+e^{-2x})。用一点代数,我们可以得到0.5*tanh(x)-0.5(函数的负数)等于e^{-2x}/(1+e^{-2x})。它的对数是
-2*x-log(1+exp(-2*x))
,它在任何地方都是有效的和稳定的。也就是说,我建议您更换:
通过这种更简单、更稳定的方法:
^{pr2}$相关问题 更多 >
编程相关推荐