tsujimotterのノートブック

日曜数学者 tsujimotter の「趣味で数学」実践ノート

日曜化学(2):3次元空間における電子雲の計算(Python/matplotlib)

2日前に公開した量子力学に関する記事なのですが、たくさんの方に見ていただいて嬉しいです。Twitter上でもたくさんの嬉しいコメントをいただきました。
tsujimotter.hatenablog.com

今日は続きとして、電子雲の可視化 をしたいと思います。


前回の記事では、水素原子におけるシュレディンガー方程式

 \displaystyle \hat{H}\Psi = E \Psi \tag{1}

を考えました。 \Psi は波動関数で  E はエネルギー。 \hat{H} はハミルトニアンという演算子で、定義は次の通りです:

 \displaystyle \newcommand{\d}{\partial}\displaystyle \hat{H} = -\frac{\hbar}{2m_e}\left(\frac{\d^2}{\d x^2} + \frac{\d^2}{\d y^2} + \frac{\d^2}{\d z^2} \right) + V({\bf r}) \tag{2}

この方程式をデカルト座標  (x, y, z) から球面座標系  (r, \theta, \phi) に直して、変数分離によって解を求めるという方法を紹介しました。

変数分離

 \Psi(r, \theta, \phi) = R(r) \, Y(\theta, \phi) \tag{3}

によって、動径波動関数  R(r) と球面調和関数  Y(\theta, \phi) に分けられるわけですが、前回の記事では特に球面調和関数  Y(\theta, \phi) について可視化を行いました。


しかしながら、球面調和関数が教えてくれるのは「どの方向に電子が多く分布しているか」という情報です。これだけだと「3次元の中でどの辺に電子が分布しているのか」という情報はわかりません。この情報を知るためには「動径方向」の分布を調べる必要があり、したがって   R(r) を計算する必要があります。


そこで、今日は動径方向の分布を考慮に入れて、実際に電子雲の計算をしてみようと思います。最終的には、Pythonのmatplotlibを用いて次のような図を描くのを目標としたいと思います。

f:id:tsujimotter:20210704033300j:plain:w300
実際、3次元での電子分布を考えるにあたっては、それなりに工夫が必要です。上記の図では「確率密度関数の積分値が90%になる領域をプロットするテクニック」を使っているのですが、なかなか実装するのが難しいのでそれについても解説したいと思います。

前回同様、Pythonのmatplotlibを使ったプログラムも掲載しますので、よかったら遊んでみてください。


GIF画像も作ってみました! 楽しいですね!

f:id:tsujimotter:20210704131328g:plain:w400

続きを読む

日曜化学:量子力学の基本と球面調和関数の可視化(Python/matplotlib)

最近、とある興味 *1 から量子力学(とりわけ量子化学)の勉強をしています。

水素原子の電子の軌道を計算すると、s軌道とかp軌道とかd軌道とかの計算が載っていて、対応する図が教科書に載っていたりしますよね。

こういうやつです:

f:id:tsujimotter:20210627223229p:plain:w400
Wikipedia「球面調和関数」より引用
Attribution: I, Sarxos


個人的な体験ですが、予備校の頃は先生の影響で「化学」に大ハマりしていました *2

ここから「Emanの物理学」というサイトの影響で「物理」に目覚め、そこからなぜか「数学」に目覚めて現在に至ります。そういった経緯もあって、化学には大変思い入れがあります。


特にこの水素原子の軌道の図は当時から気になっていて、自分で描いてみたいと思っていました。先日ようやく理解でき、実際に自分で描画できるまでになりました。以下がその画像です:

f:id:tsujimotter:20210627222555j:plain:w560

これはタイトルにもある「球面調和関数」と呼ばれる関数を可視化したものになっているのですが、この図を描くのは私の中でとても感動的な体験でした。

そこで今回の記事では上記の図の描き方を、実際に使ったプログラムを合わせて紹介したいと思います。図を描くにあたって私の中で新しく理解したことも多々ありましたので、そのノウハウも紹介したいと思います。

*1:そのうちブログでも書きたいと思っています。

*2:その頃のエピソードについては、以前記事にしたことがありました: tsujimotter.hatenablog.com

続きを読む

位数6の群の分類

私がスタッフとして携わっている日曜数学会というイベントが、今月の20日に 6周年 を迎えます。めでたいですね。また、先日私の年齢も  36 歳になりまして、つまり  36 = 6^2 というわけですね。ダブルでめでたい(?)。

せっかくなので、何か  36 にまつわる発表をしたいなと思いました。私は最近「表現者のための数学」という講座をやっていまして、その講座のテーマが群論です。準備のために群論の勉強をしているところです。

そこで、位数  36 の群の分類 をやろうと思い立ったのです。

しかし・・・



f:id:tsujimotter:20210614233324p:plain:w400

そう、挫折したのです。

群の分類というのは思った以上に難しいもので、位数の約数が増えれば増えるほど難しくなります。なんとかなるだろうと思ったのですがダメでした。

私が位数36に挑戦するのは少し早かったようです。


そこで、練習台としてもう少し簡単な位数から始めようと思いました。それが表題の 位数  6 の群の分類 です。

続きを読む

(a○+b)×(a△+b)=(a□+b)になるa,bの条件と中国剰余定理

数学ファンの鯵坂もっちょさんがツイートしていた問題が面白かったので、今日はその問題について考えてみたいと思います。


もっちょさんの問題は、任意の  a\text{○}+b, \; a\text{△}+b という形の数の積

 (a\text{○}+b) \times (a\text{△}+b)

が再び

 a\text{□}+b

と表せるような整数  a, b の条件は?という問題ですね。


この問題、見た目以上に広がりがある問題だと思います。「中国剰余定理」や「天に向かって続く数」にもつながったりします。お楽しみに!

続きを読む

フェルマー商

今日は整数論の面白概念の一つである フェルマー商 を紹介したいと思います。


まず、フェルマーの小定理という、合同式を考える上で大変有用な定理から話を始めます。

定理(フェルマーの小定理)
 a p で割り切れない任意の整数とする。このとき
 a^{p-1} \equiv 1 \pmod{p} \tag{1}

が成り立つ。


つまり、この定理は  p と互いに素な  a について、 a^{p-1} - 1 が必ず  p で割り切れると言っているわけですね。

 \bmod{p} の剰余も十分面白いのですが、今回はさらに話を進めて  \bmod{p^2} の合同式を考えてみたいと思います。

続きを読む

自由研究:レピュニットが素数で何回割れるのか問題

一つ前の記事で「LTEの補題(指数持ち上げ補題)」という有名な補題について勉強しました。せっかくなので、この補題を何か他のネタにも使えないかと考えました。

LTEの補題とは、こんな主張の補題でした:

LTEの補題(指数持ち上げ補題)
 p を奇数の素数とする。 p で割り切れない相異なる整数  x, y x\equiv y \pmod{p} なる関係を満たすとき、任意の正の整数  n に対して次が成り立つ:
 \operatorname{ord}_p(x^n - y^n) = \operatorname{ord}_p(x - y) + \operatorname{ord}_p(n)

これ自体の証明は、たとえば INTEGERSの記事 をご覧になってください。


色々考えていくうちに、レピュニットに応用できそうだ というところに思い至ったので、そのことについてまとめたのがこちらの記事です。

今回紹介するレピュニットの性質については、なかなか面白いと思います! なにせ、私はこの性質が載っているサイトを見たことがありません。

よかったらぜひ楽しんでいってください!

2021.05.26 公開後の追記:
公開後に調べていたら、どうも定理3の話は東大の入試問題になっていたようです。東大すごいですね。(記事末尾にリンクを記載しています。)

続きを読む

線形合同法(擬似乱数生成法)の周期

世の中の現象の中には「ランダム」な現象が多々あります。たとえば、サイコロを振るのは分かりやすいランダムな現象の例です。他にも天気や地震、ギャンブルなども分かりやすいランダムな現象の例です。

一方で、コンピュータの中で行われる計算は、一定のアルゴリズムにより定められた計算が順次実行されることになるため、原理的にはランダムになることはありません。

このことはコンピュータ内でランダムな現象をシミュレーションするにあたって問題になります。そこで「決められた手順によって生成されるランダムっぽく振る舞う数列」をいかにして作るかが重要になるわけですね。このような数列を擬似乱数と言います。

コンピュータによって擬似乱数を生成する手法は、これまでもさまざまな手法が提案されています。今回は、その中でも特に有名な手法である 線形合同法 について考えたいと思います。

最近はもっと良い疑似乱数の生成法が発見されているので、線形合同法はあまり使われていないと聞きます。しかし、一昔前のC言語等の乱数生成には、この方法が普通に使われていました。


簡単に線形合同法について説明しましょう。線形合同法とは次のようなアルゴリズムです。

整数の3つ組  (a, b, M) を考えます。とくに、 a > 0, \; b \geq 0, \; M > 0 とし、 a < M, \; b < M としておきます。

 x_0 0 \leq x_0 \leq M-1 であるような整数とし、整数  x_1, x_2, x_3, \ldots を次の漸化式によって定めます:

 x_{n+1} \equiv a x_n + b \pmod{M} \tag{1}

ただし、 x_n \{0, 1, 2, \ldots, M-1\} の中からとることにします。

このように数列  \{x_n\} を選ぶと、 x_0 を決める毎に異なる疑似乱数の列が得られるというわけです。


たとえば、 a = 13, b = 5, M = 24 とします。このとき、 x_0 = 8 とすると

 \begin{align} x_0 &= 8, \\
x_1 &= 13\times 8 + 5 = 109 \equiv 13 \pmod{24}, \\
x_2 &= 13\times 13 + 5 = 174 \equiv 6 \pmod{24}, \\
x_3 &= 13\times 6 + 5 = 83 \equiv 11 \pmod{24}, \\
& \vdots \end{align}

のように数列  8, 13, 6, 11, \ldots が得られます。

今回の数列は結構優秀で、数列をこのまま続けていくと、こうなるのですが・・・

 {\bf 8},13,6,11,4,9,2,7,0,5,22,3,20,1,18,23,16,21,14,19,12,17,10,15,{\bf 8}\ldots


実は、 8 からスタートして  8 に戻ってくるまで、 0 から  23 までの数がすべて1回ずつ出力されています。つまり、この数列の周期は  24 というわけですね。周期は法  M を超えることはないので、これが最大周期になります。


パラメータ  (a, b, M) によっては最大周期とならないケースも存在します。たとえば、 a = 15, b = 5, M = 24 とすると

 8,21,11,15,23,15,23,15,23, \ldots

となりますが、これは  15, 23 が繰り返されるので、周期  2 となります。


あとで述べるように、一般に線形合同法によって得られる数列は周期  \lambda を持ち、その周期は  \lambda \leq M となります。

今回のメインテーマは、パラメータ  (a, b, M) がどのような条件のとき周期最大( \lambda = M)となるのか、という問題です。周期は長ければ長いほど疑似乱数としては優秀だと考えられるので、気になる問題設定です。


線形合同法において最大周期になる必要十分条件は知られていて、以下の定理1が成り立ちます。

定理1
パラメータ  (a, b, M) に対する周期を  \lambda とする。このとき、次の必要十分条件が成り立つ:

 \newcommand{\HLNO}{\unicode[\sans-serif,STIXGeneral]{x306E}}  \lambda = M \;\; \Longleftrightarrow \;\; \begin{cases} \text{(i)}\;\;\; (b, M) = 1 \\
\text{(ii)} \;\;\; \text{任意}\HLNO{奇素数} \; p \; \text{に対して} \;\; p\mid M \; \Longrightarrow \; p \mid (a-1) \\
\text{(iii)} \;\;\; 4\mid M \; \Longrightarrow \; 4 \mid (a-1) \\
\end{cases}


実際、 a = 13, \; b = 5, \; M = 24 のケースでは、条件 (i), (ii), (iii) を満たすので最大周期  \lambda = 24 になったというわけですね。


しかし、一体どうしてこのような必要十分条件になるのか気になりますね。そこで、今回の記事は 定理1を証明 してみたいと思います。


自力で考えてみてもなかなか難しかったので、証明がどこかに載っていないか調べてみましたが、ネット上で一般的なケースについての証明は見かけませんでした。

その後、Knuth先生の "The Art of Computer Programming Volume2 Seminumerical Algorithms" に証明が載っていることを知り、今回参考にさせていただきました。


しかしながら、証明を読んでいて納得いかない部分がいくつかありました。

こちらの記事では、大枠の証明は上記に従いつつも、いくらか自分流に修正したものをまとめたいと思います。

続きを読む