6

来自 MATLAB 的 Python 新手。

我正在使用幅度尺度函数的双曲正切截断。0.5 * math.tanh(r/rE-r0) + 0.5将函数应用于范围值数组时遇到问题r = np.arange(0.1,100.01,0.01)0.0我在接近零的一侧得到了几个函数值,这在我执行对数时会导致域问题:

P1 = [ (0.5*m.tanh(x / rE + r0 ) + 0.5) for x in r] # truncation function

我使用这个解决方法:

P1 = [ -m.log10(x) if x!=0.0 else np.inf for x in P1 ]

这对我正在做的事情来说已经足够了,但它有点像创可贴的解决方案。

根据数学明确性的要求:

在天文学中,震级规模大致如下:

mu = -2.5log(flux) + mzp # apparent magnitude

其中 mzp 是每秒看到 1 个光子的幅度。因此,更大的通量等同于更小(或更负)的视震级。我正在为使用多个组件功能的源制作模型。前任。两个具有不同 sersic 索引的 sersic 函数,P1在内部组件上具有外部1-P1截断,在外部组件上具有内部截断。这样,当向每个分量添加截断函数时,由半径定义的幅度将变得非常大,因为 mu1-2.5*log( P1) 在P1渐近接近零时变得非常小。

TLDR:我想知道的是,是否有一种方法可以保留精度不足以与零区分开来的浮点(特别是在渐近接近零的函数的结果中)。这很重要,因为当取这些数字的对数时,结果是域错误。

非对数 P1 中的输出开始读取零之前的最后一个数字是5.551115123125783e-17,这是一个常见的浮点算术舍入误差结果,其中所需的结果应该为零。

任何投入将不胜感激。

@user:Dan 没有放我的整个脚本:

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]

还要注意零之前的浮点数。

4

2 回答 2

5

回想一下,双曲正切可以表示为 (1-e^{-2x})/(1+e^{-2x})。通过一些代数,我们可以得到 0.5*tanh(x)-0.5(你的函数的负数)等于 e^{-2x}/(1+e^{-2x})。这个的对数是-2*x-log(1+exp(-2*x)),它在任何地方都有效并且稳定。

也就是说,我建议您更换:

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

使用这种更简单、更稳定的方法:

r = np.arange(0.1,100.01,0.01)
#r and xvals are numpy arrays, so numpy functions can be applied in one step
xvals=(2.0 - B(rb1,drs1) ) * r / rb1 + B(rb1,drs1)
P1=2*xvals+np.log1p(np.exp(-2*xvals))
于 2014-01-18T00:06:38.237 回答
1

你可以尝试两件事。

(1)蛮力方法:找到一个可变精度浮点算术包并使用它而不是内置的固定精度。我正在解决您在 Maxima [1] 中的问题,我发现我必须大大提高浮点精度以避免下溢,但这是可能的。如果你愿意,我可以发布千里马代码。我想有一些适合 Python 的可变精度浮点库。

(2) 用泰勒级数或其他某种近似来近似 log((1/2)(1 + tanh(-x)) 以完全避免 log(tanh(...))。

[1] http://maxima.sourceforge.net

于 2014-01-18T00:50:57.707 回答