📐

Juliaは誤差伝播の法則を扱えるか？

2022/12/04に公開

Julia

数値計算 Advent Calendar 2022 4日目です. 誤差の話が続いているので便乗しました.

誤差伝播の法則は学部の必修の実験でサラッと習った記憶があります. 理論系の研究室に進みましたが, 計算値と実験値を比較するために必要になることがあり, いつも公式とにらめっこしながら計算しています. 誤差部分でいちいち√の計算をしないといけないのが大変ですね. そこで, 良さそうなパッケージを探していたところ Measurements.jl を見つけたので, ちゃんと計算できているか確かめつつ, 使ってみることにしました.

このノートでは, Measurements.jlを用いて誤差付きの計算を行います. まずは基本的な使い方を確認し, 次に誤差伝播の法則が成り立つことを確認します. 最終的に水素原子のエネルギー準位の誤差を評価することを目的とします.

環境・パッケージ

初回のみ, REPL上, またはJupyter Notebook上で, Measurements.jlをインストールしましょう. 下記のコメントアウトを外せばNotebook上でインストールできます. 使用する前には毎回using Measurementsを宣言します. また, versioninfo()でJuliaのバージョンや端末の情報を確認できます.

入力

# using Pkg
# Pkg.add("Measurements")
using Measurements
versioninfo()

出力

Julia Version 1.8.3
Commit 0434deb161 (2022-11-14 20:14 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: 8 × 11th Gen Intel(R) Core(TM) i7-1185G7 @ 3.00GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-13.0.1 (ORCJIT, tigerlake)
  Threads: 1 on 8 virtual cores

使い方

公式ドキュメントを参考に試していきます.

宣言

下記の3つはいずれも同じ意味のようです. さりげなくやっていますが, 地味にすごいです. CODATAの値も全て12.34(56)の形で提供されているので, これをパースできるのはめちゃくちゃありがたいことです.

julia> measurement(12.34, 0.56)
12.34 ± 0.56

julia> measurement("12.34(56)")
12.34 ± 0.56

julia> 12.34 ± 0.56
12.34 ± 0.56

演算・関数

基本的な演算はもちろん, 関数の利用もできるようです. 特殊関数にも対応しているらしいので, 興味がある人は公式ドキュメントを読んでください.

julia> a = 12.3 ± 0.4

julia> b = 23.4 ± 0.5

julia> a + b
35.7 ± 0.64

julia> a - b
-11.1 ± 0.64

julia> a * b
288.0 ± 11.0

julia> a / b
0.526 ± 0.02

julia> a % b
12.3 ± 0.4

julia> a ^ b
3.2e25 ± 4.7e25

julia> 2a
24.6 ± 0.8

julia> a ^ 2
151.3 ± 9.8

julia> 2 ^ a
5000.0 ± 1400.0

julia> sqrt(a)
3.507 ± 0.057

julia> sin(a)
-0.26 ± 0.39

julia> cos(a)
0.96 ± 0.11

julia> tan(a)
-0.27 ± 0.43

julia> log(a)
2.51 ± 0.033

julia> exp(a)
220000.0 ± 88000.0

値と誤差の取り出し

この時のデータは以下のように格納されています.

julia> a = 12.3 ± 0.4

julia> dump(a)
Measurement{Float64}
  val: Float64 12.3
  err: Float64 0.4
  tag: UInt64 0x0000000000000006
  der: Measurements.Derivatives{Float64}
    parent: Measurements.Derivatives{Float64}
      parent: #undef
      key: Tuple{Float64, Float64, UInt64}
        1: Float64 2.61e-321
        2: Float64 9.06030910405e-312
        3: UInt64 0x0000000000000098
      value: Float64 3.16e-322
    key: Tuple{Float64, Float64, UInt64}
      1: Float64 12.3
      2: Float64 0.4
      3: UInt64 0x0000000000000006
    value: Float64 1.0
nothing

従って, 値と誤差は以下のように.valと.errで取り出せます. 変数同士の和についても(a + b).valのように取り出せます.（同じことがMeasurements.value()とMeasurements.uncertainty()でもできます.）

julia> a = 12.3 ± 0.4

julia> a.val
12.3

julia> a.err
0.4

julia> (a + b).val
35.7

julia> (a + b).err
0.6403124237432849

誤差伝播の法則

かなり使えそうなパッケージですが, ちゃんと伝播則に基づいた計算を行っているのか確かめていきましょう. 伝播則はこちらやこちらに載っています. 掛け算の表記が少し違いますが, 式を変形すれば一致することがわかると思います.

\begin{align} (M_1 \pm u_1) + \left(M_2 \pm u_2\right) &= \left(M_1 + M_2\right) \pm \sqrt{u_1^2 + u_2^2} \\ (M_1 \pm u_1) - \left(M_2 \pm u_2\right) &= \left(M_1 - M_2\right) \pm \sqrt{u_1^2 + u_2^2} \\ (M_1 \pm u_1) (M_2 \pm u_2) &= M_1 \cdot M_2 \pm M_1 \cdot M_2 \sqrt{\left(\frac{u_1}{M_1}\right)^2 + \left(\frac{u_2}{M_2}\right)^2} \\ \frac{(M_1 \pm u_1)}{(M_2 \pm u_2)} &= \frac{M_1}{M_2} \pm \frac{M_1}{M_2} \sqrt{\left(\frac{u_1}{M_1}\right)^2 + \left(\frac{u_2}{M_2}\right)^2} \\ (M_1 \pm u_1)^{(M_2 \pm u_2)} &= M_1^{M_2} \pm \sqrt{\left(M_2\cdot M_1^{M_2-1}\cdot u_1\right)^2 + \left(\log{M_1}\cdot M_1^{M_2}\cdot u_2\right)^2} \end{align}

Measurements.jlの計算がこれらの伝播則を満たすか確かめましょう.

四則演算と累乗

まずは和における伝播則について検証する例を見せます.

julia> a = 12.3 ± 0.4
julia> b = 23.4 ± 0.5

julia> a + b                        # Measurements.jlによる計算
35.7 ± 0.64

julia> a.val + b.val                # 公式をそのまま実装して検算
35.7

julia> sqrt(a.err ^ 2 + b.err ^ 2)  # 公式をそのまま実装して検算
0.6403124237432849

ちゃんと誤差の部分も合っているようですね. まとめて他の演算でも合っているか検証しましょう.

julia> a = 12.3 ± 0.4
julia> b = 23.4 ± 0.5

julia> (a + b).val - (a.val + b.val)
0.0

julia> (a - b).val - (a.val - b.val)
0.0

julia> (a * b).val - a.val * b.val
0.0

julia> (a / b).val - a.val / b.val
0.0

julia> (a ^ b).val - a.val ^ b.val
0.0

julia> (a + b).err - sqrt(a.err ^ 2 + b.err ^ 2)
0.0

julia> (a - b).err - sqrt(a.err ^ 2 + b.err ^ 2)
0.0

julia> (a * b).err - a.val * b.val * sqrt((a.err / a.val) ^ 2 + (b.err / b.val) ^ 2)
0.0

julia> (a / b).err - (a.val / b.val) * sqrt((a.err / a.val) ^ 2 + (b.err / b.val) ^ 2)
0.0

julia> (a ^ b).err - sqrt((b.val * a.val ^ (b.val - 1) * a.err) ^ 2 + (log(a.val) * a.val ^ b.val * b.err) ^ 2)
0.0

全て問題ありませんでした. その他にも, 乱数を使って数万パターンのテストを行いましたが, 明らかに伝播則に反するような結果を与えた例はありませんでした.

BigFloatによる演算

誤差の伝播則は満足してるようですが, 数値誤差が原因となって, Measurements.jlによる計算と, 公式を愚直に実装したものの差は0になるとは限りません. しかし, Juliaには任意精度演算（BigFloat, BigInt）が標準搭載されているため, 信頼できる結果を与えてくれます. 例えば, 下記のように, そこそこ大きい桁の数値誤差が出る計算であっても, BigFloatを使えば数値誤差を限りなく小さくできます.

julia> a = 100.0 ± 0.1

julia> b = 1.0e-5 ± 0.1

julia> (a / b).err - (a.val / b.val) * sqrt((a.err / a.val) ^ 2 + (b.err / b.val) ^ 2)
-1.52587890625e-5

julia> a = big"100.0" ± big"0.1"

julia> b = big"0.00001" ± big"0.1"

julia> (a / b).err - (a.val / b.val) * sqrt((a.err / a.val) ^ 2 + (b.err / b.val) ^ 2)
2.373891936439949686868310567390489288555241458404986241569944997707252057661203e-66

関数

四則演算と累乗は既に検証していますので, これらの組み合わせも正しいと言えるでしょう. 多くの関数の実装は単純なTaylor展開にはなっていないと思いますが, ここではTaylor展開を使って誤差を検証できるか試してみようと思います. Taylor展開を使う場合の注意点として, Int64の範囲では $21!$ つまりfactorial(21)がオーバーフローしてしまう問題が挙げられます. 他の言語だと少し難しいかもしれませんが, Juliaでは標準でBigFloatが使えますので, 簡単に解決できます.

まずは指数関数の場合について計算しましょう. $n=11$ で収束しています. ちなみにfactorial(big(i))をfactorial(i)にすると $n=21$ でオーバーフローします.

入力

using Printf
a = big"1.23" ± big"0.04"

for n in union(0:12, 100)
    b = sum([a^i / factorial(big(i)) for i in 0:n])
    @printf("n = %d\t%5f ± %5f\n", n, b.val, b.err)
end

b = exp(a) # ℯ^a
@printf("exp(a)\t%5f ± %5f\n", b.val, b.err)

出力

    n = 0	1.000000 ± 0.000000
    n = 1	2.230000 ± 0.040000
    n = 2	2.986450 ± 0.089200
    n = 3	3.296594 ± 0.119458
    n = 4	3.391964 ± 0.131864
    n = 5	3.415425 ± 0.135679
    n = 6	3.420234 ± 0.136617
    n = 7	3.421079 ± 0.136809
    n = 8	3.421209 ± 0.136843
    n = 9	3.421227 ± 0.136848
    n = 10	3.421229 ± 0.136849
    n = 11	3.421230 ± 0.136849
    n = 12	3.421230 ± 0.136849
    n = 100	3.421230 ± 0.136849
    exp(a)	3.421230 ± 0.136849

ちゃんと収束していますね. 指数関数は累乗の特殊例ですので, 他の累乗についても四則演算で誤差が求められそうな気がします. さて, 次にsin(x)でうまく収束するか試しましょう. ちなみにfactorial(big(2*i+1))をfactorial(2*i+1)にすると $n=10$ でオーバーフローします.

入力

# using Printf
a = big"1.23" ± big"0.04"

for n in union(0:5, 10, 100)
    b = sum([(-1)^i * a^(2*i+1) / factorial(big(2*i+1)) for i in 0:n])
    @printf("n = %d\t%5f ± %5f\n", n, b.val, b.err)
end

b = sin(a)
@printf("sin(a)\t%5f ± %5f\n", b.val, b.err)

出力

n = 0	1.230000 ± 0.040000
n = 1	0.919855 ± 0.009742
n = 2	0.943316 ± 0.013557
n = 3	0.942471 ± 0.013364
n = 4	0.942489 ± 0.013370
n = 5	0.942489 ± 0.013370
n = 10	0.942489 ± 0.013370
n = 100	0.942489 ± 0.013370
sin(a)	0.942489 ± 0.013370

これもちゃんと収束しています. 今回は $x=0$ 周りでのTaylor展開を考えていますので, $x=0$ から離れた点では収束が遅くなります（展開の次数が増えると, 展開点に近いところから厳密なsin関数に近づくアニメーションとか見たことありますよね）. $x=1.23$ から $x=12.3$ のように10倍離れた点でのTaylor展開を計算すると, 確かに遅くなりますが, $n=21$ でちゃんと収束します. a = big"1.23" ± big"0.04"に変えて試してみてください. なお, sinやcosではa = mod(big"12.3" ± big"0.4", 2π)とすれば $x=0$ から離れた値でも収束を早められます.

水素原子のエネルギー準位

これから紹介する相対運動のSchrödinger方程式はあまり見慣れないと思います. なぜなら, 量子化学や固体物理ではBorn-Oppenheimer近似を用いて原子核と電子の運動を分離して扱うことが標準的であり, そもそも考えている座標系からして異なるためです. 量子力学の教科書に載っている水素原子の解説でも, 陽子の質量 $m_\mathrm{p}$ を無限大に近似（ $m_\mathrm{p}\rightarrow\infty$ ）するため, $m_\mathrm{p}$ を扱うことはあまり多くないでしょう.

また, 質点・点電荷のモデルであることに注意しておきます. 水素原子を陽子と電子の2体系と考えます. もちろん, 陽子はクォーク3つからできた多体系ですが, その内部構造や大きさは考慮せず, 陽子と電子を等しく質点として扱うことにも注意しておきます. このことは, 点電荷のポテンシャルを用いるモデルでは暗に仮定されていることです.

Schrödinger方程式

水素様原子の時間に依存しないSchrödinger方程式は次式です.

\left[ -\frac{\hbar^2}{2m_\mathrm{e}}\nabla_\mathrm{e}^2 -\frac{\hbar^2}{2m_\mathrm{p}}\nabla_\mathrm{p}^2 - \frac{Ze^2}{4\pi\varepsilon_0 r} \right] \psi_\mathrm{tot}(\pmb{r}_\mathrm{e},\pmb{r}_\mathrm{p}) = E_\mathrm{tot} \psi_\mathrm{tot}(\pmb{r}_\mathrm{e},\pmb{r}_\mathrm{p})

ただし, 電子と陽子について, それぞれの質量は $m_\mathrm{e}$ と $m_\mathrm{p}$ , 位置は $\pmb{r}_\mathrm{e}$ と $\pmb{r}_\mathrm{p}$ , ラプラシアンは $\nabla_\mathrm{e}^2$ と $\nabla_\mathrm{p}^2$ , 電荷は $-e$ と $Ze$ しています（ $Z=1$ ですが, 一般化できるように $Z$ としておきます）. この方程式から重心運動を分離し, 相対運動のSchrödinger方程式を得ます.

\left[ -\frac{\hbar^2}{2\mu}\nabla^2 - \frac{Ze^2}{4\pi\varepsilon_0 r} \right] \psi(\pmb{r}) = E \psi(\pmb{r})

ただし, 換算質量を $\frac{1}{\mu}:=\frac{1}{m_\mathrm{e}}+\frac{1}{m_\mathrm{p}}$ , 相対運動の座標 $\pmb{r}:=\pmb{r}_\mathrm{e}-\pmb{r}_\mathrm{p}$ としています. 詳細は原田(2007), 高柳(2000), Jensen(2007)などに書いてあります.

厳密解

上記の相対運動の方程式の固有値は次式で与えられることがわかっています（高柳(2000)の式(2.29)を元に一部変更）.

E_n = -\frac{\mu e^4 Z^2}{2n^2(4\pi\varepsilon_0)^2\hbar^2} = -\frac{Z^2}{2n^2} \frac{\mu}{m_\mathrm{e}} E_\mathrm{h}

ただし, $E_\mathrm{h} = \frac{m_\mathrm{e} e^4}{(4\pi\varepsilon_0)^2 \hbar^2}$ はHartreeエネルギーという, 原子単位の一つです. 原子単位系から他の単位系に変換すると実験誤差が混入する様子は後で確かめましょう. また,

\begin{align} \frac{m_\mathrm{e}}{\mu} = \frac{m_\mathrm{e}}{m_\mathrm{e}} + \frac{m_\mathrm{e}}{m_\mathrm{p}} = 1 + \frac{1}{m_\mathrm{p}/m_\mathrm{e}} \end{align}

となるので, 陽子/電子の質量比 $m_\mathrm{p} / m_\mathrm{e}$ の実験値を利用すればよいことがわかります. なお, 大体の量子力学の教科書は $m_\mathrm{p}\rightarrow\infty$ と近似してしまいますし, Born-Oppenheimer近似での電子状態計算の結果を原子単位系で扱う限りでは実験値の入る余地はありません. 原子単位系は原子・分子系のハミルトニアンを無次元化できる単位系だからです.

原子単位・無次元化

具体的な計算をする前に, 変数を無次元化しておきます. 普通は方程式の時点で無次元化しますが, 今回は計算の直前まで無次元化しないでおきました. 原子単位系や, この辺りの事情はIUPAC GreenBook 3.9.2が一番わかりやすいです.

\begin{align} m_\mathrm{e}^* &:= m_\mathrm{e}/m_\mathrm{e} = 1 \\ m_\mathrm{p}^* &:= m_\mathrm{p}/m_\mathrm{e} \\ \mu^* &:= \mu/m_\mathrm{e} = \frac{1}{\frac{m_\mathrm{e}}{m_\mathrm{e}} + \frac{m_\mathrm{e}}{m_\mathrm{p}}} = \frac{1}{\frac{1}{m_\mathrm{e}^*} + \frac{1}{m_\mathrm{p}^*}}\\ E_n^* &:= E_n / E_\mathrm{h} \end{align}

計算ではこれらのアスタリスクがついた無次元量を扱っていきます. ただし, アスタリスクは省略されるため, 勘違いしたり, 間違いを犯す可能性がありますので, 分野外の方にはこのように丁寧に説明する必要があります. 丁寧過ぎて却ってわかりにくくなっているかもしれませんが.

実験値

全て2018年CODATA推奨値を用います. $E_\mathrm{h}$ からeVへの換算係数は

E_\mathrm{h} = 27.211~386~245~988(53)~\mathrm{eV} \\

です. また, 陽子/電子の質量比 $m_\mathrm{p} / m_\mathrm{e}$ には

m_\mathrm{p} = 1836.152~673~43(11)~m_\mathrm{e}

を用います.

具体的な計算

最小固有値（基底状態のエネルギー）は次の値になります.

\begin{align} E_1 = -0.499~727~839~712~381(16)&~E_\mathrm{h} \\ E_1 = -13.598~287~264~287(26)&~\mathrm{eV} \end{align}

単位を $E_\mathrm{h}$ から $\mathrm{eV}$ に変換すると, 有効桁数が1桁減ってしまいましたね. 同様のことは他の単位系でも起こるので, 計算結果を無意味に単位変換しない方が良いです. 学会発表などでは聴衆に合わせて単位変換することになると思いますが, 元の計算結果はどこかに残しておきましょう. LaTeXならコメントアウトして書いておけばよいですし, PowerPointなら発表者用のメモに書いておくとよいと思います.

julia> mₑ = 1
1

julia> mₚ = measurement("1836.15267343(11)")
1836.15267343 ± 1.1e-7

julia> hrev = measurement("27.211386245988(53)")
27.211386245988 ± 5.3e-11

julia> μ = 1 / (1 / mₑ + 1 / mₚ)
0.999455679424763 ± 3.3e-14

julia> E = -μ / 2
-0.499727839712381 ± 1.6e-14

julia> E * hrev
-13.598287264287 ± 2.6e-11

まとめ

Juliaは誤差伝播の法則を扱えました. それもめちゃくちゃすごかったです. Juliaには任意精度演算（BigFloat, BigInt）が標準搭載されているため, 計算の際に生じる数値誤差の影響を除去できるなど, 信頼できる結果を与えてくれます. この記事を書く過程で, 乱数を使って数万パターンのテストを行いましたが, 明らかに伝播則に反するような結果を与えた例はありませんでした. 割り算の誤差(a / b).errは他の計算よりは数値誤差の影響が表れやすいようですが, BigFloatで解決できますので, 不安なら例に示した通りにBigFloatを使えば問題ないと思います. 今後はMeasurements.jlを安心して使います.

参考文献

Measurements.jl

https://juliaphysics.github.io/Measurements.jl/stable/

誤差の伝播則

水素原子

原子単位およびCODATA推奨値

IUPAC GreenBook 3.9.2
Eite Tiesinga, Peter J. Mohr, David B. Newell, and Barry N. Taylor, CODATA recommended values of the fundamental physical constants: 2018, Rev. Mod. Phys. 93, 025010 (2021) https://doi.org/10.1103/RevModPhys.93.025010

追伸

JuliaLang Advent Calendar 2022 15日目でMeasurements.jlが紹介されました！この記事ではブラックボックス的に誤差伝播の公式やパッケージを扱いましたが, こちらの記事では線形誤差伝播の理論についても丁寧に解説されています. また, 他のパッケージとの組み合わせについても解説されており, 非常に参考になります. ご一読ください.