動径と緯度の波動関数をシミュレーションしました(Excelファイル添付)

 wikiによりますと動径R(r)がこうで、緯度Θ(θ)がこうなんだそうで。
 

動径

動径のほうから微分方程式を差分方程式に離散化してみますとこんな感じに。
 
R₊とR₀とR_-は差分化したときのフロントとセンターとバックです。
以降、変数はなるべくwikiのとおりとします。

ρが中心からの距離を規格化・無次元化したもので
dρはまあΔρと考えてください。
nは主量子数、lは方位量子数です。

境界条件の注意点

 
数値計算で解く際の注意点ですが
最初のR₀とR_-がまだ定まってないんですね。
そこで、以前井戸型ポテンシャルでやったとき同様、偶関数か奇関数であることを利用して

投げ落とすか投げ上げるかのイメージでやります。
 
l=0の時はゼロ距離ρ=0でR≠0なので上の図を参照
最初のR≠0、 最初のR₀-R_-≠0 
初期値問題とのアナロジーだと、真横ではなく下に投げ落とすイメージです。
実際に波動関数が折れ曲がる心配はありません。動径関数なので負の距離がありませんから。


一方、l≠0の時は遠心力が働いているのでゼロ距離にRはおらず(R=0＠ρ=0)
 下の図の右側のような格好にします。
最初のR=0、最初のR₀-R_-≠0にします。
(境界値ではなく初期値問題だと、地上から投げ上げるイメージ)
  

規格化

中心からの距離rをボーア半径a₀で割ったものが0から40くらいまでをプロットできれば実際問題としては十分なので
だいたいそのあたりまでのrmsを取ります。
刻み幅は意識しなくても構いません。ただ単に離散した数値の2乗平均ルート(rms)を取ればいいので
rms=sqrt(sumsq(塊)/count(塊))で求まります。
波動関数を2乗して存在確率を一旦出してやった上で、rms=sqrt(average(存在確率の塊))でも可です。

標準偏差stdevとの違いを意識したければ、rms²=stdev²(塊) +average²(塊)
 
を使ってもいいかもしれませんが、忘れた頃に何を計算したのか思い出せるかどうかは自己責任でお願いしますｗｗｗｗ
stdev関数があるならrms関数もあればいいのに・・・

このrmsで、算出した個々の波動関数の値を割れば、規格化完了です。

規格化と固有値の関係

 
量子数であるnやlが適正な値でないときは
遠距離での「規格化されてない」波動関数がバカでかい数値としてはじき出されます。
よって、規格化すると分母が大きいので波動関数全体がほとんどゼロになります。
つまり粒子が存在できないということです。
 

波動関数はnやlが本来整数のときに存在するのですが
差分法という手法や精度の粗さなどが原因なのか、シミュレーションでは整数からちょっとズレた値で波動関数が有限の値を取るようです。

nやlが変数となっている関数として波動関数のrmsを定義し、二分法(手動ｗｗｗｗ)を用いて関数rmsがほぼ零になるポイントを割り出しています。

二分法

 
二分法は工事中です。
なんかうまく収束しなかったので井戸型ポテンシャルのときのファイルを参照してみたんですがいまいちわからず。
もしかしたらz状にトレースしてる循環参照の脆弱性が現れたのかもしれません。

「rmsを最小にする」方法と、端っこを0にする方法の2種類があると思います。
ただ、rmsは0が最小なので、微分か差分しないと二分法には使えないのです。

また、緯度関数のほうは「端っこ」で評価できないみたいです。wikiカンペしましたすんません(・ω・｀；)

 

緯度

次に、緯度依存の波動関数を差分方程式に離散化するとこうなります。
 
動径R同様、P₊とP₀とP_-は差分化したときのフロントとセンターとバックです。

緯度θを変数変換してz=cosθとし、
関数もΘ(θ)からP(z)にしています。P(z)=Θ(θ)です。
要は、一度極座標にしたものを直交座標(円筒or円柱座標)に戻したような感じで、zというのは高さのことですね。(厳密には動径rで割った高さz/r)
lは先ほどと同様、方位量子数で、mは磁気量子数です。

境界条件の注意点

 
 緯度関数を解く際も初期値のような境界条件に注意してください。
方位量子数lと磁気量子数mとの和だか差の偶奇が、関数の偶奇に直結していますので(カンペ)

l+mを2で割ったあまりmod(l+m,2)が0のときは下図の左を(高所から真横に投げ落とす偶関数)
：最初のP≠0、最初のP'=0
あまりが1だった場合は（mod(l+m,2)=1）下図の右(地上から斜め上に投げ上げる奇関数)
：最初のP=0、最初のP'≠0

を、採用してください。

 

負の変数での関数値

 
zが負のときの計算値は、zが正の計算値を、符合を変えたり変えなかったりして流用してます。
-1の整数乗を作るのに、mod(l+m,2)の小数点以下をround関数で切り捨てています。-1の整数でない実数乗は基本的に多価になりますからね。
また、緯度θは-90度から90度までの180度が範囲です。360度ではありません。

疑問

 
規格化の際にrmsに√2をかけるとなぜかwikiと同じ値になりました。
理由はいまいちわかってませんが、360度じゃなくて半分の180度ってところと、存在確率が波動関数の2乗ってところが関係しているのかもしれません



それとDL用Excelファイルです。よかったら遊んでください
delおしっぱで再計算～the　animation～シリーズです。

デフォルトでは時刻取得関数で固有値一覧表を1つずつ移動していく状態になっています。
now()-today()に速さ50000だかをかけて動きが見えるようにし
round関数で整数にしてから、表の行数である20で割ったあまりをモジュロ関数で算出し
最後に1を足すことで0～19ではなく1～20にしています。

ここに手動で1～20の整数を入れると好きに固有値が見れます。

また、連続固有値探索モードとして、n、m、lそれぞれに循環参照を入れてもいいかと思います。
スイッチが0だったら初期値を返し、スイッチが1だったら増分だけ自分自身に加えるセルを用意してます。
初期値と増分を適当に調整しながら最小のrmsを探せ！(おもちゃのCMかな？)
割りと厳密に近づけないと、規格化された波動関数が有限にならないことがわかるかと思います。(特に動径のほう)



放送大学を見ると録画した他の回も見たくなって見ただけで俄然やる気がでてきます。
たった3回で水素の説明を終えてしまわれたのでだいぶ励みになりましたｗｗｗ

にほんブログ村