0

I want to calculate and plot the probability density of a wave function in Julia. I wrote a small snippet of Julia code for evaluating the following function:

P = A^2 e^{-(\sqrt(Cm)/\hbar)x^2}

The Julia (incomplete) code is:

set_bigfloat_precision(100)
A = 10
C = 5
m = BigFloat(9.10938356e-31)
ℏ = BigFloat(1.054571800e-34)
t = exp(-(sqrt(C * m) / ℏ))

The last line where I evaluate t gives 0.000000000000.... I tried to set the precision of the BigFloat as well. No luck! What am I doing wrong? Help appreciated.

4

1 回答 1

6

Chris Rackauckas 在评论中指出您输入了错误的公式。我认为无论如何回答这个问题已经足够有趣了

让我们分解一下,这样我们就可以看到我们在筹集什么:

A = 10
C = 5
m = BigFloat(9.10938356e-31)
h = BigFloat(1.054571800e-34)

z = -sqrt(C * m)/h
t = exp(z)

所以 z =-2.0237336022083455711032042949257e+19 非常粗略z=-2e19)

如此粗略t=exp(-2e19) (即t=1/((e^(2*10^19)))这是一个非常小的数字。

考虑到 exp(big"-1e+10") = 9.278...e-4342944820exp(big"-1e+18") = 2.233...e-434294481903251828

是的,朱莉娅说: exp(big"-2e+19) = 0.0000

exp(big"-2e+19)是一个很小的数字。这将我们置于我希望的背景中。非常小的数字。


所以 julia 对于 BigFloats 依赖MPFR你可以在线尝试 MPFR。精度为 8192,exp(-2e10)=0 所以结果相同。

现在,我们关心的不是精度。而是指数的范围。

MPFR 使用有点像 IEEE 风格的浮点数,其中精度是尾数的长度,然后你有一个指数。2^exponent * mantissa

所以指数的范围是有限制的。

请参阅:MPFR 文档:

函数:mpfr_exp_t mpfr_get_emin (void) 函数:mpfr_exp_t mpfr_get_emax (void)

返回浮点变量允许的(当前)最小和最大指数。浮点变量的最小正值是 2 的二分之一乘以最小指数,最大值的形式为 (1 - epsilon) 乘以 2 乘以最大指数,其中 epsilon 取决于所考虑变量的精度.

现在 julia 确实将这些设置为相当默认的 MPFR 编译允许的最大范围。我一直在挖掘 MPFR 源,试图找到它的设置位置,但找不到它。我相信这与 Int64 可以容纳的最大故障有关。

Base.MPFR.get_emin() = -4611686018427387903 =typemin(Int64)>>1 + 1

您可以调整它,但只能向上调整。

所以无论如何

0.5*big"2.0"^(Base.MPFR.get_emin()) = 8.5096913117408361391297879096205e-1388255822130839284

0.5*big"2.0"^(Base.MPFR.get_emin()-1) = 0.00000000000...


现在我们知道了

exp(x) = 2^(log(2,e)*x)

所以我们可以exp(z) = 2^(log(2,e)*z)
log(2,e)*z = -29196304319863382016
Base.MPFR.get_emin() = -4611686018427387903

因此,由于指数(大约 -2.9e19)小于允许的最小指数(大约 -4.3e17)。 发生下溢。

因此,您对为什么得零的答案。

用 Int128 指数重新编译 MPFR 可能(或不可能),但 julia 没有。

也许朱莉娅应该抛出一个下溢异常。Free 鼓励将其报告为 Julia Bug Tracker 上的一个问题。

于 2016-08-08T09:31:28.093 回答