2016-08-08 1 views
0

私はジュリアの波動関数の確率密度を計算してプロットしたいと思います。私は、次の機能を評価するためのジュリアコードの小さなスニペットを書いた:ジュリアのexp関数が0に評価する

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

ジュリア(不完全)のコードは次のとおりです。

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

私はt0.000000000000...を与える評価し、最後の行。私はBigFloatの精度を設定しようとしました。いいえ、運がない!私は間違って何をしていますか?ヘルプは高く評価しました。

+0

数式の平方根は少し曖昧です。分数の分子と分母の両方を調べるように見えます。しかし、あなたのコードでは、分子の平方根しか取っていません...これは問題ですか? –

+0

通常、 'm = BigFloat(9.10938356e-31)'のようなものを使用したくないことに注意してください。 '9.10938356e-31'はFloat64であり、丸めてからBigFloatに変わりますが、すでにFloat64丸め誤差が発生しています。代わりに、 'parse(BigFloat、" 9.10938356e-31 ")'を使うのが普通です。私はここであなたのエラーだとは思わない(あなたが小数点を入れているため)が、近い将来にいくつかのエラーを停止することがあります。 –

+0

@ColinTBowersあなたは平方根について正しいと分かります。編集されました。助けてくれてありがとう: –

答えて

6

コメントの中で、Chris Rackauckasは、あなたが数式を間違って入力したと指摘しています。私は、我々は我々が調達されているかを見ることができ、それを打破することができます

とにかく質問に答えるために十分面白かった考え出し:そうおおよそので、非常に大まか z=-2e19)

z =-2.0237336022083455711032042949257e+19 ので

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

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

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)はごく少数です。 それは文脈の中で私たちを望みます。非常に小さい番号。


だから、ジュリアはあなたがMPFR onlineを試すことができます ビッグフロートのためMPFRに依存します。精度8192でexp(-2e10)=0 同じ結果です。

今、私たちが気にする精度ではありません。 しかし、むしろ指数の範囲。

MPFRは、IEEEスタイルの浮動小数点数を使用しています。ここで、precisionは仮数の長さで、指数部があります。 2^exponent * mantissa

したがって、指数の範囲には制限があります。

参照:MPFR docs:

機能:mpfr_exp_tのmpfr_get_emin(無効) 機能:mpfr_exp_t mpfr_get_emax(無効)浮動小数点変数に許さ

リターン(現在の)最小と最大の指数。 浮動小数点変数の最小正の値は、半分の2が最小の指数に上げられ、最大値は(1 - ε)倍2が最大指数まで上げられた形式を持ちます。ここで、εは考慮される変数。

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)

だから我々はできるしたがって
log(2,e)*z = -29196304319863382016
Base.MPFR.get_emin() = -4611686018427387903

指数(粗い-2.9e19)は指数(おおよそ-4.3e17)許容最小値未満であるからです。 アンダーフローが発生します。

このように、ゼロになる理由についてのあなたの答え。

Int128指数でMPFRを再構築することは可能かもしれませんが、できません。

おそらく、ジュリアはアンダーフロー例外をスローする必要があります。 無料でJulia Bug Trackerの問題として報告することをお勧めします。

+1

あなたはとてもエレガントに説明しました。非常に親切にありがとう!私はgitの問題をオープンしました - [非常に小さな数値のためにexpを評価するとアンダーフローが発生する](https://github.com/JuliaLang/julia/issues/17893) –