水素原子の電子雲ジェネレータ＆ビューア完成！

退屈なおうち時間の合間に楽しんでもらおうと思って、ようやく完成させました！
水素原子様の電子雲、ジェネレータとビューアです！
(何十年越しの夢だったろう…自粛期間はすぎましたが、とりあえずまた計画を立ち上げられてよかった…！)




ExcelファイルDLはこちらから！ 

ｓ

ラゲール備忘録

差分方程式化。

ずっと前だったから忘れてたけど、なるほどルジャンドルよりラゲールの方が簡単な気はした。
先にラゲールに手を出して正解だったようだ。


変数ρ(r)の範囲は0から40まで。刻み幅は0.02だとまだちょっと発散するから
0.01にすれば収まりそう。4000行使うことになるが
そうすれば、入力規則で量子数を整数に縛れる。




・関数直交判定ツール
・ルンゲクッタによる容量圧縮
・入力規則のアピール
・変数関数変形動画

ルジャンドル陪多項式の検算

まずいきなり訂正です。
またしても同じ式を訂正することになってしまいました…。
昨日の訂正日記


sqrt(sumsq()*dz/count())

と書いていましたが、正しくは

sqrt(sumsq()*dz)

で、個数で割って平均を取るというのが間違っておりました。
てっきり実効値とまったく同じものと思いこんでいたらしく、いつからそう思っていたのか謎です


この間違いが見つかったのは今さっき解析解と照らし合わせていたときのことで、正しく計算しますと、このように解析解と数値解はほぼ一致します。
量子数lとmが大きい場合は、おそらく刻み幅dzを細かくとらないと一致しないのではないかと思います

またしても訂正です。

さっきの日記の規格化の式が間違っておりました。



=sqrt(average())
とか
=sqrt(sum()/count()) とかでいいです。
（ぶっちゃけ=sqrt(sumsq()/count())


ではなく


=sqrt(average())
とか
=sqrt(sum()*dz/count()) とかでいいです。
（ぶっちゃけ=sqrt(sumsq()*dz/count())

dzをしっかり掛け算して、積分(あるいは和を取る)してください

ルジャンドル陪多項式の偶奇と規格化

きのうの続きですが、訂正があります。

この図の2番目の関数Pの初期値ですが、if(E7=0,0.01,1)ではなく
正しくはif(E7=0,1,0.01)でした。



さて、z=0からz=1まで入力していきますと
たとえばl=m=1ですと、以下のような感じになります。



これを、負のzにも拡張しましょう。
偶関数の場合は、負の変数zでの関数Pの振る舞いは
P(-z)=P(z)ですし
奇関数の場合は、
P(-z)-P(z)なので、

z=1の下に以下のように符号反転したzを並べて

先ほどと同様にl+mの偶奇に合わせて、関数の偶奇を決めます。
つまり、l+mの偶奇のセルであるE7セルに応じて、zがプラスの0.05の関数Pの値を符号反転するかそのままか分岐させて、順番にコピペしていきます。
(絶対参照などにご注意ください)



そうして、いちおうの関数が出来上がってきました。
たとえば、l=m=1の偶関数ならこうで

l=1、m=0の奇関数ならこのような感じになります




最後に、規格化をします。
これは、ある範囲の中に、粒子のいる確率が総合して１００％になるようにするものです
5月20日の日記を参考にします。

（絶対値の）2乗を取って範囲内で積分したものが１になればよいので、

このような列を設け、みなさん大好きなsum関数で和をとり
データの個数で割って平均を取ってから、その平方根を取ります。(電圧の実効値RMSと一緒です)
=sqrt(average())
とか
=sqrt(sum()/count())
とかでいいです。
（ぶっちゃけ=sqrt(sumsq()/count())もありですが）

この値が、規格化定数となるため、算出した関数をこの定数で割ると、規格化されます。

ルジャンドル陪多項式をExcelに実装してみよう！

おとといの続きです。
ルジャンドル陪多項式を数値的に算出してます。


dz=0.05と定め、とりあえず量子数lとmをl=m=0としておきます。


z=0から始めて、z=1までdzだけ増えるようにしてみます。



それから、関数の偶奇を決める項目を設けます。
このルジャンドル陪多項式は偶関数か奇関数しかありません。
偶関数だったら0、奇関数だったら1となるように、量子数mの下に0か1の数値を入れます



関数Pの偶奇が反映されるように初期値を2つ決めます。
位置に関する2階微分方程式なので初期条件では変なのですが、境界条件は計算しづらいので
偶関数か奇関数しかないことを利用して、初期条件のように計算しています。

実は、量子数lとmの和の偶奇が関数の偶奇に関係してくるので
=mod(l+m,2)と入力してしまいます
2で割った余りを意味しています。


さらに、1-z^2が頻出するので、関数Pの右隣に1-z^2の列も作っておくと便利です

そして、3つ目の関数Pから、計算開始です！

この式を入れますが、

3番目のPの式に代入する1-z^2やzに3番目ではなく2番目のzの値を参照しています。

僕も今これに気づいたのですが、3番目のPに3番目のzを参照するより、2番目のzを参照するほうがグラフの見た目が綺麗だということに気が付きました。
そして、こうすることで、z=±1のときの1-z^2、つまり0除算を直前で回避することができるようです。


これをz=1まで続けると、ルジャンドル陪多項式の半分ができあがりです。

相対・絶対・複合参照にお気を付け下さい

ルジャンドル陪多項式を数値計算で算出してみよう！

水素原子様の波動関数の一部である、球面調和関数の
緯度方向に関する微分方程式


を解くと、ルジャンドルの陪多項式が得られます。


これを、数値的に解いてみることにします。


まず、上述の微分方程式を差分方程式に変えます。
zの関数P(z)のzによる1階微分は

このように色々定義できますが、2階微分との絡みがあるため、今回は3つ目を採用しようと思います

また、2階微分のほうは

このように差分化されますね。


これを踏まえて微分方程式を差分方程式として解いていってみましょう

元の式はこのようになり

今回はP1とP0を初期値として与えてP2を求めるようにしたいので、この式をP2の式に変形します。





今日はおねむなんでこの辺で。
発散項どうにかならんかなぁ
そもそもの微分方程式の両辺に1-z^2をかけちゃったらどうなるやろか

球面調和関数の世界地図

3日くらい坊主にして申し訳ない！

球面調和関数の世界地図の進捗はこれです！
l=3、m=2の例
上下が緯度、中央が赤道で、左右が経度

まだルジャンドる前の緯度方向微分方程式を差分に

変数θとθの関数Θが紛らわしいので、ここからはFという名前にします


内側から順に、微分を差分にしていきます。


差分化はここで完了です。

あとは、F2の式に整理するだけです。


以上です

ルジャンドル陪関数の規格化

ルジャンドル陪関数P(z)を微分方程式から数値計算で算出したんですよ。

wikiに書いてある方法だと、一旦zの関数P(z)として出すんですけど
(zは電子の回転軸z軸の変位みたいなもんです)

規格化するときに解析的な書き方だと、z=cosθとしてからP(cosθ)の2乗を取って
さらにsinθを掛け算して積分するんです。



どうにかしてdzかdθかどっちかに統一したくてですね

じゃあ規格化の積分のときにdθをdzに置換すればいいじゃないかってなりまして

z=cosθなので、dz/dθ=-sinθですし、
θ=0の場合はz=1、
θ=πの場合はz=-1なので


この当たり前の式たったこれだけになるんですよ。
なぜか掛け算してたsinθがきれいさっぱり消えるんです。



＝＝＝＝＝＝＝
ただ、こういうアプローチもいいんですが、
ルジャンドル陪関数を導出する微分方程式

が、z=±1で発散してしまうので、できればジワジワ近づけたいんですよね。
でも微小量を変化させるのは面倒くさいなあと。

だったらdθ使えばいいんじゃね？
ってなりましてね、
さすが人類の叡智の結晶と言いましょうか
原子の緯度の基準を赤道ではなく北極に取った神采配と言いましょうか
dθは固定のまま、θ=0,π付近ではじわじわと発散ポイントに向かってくれるわけですよ


じゃあその、ルジャンドル陪関数を求める微分方程式の
変数がzじゃなくてθのやつはどこにあんの？


探しだしましたよ昨日の寝る前に！(スクショしといてよかった～)
球面調和関数のページの「Y(θ,φ)の意義」の項目のですね
折りたたまれてる「証明」のところを開くと


(kは方位量子数lに読みかえ可能です)
っていう元々の微分方程式があるんですよ！

数値的にはこちらを解いても何の問題もないため、これを直接解くことで
発散ポイントであるz=±1(原子の北極と南極)を精度よく算出することができるはずなんです

そのあと、規格化の際にもθで規格化が可能なのです

球面調和関数の2重積分

あんまりブランクが長いんで、勘違いしてたことがありましてね

たとえばここのdzx軌道の波動関数で、球面調和関数Θ(θ)Φ(φ)の極座標表現が

Θ(θ)Φ(φ)∝sinθcosθcosφ

って傾向があって、規格化定数を求めるとしますよ


波動関数の2乗を取って、sinθを掛け算して2重積分するじゃないですか。

僕は何を血迷ったのか、θに関する積分とφに関する積分が積になっているので、部分積分が必要なんじゃないかとかわけのわからない勘違いをし始めましてね
そもそも2重積分の2要素って、足し算じゃなく掛け算でしか結合できないのに、何を考えているんだ自分？ってなったんですよ。


しかも、わざわざ変数分離してんだから、これは2重積分なんてもんじゃなくて
θに関する積分とφに関するただの積分の積に分けることが可能じゃないですか。
なにやってだアホーって感じですよ




今はちょっと、ルジャンドル陪関数が数値的に解けてきたので、これを規格化する準備をしています。
解析的には変数をtという回転軸方向の変位に一旦直して、半分くらい円筒座標系みたいな感じにしてから、tをcosθに置き換えて、最終的にcosθの式にしてからcosθで積分して規格化するんですけど
数値計算だったら、そもそも微分方程式を解く際に変数はtである必要なくね？
って思ったので、その辺も考慮に入れる予定です。変数をcosθにしたまま、dθで微分した微分方程式でも解けるんじゃないかと企て中です

ルジャンドル陪関数の偶奇と量子数の偶奇

ルジャンドル陪関数P(x)、変数のxを水素原子では極座標の緯度方向θに置換して
P(cosθ)とかにしやがるから

xは奇関数だけど、cosθは偶関数なので、一両日くらい混乱してた。

解く際はP(x)とかやって、xとか√(1-x^2)とかが大量に出てきて
xがcosθだと√(1-x^2)はsinθになるんだけど

式の次元だけ見ると
√(1-x^2)は2乗してまたその1/2乗するからこれは奇関数やな
とか思ったら大間違いで、グラフをよく確認するとこれ偶関数なんよね

だから一度勘違いした結果、表を参照して磁気量子数mと方位量子数lを用いて
2m+lの偶奇がルジャンドル関数の偶奇を決めるとか思ってしまったけど


ここでなんか知らんけどwiki見返して
2m+lなんて項目あったっけ？とかわれながらわけわかんない着目をしたおかげで
√(1-x^2)が偶関数であることに気づくことができた。

そうすると、cosθとかで表してたときにはmの偶奇だけにP(cosθ)の偶奇が依存してたのが
P(x)の場合は2m+lなどではなく、l+mの偶奇に関数P(x)の偶奇が依存しているとわかった。


数値計算する際にこの偶奇情報はかなり大事だったりする。



あとは、これをこのP(x)にφの関数をかけてから2乗してさらにsinθをかけてからdθで積分してからさらにdφでも積分して
規格化定数を得るわけだけど
微分方程式を変数xで解いておきながら、数値積分するときはθでやれとかモヤっとするよな



なんかもう問題が次々割り込んできてなかなか到達できる気がしない(笑)



ところで、ルジャンドル陪関数を微分方程式から数値的に導出するのはほぼバグが取れて
オイラー法でも精度上げればなんとかなることがわかったのでよかった。

何が主なバグだったんだったかな、
ただもう1つの要因として、-m^2って項が微分方程式にあって
Excelだとこれが曲者で、変な計算結果を吐き出しかねないので、きちんと(-1)*m^2って記述したのも大きいのかもしれない


あ、思い出した。この-m^2の分母に1-x^2があったから、x=±1で発散するんよね
あろうことかこの発散ポイントから解くことを始めてしまったからそりゃ終始発散もするわ
って納得したんだった

久しぶりに差分方程式に手を出した

発散した。
非対称に。
途方に暮れた。




僕としては理屈を理解してるオイラー法でちゃっちゃと確かめたいのに
オイラー法はどちらかというと精度が悪いらしく
久々なのもあってか、何が原因でうまくいかないのかわからず、とりあえず今日一日途方に暮れた。


特に、固有値のある微分方程式なんかは、発散するのがデフォみたいなところがあるし
オイラー法の精度だと、解析計算での固有値と数値計算での固有値が普通にズレるため
試行錯誤が必要だし
球面調和関数(今日はまだルジャンドル多項式)の場合、量子数が2つあるので
どうやって量子数を探るのかがまずわからない


こうなると、もう最初からルンゲクッタ法やっとけよって気もしてくる


そもそも久しぶりなので、ちゃんと式を記述できているかどうかが結構曖昧だったりする。
だからこそ、オイラー法でもできたよ！っていう経験がほしいっていうのもある。


ルジャンドル陪関数は…たしかまだやってない…と思うんだけど
どうだったかな。
ノウハウとか知らんし、やってないはずだ。たぶん


気になる点とりあえず1つ
二階微分の差分表現(y2-2y1+y0)/dt^2

に対して、1階微分の差分表現をどうするか

(y2-y1)/dt
にするか

(y1-y0)/dt
にするか

間を取って
(y2-y0)/dt/2
にするか



そういえば極座標の緯度φじゃなく経度θに関する微分方程式はどう解いたんだったっけ
そもそも数値的に解く必要性を今は感じないんだけど
昔どうやってたか、謎だ。

特に気になるのは境界条件
ぐるっと一回りして元に戻るはずなんだから、数値計算したときに発散するわけないんだよな
別に井戸型ポテンシャルと同じように解いてよかったんだっけか？
発散しなければ固有値がおｋみたいな？
なんだろう？まったく記憶にない



でもブログにはアップされてるんだよなぁ。

球面調和関数の規格化例

この計算は楽しかった。

例えるならそう、フルパワーグリッドマンのような変形合体。
複素数はロマンだよねぇ。
各々の三角関数をアシストウェポンとすると、それらが一旦複素の世界に変形して
掛け算されることでグリッドマンに合体し、
三角関数に変形されなおす。




数式エディタで直に式展開するスタイル、久々にやったけどなまってなくてよかった。
こっちのほうが紙媒体より早く正確にできる場合は確かにある。

それでも紙媒体に数式を書くのは
閃きと、計算の自由度を得たいためかな。


以前僕の先生が、計算してると呼吸を忘れるって言ってた気がするんだけど
今更冗談交じりにググってみたら「電子メール無呼吸症候群」ってのがガチであるみたいで驚いた。

僕が学生のころはネットの黎明期で、まだポケベルも使われてたころだから
先生が電子メール無呼吸症候群を知っていたかどうかは定かではない。


でも久々に数式エディタで長時間作業してから一息ついたときに
やたらあくびがでることに気づいた。

やっぱりこの計算の趣味はやるだけで体のあちこちをむしばむ可能性があるようだ。
実際、腰痛や逆流性食道炎にはなったし、いつも左足がしびれる。

VDTも苦手なんだよな。

スマホばっかりやったらそれはそれで腰が痛くなるし


寿命が縮む云々より、その寿命まで計算力をいかに劣化させないかが問題だ





＝＝＝＝＝＝＝＝＝
追記
ああそうそう、球面調和関数の規格化の際、緯度φと経度θで積分するけど、波動関数としての球面調和関数の絶対値の2乗に、なぜかsinθを掛け算してから積分するらしい。

これもだけど、そもそも極座標のラプラシアンとか、最近原理をあまり理解しないで使っていることが多くなってきた。
今後使うことになるだろう、オイラー法に変わるルンゲクッタ法なんかもあまり理解してない
もちろんラゲールやルジャンドルの多項式も、最初からノータッチだ。触れる気がしない


球の体積を求める際にも、なんか微小量としてsinθを掛け算していたような…？

この辺はもうどうしようもない老いの部分として、開き直った方がいい気がする


ただ、極座標のラプラシアンとかについては、もしかしたらテンソルの運用を身につける際に簡単に思えてくるのかもしれない

ガチで自分の一生の間にどこまでできるんだろう…昔は冗談半分で笑っていたけど、だんだん笑えなくなってきた。
まあ僕は一線を走るような人ではないので、どんな途中でもいいんだけどね
最後まで楽しくはありたい

俺たち、取説をどこで読み間違えちまったんだろうな…

極座標のラプラシアン、ズルして導出せずにどっかのサイト参照したんですが
そのサイト、ちゃんと正解を書いているのに読み間違えて


こんな風に書いてしまいました。

正しくはこっちです。
 


緯度方向の波動関数Θ(θ)をたまたま計算していたら、wikiと合わなかったので
おっかしいな～？って思って調べたんでした。
危なかった
具体的にはΘ(θ)=3cos^2θ-1のところで違和感に気づきました。




ところで、昨日はですね
緯度方向の波動関数Θルジャンドルの陪関数のことですが、
一周すると元に戻るはずなので赤道を原点にすると、
偶関数か奇関数になるはずだよな
っていう妄想をしていたんですね

もしかして、北極(か南極)を原点にしても偶関数か奇関数になったりしますかね？
水素原子のwikiではθの取り方は北極を原点にしているので、数式を眺める限りだと
純粋な偶関数か奇関数っぽく見えますね