「物理学」カテゴリーアーカイブ

水素様原子の波動関数のコード

非相対論的な水素様原子の電子の波動関数を求めるfortranコードを書きます。

解く問題

電荷\(+Ze\)の原子核の周りに存在する、電荷\(-e\)を持つ電子が満たす時間依存しないシュレーディンガー方程式は、
\(
\displaystyle \left[-\frac{1}{2}\Delta-\frac{Z}{r}\right]\psi({\bf r})=E\psi({\bf r})
\)

です。
\(
\left. \psi({\bf r})\right|_{|{\bf r}|\to \infty}\to 0
\)

の境界条件を満たす正規直交化された解は、
\(
\begin{align}
\psi({\bf r})&=R_{n,l}(r)\cdot Y_{l,m}(\theta,\varphi)\\
&=\sqrt{\left(\frac{2}{n}\right)^3\frac{(n-l-1)!}{2n[(n+l)!]}} e^{-Zr/n} \cdot \left(2\frac{Zr}{n}\right)^l\cdot L_{n-l-1}^{2l+1}\left(2\frac{Zr}{n}\right)\cdot Y_{l,m}(\theta,\varphi)
\end{align}
\)

と表されます。ここで、

  • \(n,l,m\)  主量子数,方位量子数,磁気量子数
  • \(r,\theta,\varphi\)  動径方向,z軸方向の角度,xy平面の角度
  • \(L_n^{(\alpha)}(x)\)  ラゲール陪多項式[1]
  • \(Y_{l,m}(\theta,\varphi)\)  球面調和関数[2]

を意味します。

動径関数\(R_{n,l}(r)\)は微分方程式
\(
\displaystyle \left[-\frac{1}{2}\left(\frac{d^2}{dr^2}+\frac{2}{r}\frac{d}{dr}\right)+\frac{l(l+1)}{2r^2}-\frac{Z}{r}\right]R_{n,l}(r)=E_n R_{n,l}(r)
\)

を満たす正規直交化された関数で、固有値は
\(
\displaystyle E_n=-\frac{1}{2n^2}
\)
です。

波動関数を出力するFortranコードはこれ。

計算結果


3次元の波動関数を表示するには3つの座標と波動関数の値の計4つが必要です。
しかし、4次元空間を分かりやすく描画する方法はありません。

そこで、ここでは\(\mathrm{Re}{\psi(x,0,z,t)}\)のように、\(y=0\)にしてグラフを描画します。
\(z\)軸周りの構造がうまく表現できませんが、良しとします。

\(n=1\)

\(n=2\)

\(n=3\)

”color”と書いてあるのは、波動関数の形がはっきり映るようにカラーバーの範囲を調整したものです。
なので、色が飽和しているところはそのカラーバーの範囲で表現することはできない領域です。

\(n=20\)まで同じような画像を作ったのですが、画像の容量が大きくなりすぎました。
なので、見たい方は以下のリンクから元ファイルをダウンロードしてください。

\(n=1\)  Hatomc_n1.png(48KB)

\(n=2\)  Hatomc_n2.png(223KB)

\(n=3\)  Hatomc_n3.png(612KB)

\(n=4\)  Hatomc_n4.png(1MB)

\(n=5\)  Hatomc_n5.png(2MB)

\(n=6\)  Hatomc_n6.png(3MB)

\(n=7\)  Hatomc_n7.png(4MB)

\(n=8\)  Hatomc_n8.png(3MB)

\(n=9\)  Hatomc_n9.png(6MB)

\(n=10\)  Hatomc_n10.png(6MB)

\(n=11\)  Hatomc_n11.png(9MB)

\(n=12\)  Hatomc_n12.png(9MB)

\(n=13\)  Hatomc_n13.png(10MB)

\(n=14\)  Hatomc_n14.png(14MB)

\(n=15\)  Hatomc_n15.png(15MB)

\(n=16\)  Hatomc_n16.png(16MB)

\(n=17\)  Hatomc_n17.png(18MB)

\(n=18\)  Hatomc_n18.png(19MB)

\(n=19\)  Hatomc_n19.png(21MB)

\(n=20\)  Hatomc_n20.png(21MB)

n=1~n=20までのファイル hatom.zip(164MB)

例えば、n=13, n=20はこんな感じです(荒くしてます)。

特殊関数に関するプログラミングは特殊関数へ。

——————————————-

How to plot in gnuplot?
Just for note for myself.

You need 2 files which named : multi.plt and atommulti.plt

Before ploting, you must make directory “dat” and move every data files to “dat”.

and in gnuplot, type

call "multi.plt" 2

then, you can get figure.

multi.plt

atommulti.plt

参考


[1]ラゲール陪多項式\(L_m^{(\alpha)}(x)\)

  • 定義域

    \(\;\;\;\;
    0 \lt x \lt \infty
    \)

  • 微分方程式

    \(\;\;\;\;
    \displaystyle
    \left[x\frac{d^2}{dx^2}+(\alpha+1-x)\frac{d}{dx}+n\right]L_m^(\alpha)(x)=0
    \)

  • 具体例

    \(\;\;\;\;
    \begin{align}
    & L^{(0)}_0(x)=1 \\
    & L^{(0)}_1(x)=-x+1 \\
    & L^{(0)}_2(x)=\frac{1}{2}(x^2-4x+2) \\
    & L^{(0)}_3(x)=\frac{1}{6}(-x^3+9x^2-18x+6) \\
    & L^{(1)}_0(x)=1 \\
    & L^{(1)}_1(x)=-x+2 \\
    & L^{(1)}_2(x)=\frac{1}{2}\left(x^2-6x+6\right) \\
    & L^{(1)}_3(x)=\frac{1}{6}\left(-x^3+12x^2-36x+24\right)
    \end{align}
    \)

  • 漸化式

    \(\;\;\;\;
    \begin{align}
    L_n^{(\alpha)}(x)&=\binom{n+\alpha}{n}a_0(x) \\
    &\hspace{1em}a_{m-1}(x)=1-\frac{n-m+1}{m(\alpha+m)}\cdot x\cdot a_m(x) \\
    &\hspace{5em}(m=n,n-1,n-2,\cdots,1,\;\;a_n(x)=1)
    \end{align}
    \)
    [3]
    1階微分
    \(\;\;\;\;
    \begin{align}
    \frac{d}{dx}L_n^{(\alpha)}(x)=-L_{n-1}^{(\alpha+1)}(x)
    \end{align}
    \)
    2階微分
    \(\;\;\;\;
    \begin{align}
    \frac{d^2}{dx^2}L_n^{(\alpha)}(x)=L_{n-2}^{(\alpha+2)}(x)
    \end{align}
    \)

  • 直交性

    \(\;\;\;\;
    \displaystyle
    \int_{0}^{\infty} e^{-x}x^{\alpha}L^{(\alpha)}_n(x)L^{(\alpha)}_{n’}(x) dx =\frac{\Gamma(\alpha+n+1)}{n!}\delta_{n,n’},\;\;(\alpha>-1)
    \)

[2]球面調和関数\(Y_{l,m}(\theta,\varphi)\)
ルジャンドル多項式(\(P_l^m(x)\),[3])を用いて
\(
Y_{l,m}(\theta,\varphi)=(-1)^{(m+|m|)/2}\sqrt{\frac{2l+1}{4\pi}\frac{(l-|m|)!}{(l+|m|)!}}P_l^{|m|}(\cos\theta)e^{im\phi}
\)

[3]ルジャンドル陪多項式\(P_l^m(x)\)

  • 定義域

    \(\;\;\;\;
    1 \lt x \lt -1
    \)

  • 微分方程式

    \(\;\;\;\;
    \displaystyle (1-x^2)\frac{d^2}{dx^2}P_l^m(x)-2x\frac{d}{dx}P_l^m(x)+\left[l(l+1)-\frac{m^2}{1-x^2}P_l^m(x)\right]=0
    \)

  • 具体例

    \(\;\;\;\;
    \begin{align}
    & P_0^0(x)=1 \\
    & P_{1}^{-1}(x)=\frac{1}{2}(1-x^2)^{1/2} \\
    & P_{1}^{0}(x)=x \\
    & P_{1}^{1}(x)=-(1-x^2)^{1/2} \\
    & P_{2}^{-2}(x)=\frac{1}{8}(1-x^2) \\
    & P_{2}^{-1}(x)=\frac{1}{2}x(1-x^2)^{1/2} \\
    & P_{2}^{0}(x)=\frac{1}{2}(3x^2-1) \\
    & P_{2}^{1}(x)=-3(1-x^2)^{1/2} \\
    & P_{2}^{2}(x)=3(1-x^2)
    \end{align}
    \)

  • 漸化式

    用いている式
    \(\;\;\;\;
    \begin{align}
    & P_{|m|}^{|m|}(x)=(-1)^{|m|}(2|m|-1)!!(1-x^2)^{\frac{|m|}{2}} \\
    & P_{l}^{-m}(x)=(-1)^m\frac{(l-m)!}{(l+m)!}P_l^m(x) \\
    & P_{|m|}^{|m|+1}(x)=(2|m|+1)xP_{|m|}^{|m|}(x) \\
    & P_{|m|}^{|m|+q}(x)=\left(\frac{2|m|-1}{q}+2\right)xP_{|m|}^{|m|+q-1}(x)-\left(\frac{2|m|-1}{q}+1\right)P_{|m|}^{|m|+q-2}(x)
    \end{align}
    \)[1]
    1階微分
    \(\;\;\;\;
    \begin{align}
    m=\pm 1&\; \\
    &:\hspace{1em} \frac{d}{dx}P_l^m(x)=\frac{lx}{(x^2-1)}P_l^m(x)-\frac{l+m}{(x^2-1)} P_{l-1}^m(x) \\
    m\ne\pm 1\; \\
    &:\hspace{1em} \frac{d}{dx}P_l^m(x)=c_{2}P^{m+2}_{l-1}(x)+c_0 P^{m}_{l-1}(x)+c_{-2} P^{m-2}_{l-1}(x) \\
    &\hspace{1em} c_{2}=\frac{1}{4(m+1)}\\
    &\hspace{1em} c_{0}=\frac{l+m}{2}\left(1+\frac{l}{1-m^2}\right)\\
    &\hspace{1em} c_{-2}=-\frac{(l+m)(l+m-1)(l+m-2)(l-m+1)}{4(m-1)}
    \end{align}
    \)
         ※\(m=\pm 1\)の時、\(x=\pm 1\)で発散します。
    2階微分
    \(\;\;\;\;
    \begin{align}
    m=\pm 1,\pm 3&\; \\
    &:\hspace{1em} \frac{d^2}{dx^2}P_l^m(x)=\frac{1}{1-x^2}\left[2x\frac{dP_l^m(x)}{dx}+\frac{(l+1)(l+m)}{2l+1}\frac{dP_{l-1}^m(x)}{dx}-\frac{l(l-m+1)}{2l+1}\frac{dP_{l+1}^m(x)}{dx}\right] \\
    m\ne\pm 1,\pm 3&\; \\
    &:\hspace{1em} \frac{d^2}{dx^2}P_l^m(x)=c_{2}\frac{dP^{m+2}_{l-1}(x)}{dx}+c_0 \frac{dP^{m}_{l-1}(x)}{dx}+c_{-2} \frac{dP^{m-2}_{l-1}(x)}{dx}
    \end{align}
    \)
         ※\(m=\pm 1,\pm 3\)の時、\(x=\pm 1\)で発散します。

  • 直交性

    \(l\)に対する直交性
    \(\;\;\;\;
    \displaystyle
    \int_{-1}^{1} P_{m}^{l}(x)P_{m}^{l’}(x) dx =\frac{2}{2l+1}\cdot \frac{(l+m)!}{(l-m)!}\delta_{l,l’}
    \)

[3] Abramowitz & Stegun (1965), ISBN 978-0486612720 (Numerical Method, 22.18, Page 789)
[4] Abramowitz & Stegun (1965), ISBN 978-0486612720 (Numerical Method, 22.2, Page 774,775)

時間依存しないシュレーディンガー方程式と変分原理 2/2 (計算例)

計算例


具体的に適当なポテンシャルの形を考えて解いてみましょう。

プログラムは
時間依存しないシュレーディンガー方程式と変分原理 1/2
にあります。

おさらいですが、間接法による変分原理は、

  • 固有値の上限を与える(変分原理によるエネルギー固有値は必ず厳密なエネルギー固有値よりも大きくなる)
  • 基底関数の数を増やすほど精度は上がる(完全系の場合)

ということを想定しています。
そして、時間依存しないシュレーディンガー方程式と変分原理 1/2 (計算手法の説明)
では、これから束縛状態のみを考えるので、sin関数系の場合の変分原理について計算を行っていきました。
上に想定した変分原理の想定は、

正しい境界条件の下で誤差なく行列要素が求められ、誤差なく対角化が行われた

場合にのみ成り立ちます。sin関数系の場合の境界条件は固定端条件、つまり区間[a,b]で
\(
\psi(x=a)=\psi(x=b)=0
\)

ということです。

ここでは、解析的な解が存在する

  1. 井戸型ポテンシャル
  2. 三角量子井戸
  3. 調和ポテンシャル
  4. F.モースによるポテンシャル[3]
  5. \(1/(\cosh^2)\)型ポテンシャル[3]
  6. クーロンポテンシャル(動径方向, \(l=0\))

について比較していきます。

シュレーディンガー方程式は全て原子単位系(\(m=1, \hbar=1\))で考え、
\(
\displaystyle \left[-\frac{1}{2}\frac{d^2}{dx^2}+V(x)\right]\psi(x)=E\psi(x)
\)

の形を考えます。

無限に深い井戸型ポテンシャル


ポテンシャル

\(
V(x)=\left\{
\begin{eqnarray}
0 \hspace{2em}&(& 0\lt x \lt L,)\\
\infty \hspace{2em}&(& x\le 0 \ or\ L\le x)
\end{eqnarray}
\right.
\)

※井戸の区間が\([a,b]\)であっても変数変換\(y=x-a,\;\;L=b-a\)によって上記ポテンシャルの形に変換されます。

固有値

\(
\displaystyle E_n=\frac{(n+1)^2\pi^2}{2 L^2},\;\;(n=0,1,2,\cdots)
\)

固有関数

\(
\psi_n(x)=\left\{
\begin{eqnarray}
\sqrt{\frac{2}{L}}\sin\left(\frac{(n+1)\pi}{L}(x+\frac{L}{2})\right) &(& 0\lt x \lt L,)\\
0 \hspace{6em}&(& x\le 0 \ or\ L\le x)
\end{eqnarray}
\right.
\)

解法

区間\(0\lt x \lt L,\)ではシュレーディンガー方程式の\(V(x)=0\)の解、sin,cosの形をしてなければなりません。
また、\(\psi(x=0)=0,\;\psi(x=L)=0\)では境界条件より、波動関数の値はゼロです。
よって、解はsin,cosの形で、そのゼロ点がちょうど\(x=0,x=L\)にくるような関数となります。

数値計算解と解析解との比較

数値計算の設定:
Gauss-Lobatto求積法の次数(Ngl):12
計算区間(a,b):[0,\(\pi\)]
基底関数の数(N):200

エネルギー固有値

量子数\(n\) 数値計算解 解析解
0 0.5000000000000000 0.5000000000000000
1 2.0000000000000000 2.0000000000000000
2 4.5000000000000000 4.5000000000000000
3 8.0000000000000000 8.0000000000000000
4 12.5000000000000000 12.5000000000000000
10 60.500000000000000 60.500000000000000
20 220.50000000000000 220.50000000000000
30 480.50000000000000 480.50000000000000
40 840.50000000000000 840.50000000000000
50 1300.5000000000000 1300.5000000000000
考察

行列は元々対角です。なので、対角化する前から対角行列です。
誤差は入りようがありません。

三角量子井戸


ポテンシャル

\(
V(x)=\left\{
\begin{eqnarray}
\infty \hspace{2em}&(& x\lt 0)\\
\alpha x \hspace{2em}&(& 0 \lt x)
\end{eqnarray}
\right.
\)

固有値

\(
\displaystyle E_n=-\left(\frac{\alpha^2}{2}\right)^{1/3}a_n
\)

ここで、\(a_n\)はAiry関数[4]のゼロ点です。

固有関数

\(
\displaystyle \psi_n(x)=C\cdot \text{Ai}\left(\left(\frac{2}{\alpha^2}\right)^{1/3}(\alpha x-E_n)\right)
\)

※積分区間がエアリー関数のゼロ点から無限までであり、これは解析的に実行できません。なので規格化定数をあらわに書くことはできません。

解法

区間\(0\lt x \lt \infty,\)ではシュレーディンガー方程式の\(V(x)=\alpha x\)の解、Airy関数の形をしてなければなりません。
また、\(\psi(x=0)=0, \psi(x\to \infty)=0\)の境界条件を満たします。
よって、解はAiry関数の形で、\(x\)が漸近で減衰していく方の解(\(Ai(x),Bi(x)\)のAiの方)で、そのゼロ点がちょうど\(x=0\)にくるような関数となります。

数値計算解と解析解との比較

数値計算の設定:
Gauss-Lobatto求積法の次数(Ngl):12
計算区間(a,b):[0,60]
基底関数の数(N):200
\(\alpha\):1

エネルギー固有値

量子数\(n\) 数値計算解 解析解
0 1.8557572021449626 1.8557570814892386
1 3.2446077578572123 3.2446076240031596
2 4.3816713853293896 4.3816712392861303
3 5.3866139387958034 5.3866137807905003
4 6.3052631766282223 6.3052630065857747
10 10.8669422991860358 10.8669420487522821
20 16.8461591354108187 16.8461586901918032
30 21.8969186983512394 21.8969179157570224
40 26.4182318860569758 26.4182304793452509
50 30.5803380635310234 30.5803354172938597

考察

上記計算では基底状態から励起状態までほとんど精度が変わっていません。
また、変分原理が想定したとおり、全ての数値計算解は、厳密な解よりも高い結果を与えています。
厳密な解よりも高い結果を与えていることから、
行列要素は高精度で計算されているが、基底関数が足らない、もしくは計算区間が足らない
事が考えられます。

精度がなぜ良くないのかを考えましょう。
\(x=0\)の境界条件は\(\sin\)関数で良く表現されます。が、\(x\to\infty\)側は指数関数で減少していきます。
(一次の)指数関数で減衰していくことを\(\sin\)関数系ではうまく表現できないのではないか、と思います。
実際、Airy関数の指数関数での減衰は漸近で
\(
\displaystyle \frac{e^{-\frac{2}{3}x^{3/2}}}{x^{1/4}}
\)

と表されます。
この後、調和ポテンシャルとクーロンポテンシャルの場合も考えますが、
調和ポテンシャル(漸近で\(\exp(-x^2)\)、基底状態の精度14,15桁)

三角量子井戸(漸近でおおよそ\(\exp(-x^{2/3})\)、基底状態の精度7,8桁)

クーロンポテンシャル(漸近で\(\exp(-x)\)、基底状態の精度2,3桁)
の順に精度が悪くなっていきます。
指数関数的な減衰が早い≒無限に深い量子井戸
を意味するはずなので、恐らくこうなのではと思います。

調和ポテンシャル



\(
\displaystyle V(x)=\frac{1}{2}x^2
\)

固有値

\(
\displaystyle E_n=n+\frac{1}{2},\;\;(n=0,1,2,\cdots)
\)

固有関数

\(
\displaystyle \psi_n(x)=\frac{\pi^{-1/4}}{\sqrt{2^n n!}}e^{-\frac{x^2}{2}} H_n\left(x\right)
\)
ここで\(H_n(x)\)はエルミート多項式を表します。

解法

1次元調和振動子の直接解法

数値計算解と解析解との比較

数値計算の設定:
Gauss-Lobatto求積法の次数(Ngl):12
計算区間(a,b):[-30,30]
基底関数の数(N):200

量子数\(n\) 数値計算解 解析解
0 0.5000000000000038 0.5
1 1.5000000000000340 1.5
2 2.4999999999998268 2.5
3 3.4999999999996776 3.5
4 4.4999999999999165
10 10.4999999999999680 10.5
20 20.4999999999997691 20.5
30 30.5000000000065015 30.5
40 40.5000064794500432 40.5
50 50.5426649449273313 50.5
考察

変分原理の典型的な傾向を見せています。というのは高い励起状態になればなるほど精度が悪くなっているためです。
高励起状態になるにつれて精度が悪くなっていくのは基底関数が足らない為でしょう。計算区間に関してはポテンシャルは境界で高いまま(\(V(a)=V(b)=450\))であり、求められたエネルギー固有値の値\(E_{50}\approx 50\)と比べても十分大きいです。なので、計算区間は足りています。

また、変分原理の固有値の上限を与えるに関しては必ずしも満たしているわけでは無いことに気が付きます。これは行列要素を求める際の誤差や丸め誤差に起因するものであると考えられます。なので、原理が破られているのではなく、数値計算上の問題に起因しています。\(n=30,40,50\)では厳密解よりも大きな値を与えています。これは、数値計算上の誤差よりも変分原理の原理的な値の方が大きくなっているためです。

F.モースによるポテンシャル[3]



\(
\displaystyle V(x)=A\left(e^{-2\alpha x}-2e^{-\alpha x}\right)
\)

固有値

\(
\displaystyle E_n=-A\left[1-\frac{\alpha}{\sqrt{2A}}\left(n+\frac{1}{2}\right)\right]^2
\)

ここで、\(n\)は正の整数で、ゼロから始まり、
\(
\displaystyle\frac{\sqrt{2A}}{\alpha}\gt n+\frac{1}{2}
\)

を満足する最大値\(n_{\text{max}}\)まで。

固有関数

\(
\begin{align}
\psi(x)&=e^{-\xi/2} \xi^s w(\xi) \\
& w(\xi)=F(-n,2s+1,\xi) \\
& \xi=\frac{2\sqrt{2A}}{\alpha}e^{-\alpha x} \\
& s=\frac{\sqrt{-2E}}{\alpha}
\end{align}
\)
ここで\(F(a,b,x)\)は合流型超幾何関数です。

この問題では離散スペクトルは有限個だけ存在します。もしも
\(
\displaystyle\frac{\sqrt{2A}}{\alpha}\gt n+\frac{1}{2}
\)

であれば、離散スペクトルは一般に存在しません。

解法

[3]を参照してください

数値計算解と解析解との比較

数値計算の設定:
Gauss-Lobatto求積法の次数(NGL):12
計算区間(a,b):[-3,80]
基底関数の数(N):200
\(A\):5
\(\alpha\):1

エネルギー固有値

量子数\(n\) 数値計算解 解析解
0 3.5437173567132851 -3.5438611699158100
1 1.3814252908344455 -1.3815835097474309
2 0.2185797147047042 -0.2193058495790518
3 -0.0049226123802256
4 0.0081675032299395
10 0.0527783252864115
20 0.2666761210195411
30 0.6404251854103465
40 1.1703593933291645
50 1.8543552927405405
考察

厳密解は3つしか存在しません。この3つの状態に対しては数値計算解の方が厳密解よりも大きくなり、想定通りです。
しかし精度があまり良くありません。これはポテンシャルが\(V(x\to\infty)=0\)であり、波動関数が速やかにゼロに向かわないためだろうと考えられます。
また、基底状態でもあまり良い結果が得られていません。上記ポテンシャルの閉じ込めが強くなく、多くの基底関数が必要となり、数値計算で用いた200個では十分ではないのでしょう。

4つ目の解に注目しましょう。厳密解は存在しないのに数値計算ではあたかも解が存在するように見えます(\(E\le 0\))。また、変分原理は上限を与えるため、4つ目の解が存在すると考えてしまうかもしれません。しかし、これは違います。明らかにするために4つ目の解の波動関数の形(青色の線)を拡大してみますと、

となり、物理的に意味を成していないことが分かります。なぜなら、波動関数は\(x\to\infty\)に向かうにつれて徐々に小さくなっていかなければなりません。もしもこれが正しければ、\(x=80\)で節を持っている解ということになり、これは束縛状態の解としては適しません。

エネルギー固有値だけを考えては求められないものだったということですね。

さらに、束縛状態以上ではすべて連続状態です。それにもかかわらず\(E\gt 0\)でも離散的に固有値が得られているのは基底関数に固定端の境界条件を課しているためです。これらに物理的な意味は無いので評価する際には気を付けましょう。

\(1/(\cosh^2)\)型ポテンシャル[3]



\(
\displaystyle V(x)=-\frac{V_0}{\cosh^2 \alpha x}
\)

固有値

\(
\displaystyle E_n=-\frac{\alpha^2}{8}\left[-(1+2n)+\sqrt{1+\frac{8V_0}{\alpha^2}}\right]^2
\)

で、\(n\)は
\(
\displaystyle n\lt s = \frac{1}{2}\left(-1+\sqrt{1+\frac{8V_0}{\alpha^2}}\right)
\)

から有限個の準位が決まります。

固有関数

\(
\begin{align}
\psi(x)&=(1-\xi^2)^{\xi/2}F\left[\varepsilon-s, \varepsilon+s+1, \varepsilon+1, \frac{1-\xi}{2}\right]\\
& \varepsilon=\frac{\sqrt{-2E}}{\alpha}\\
& \xi=\tanh \alpha x
\end{align}
\)

ここでの\(F(a,b,c,x)\)は超幾何関数です。

解法

[3]を参照してください

数値計算解と解析解との比較

数値計算の設定:
Gauss-Lobatto求積法の次数(NGL):12
計算区間(a,b):[-30,30]
基底関数の数(N):200

エネルギー固有値

量子数\(n\) 数値計算解 解析解
0 3.8286053565505553 -3.885009803959261
1 1.9750293441823026 -1.9750294118777854
2 0.6459145040771800 -0.7050490197963089
3 0.0750649726143545 -0.0750686277148326
4 0.0076704714372625
10 0.1139873400292487
20 0.4849841979196881
30 1.1689850012563472
40 2.1388129318397651
50 3.3888397726738235
考察

エネルギー固有値の上限を確かに与えています。そのほかはおおよそF.モースによるポテンシャルと同じような結果を与えています。

精度に関して、エネルギー固有値の精度は\(n=0,2\)の時、1~2桁程度の精度ですが、\(n=1,3\)の時、は6~7桁一致と著しく精度が上がっています。計算区間を変えたりしたのですが、この傾向は変わりません。
なぜなのか…詳しくは分かりません。偶関数、奇関数の問題だと思います。すなわち、奇関数の場合はなにか\(x\gt 0\)の領域に誤差があっても同じ量の誤差で\(x\lt 0\)で打ち消し合ってくれるのですが、偶関数の場合は打ち消し合いが起きないのではないか、ということです。
ただ、これが行列要素の精度からくるものなのかは詳しく見ていないので分かりません…

クーロンポテンシャル(動径方向, \(l=0\))


ポテンシャル

\(
\displaystyle V(x)=-\frac{1}{x}
\)

※ここでは、一階微分の存在しない形にした時のシュレーディンガー方程式を考えています。

固有値

\(
\displaystyle E_n=-\frac{1}{2n^2},\;\;(n=1,2,\cdots,)
\)

です。

固有関数

\(
\begin{align}
\psi(x)=C \cdot e^{-x/n} x^l L_{n-1}^{(1)}(x)
\end{align}
\)

ここで,\(C\)は規格化定数, \(L_n^k(x)\)はラゲール倍多項式[5]で
\(
\begin{align}
L_0^k(x)&=1 \\
L_1^k(x)&=-x+k+1 \\
L_2^k(x)&=\frac{1}{2}[x^2-2(k+2)x+(k+1)(k+2)] \\
L_3^k(x)&=\frac{1}{6}[-x^3+3(k+3)x^2-3(k+2)(k+3)x+(k+1)(k+2)(k+3)]
\end{align}
\)
です。もちろん規格化定数をあらわに決めることが出来ますが、通常、この問題は3次元のシュレーディンガー方程式を解いた時の動径方向として出てくるので、その時の規格化は
\(
\displaystyle \int_0^\infty x^2 \psi(x)^2 dx =1
\)

で規格化されます。しかし、本稿の規格化は\(\displaystyle \int_0^\infty \psi(x)^2 dx =1\)
で規格化しているので、エネルギー固有値のみの比較を行います。

解法

\(x=0, x\to infty\)の漸近形を考えて、本当の解を
(\(x=0\)の漸近形)×(\(x=\infty\)の漸近形)×(未知関数)
と仮定します。これをシュレーディンガー方程式に代入すると未知関数がラゲール陪多項式だと分かります。

数値計算解と解析解との比較

数値計算の設定:
Gauss-Lobatto求積法の次数(NGL):12
計算区間(a,b):[0,80]
基底関数の数(N):200

エネルギー固有値

量子数\(n\) 数値計算解 解析解
1 0.4970775696080274 -0.5000000000000000
2 -0.1246297940063756 -0.1250000000000000
3 0.0554455545682885 -0.0555555555555556
4 0.0312035408553899 -0.0312500000000000
5 0.0199704293764349 -0.0200000000000000
6 0.0134298208377274 -0.0138888888888889
7 0.0068151982214300 -0.0102040816326531
8 0.0018616612074302 -0.0078125000000000
10 0.0251247065596211 -0.0050000000000000
20 0.2439421094228725 -0.0012500000000000
考察

クーロンポテンシャルの場合、\(x=\infty\)では波動関数は指数関数で減衰していきます。
このため、減衰がゆっくりになってしまい、sin関数でうまく表現できないために精度が落ちているのだと推測できます。
また、計算区間を十分に大きく取ると励起状態の表現はうまくいきますが、反対に基底状態の表現ができなくなります。なぜなら、多くの区間で波動関数はゼロですが、\(x=0\)の近傍だけに大きなピークを持つため、sin関数ではうまく表現できません。その兆候は上の数値計算でも見えてます。基底状態は2桁程度ですが、量子数の増加と共に1桁くらい精度が上がっています。
クーロンの場合、ラゲール関数を基底関数としてとる方が良いでしょう。

参考文献


[1]P. J. Davis, P. Rabinowitz著, 森正武訳, 「計算機による数値積分法」, Gauss-Lobatto求積法 p.88
[2]P. J. Davis, P. Rabinowitz著, 森正武訳, 「計算機による数値積分法」, ゼロ点に挟まれる区間の積分 p.131
[3]ランダウ=リフシッツ,「量子力学1」第3刷, p81-84
[4]Abramowitz and Stegun, HANDBOOK OF MATHEMATICAL FUNCTIONS, p.446,478 http://people.math.sfu.ca/~cbm/aands/page_446.htm

Airy関数のゼロ点は
エアリー関数 Ai, Bi(ゼロ点)-Keisan
より22桁の精度で求めたものを用いています。
ちなみにそれらの値(Airy関数のゼロ点)は、

0 -2.338107410459767038489
1 -4.087949444130970616637
2 -5.52055982809555105913
3 -6.78670809007175899878
4 -7.944133587120853123138

10 -13.6914890352107179283
20 -21.2248299436420969552
30 -27.58838780988244481195
40 -33.28488468190140187962
50 -38.52880830509424882263

です。ゼロ点を0,1,2,…と、今回の計算の都合上0からインデックスを付けています。
もしかすると最後の1~2桁はあってないかもしれませんが、本稿では16桁まであっていればいいので、問題ありません。

[5]
Associated Laguerre Polynomial -wolfram mathworld式(22)-(25)

時間依存しないシュレーディンガー方程式と変分原理 1/2 (計算手法の説明)

時間依存しないシュレーディンガー方程式を数値的に解くための、変分原理について述べていきます。
変分原理はエネルギー固有値の上限を与えます(厳密な解よりも変分原理で求めた解は必ず大きくなる)。
基底関数を適切に選べば高精度の結果が与えられます。ですが、基底関数をどのように選べばよいか?の指標は無く、経験則に頼らざるを得ない状況が出てきます。

ここでは、変分原理によってなぜエネルギー固有値が求まるのかを述べて行きます。

もしも変分原理より正しい解を得たければ、shooting methodと呼ばれる方法も推奨します。
恐らくこちらの方が色々なポテンシャルに対して有効ですが、
・計算速度が遅い
・発散の可能性がある
・変分原理よりも経験値が必要
です。

まとめ


区間[a,b]で定義された時間依存しないシュレーディンガー方程式
\(
\displaystyle \left[-\frac{1}{2}\frac{d^2}{dx^2}+V(x)\right]\psi(x)=E\psi(x)
\)

の解である固有関数を
\(
\displaystyle \psi(x)=\sum_n c_n \varphi_n(x) \ \ …(6)
\)

の形で固定端条件(\(\psi(a)=\psi(b)=0\))の下で探します。
ここで、\(\varphi_n(x)\)は区間\([a,b]\)で正規直交化されたsin関数系
\(
\displaystyle \varphi_n(x)=\sqrt{\frac{2}{b-a}}\sin\left(\frac{n\pi}{b-a}(x-a)\right),\;\;\;(n=1,2,\cdots,N)
\)

です。この時、行列要素\(H_{ij}=\langle\varphi_i|-\frac{1}{2}\frac{d^2}{dx^2}+V(x)|\varphi_j\rangle\)は
\(
\begin{align}
H_{ij}&=K_{ij}+V_{ij} \\
&\hspace{1em}K_{ij}=\frac{1}{2}\left(\frac{j\pi}{b-a}\right)^2\delta_{ij} \\
&\hspace{1em}V_{ij}=\frac{1}{j}\sum_{m=2}^{M-1}\omega_m^{(GL)} \cos\left(\frac{\pi}{2}x_m^{(GL)}\right)
\sum_{n=0}^{j-1}(-1)^n \sin\left(\frac{i\pi}{2j}(x_m^{GL}+2n+1)\right)V(X_n^{(m)}) \\
&\hspace{2em}X_{n}^{(m)}=\frac{h}{2}(x_m^{(GL)}+2n+1)+a,\;\;h=\frac{b-a}{j}
\end{align}
\)
と求められます(\(j\ge i\)を推奨)。ここで、\(\omega_m^{(GL)},\ x_m^{(GL)}\)はM点のGauss-Lobatto求積法の区間[-1,1]における重みと求積点で、
\(
\displaystyle \omega_m=\frac{2}{M(M-1)[P_{M-1}(x_m)]^2},\ \ x_m\ne \pm 1
\)

であり、求積点は関数
\(
\displaystyle \frac{P_{M-1}(x)+P_{M}(x)}{x-1}
\)

のゼロ点です。\(P_{n}(x)\)は\(n\)次Legendre多項式を表します。

  1. 解説
  2. 時間依存から時間依存しないシュレーディンガー方程式へ
  3. 間接法
    1. 行列要素の計算
    2. K_{nm}の計算
    3. V_{nm}の計算
    4. Gauss-Lobatto求積法による積分の計算
    5. ゼロ点に挟まれる区間の積分
  4. 数値計算プログラム
  5. 正規直交系を用いない場合
  6. 研究分野での数値計算法について
  7. 参考文献

計算の具体例は
時間依存しないシュレーディンガー方程式と変分原理 2/2 (計算例)をご覧ください。


スポンサーリンク



解説


対象とする問題を導出します。何はともあれ、時間依存シュレーディンガー方程式から出発します。

時間依存から時間依存しないシュレーディンガー方程式へ


時間依存シュレーディンガー方程式は
\(
\displaystyle i\frac{\partial \Psi(x,t)}{\partial t}=\hat{H}(x)\Psi(x,t)
\)

と記述されます。ここで、簡単ために原子単位系(\(\hbar=1, m=1\))を考えます。
変数分離の形で解を仮定します。
\(
\displaystyle \Psi(x,t)=\psi(x)U(t)
\)

代入して移行すると以下の式を得ます。

\(
\begin{align}
i\psi(x)\frac{\partial U(t)}{\partial t}=U(t)\hat{H}\psi(x) \\
i\frac{1}{U(t)}\frac{\partial U(t)}{\partial t}=\frac{1}{\psi(x)}\hat{H}\psi(x) = E
\end{align}
\)

2番目の式の右辺\(E\)は何らかの定数です。
左辺は時間のみを含む式、真ん中の式は位置のみを含む式なので、いつどんな時でもこの関係式が成り立っていないといけないという条件からこれらはなんらかの定数でなければなりません。その定数を\(E\)と置きました。

すなわち、2つの式
\(
\begin{align}
i\frac{\partial U(t)}{\partial t}=E U(t)\;\;\; \cdots(1)\\
\hat{H}\psi(x) = E\psi(x)\;\;\; \cdots(2)
\end{align}
\)
が導けます。

\(t\)に関する式は\(\hat{H}(x)\)に依らず解けてしまい、
\(U(t)=e^{-i E t}\)(1)
です。位置に関する式は
\(
\hat{H}\psi(x) = E\psi(x)
\)
(2)
を満たさなければなりません。この(2)が時間依存しないシュレーディンガー方程式と呼ばれているものです。
これを数値計算で解いていくことを考えます。

左から\(\psi^*\)を掛けて全空間で積分しますと、
\(
\begin{align}
\int \psi^* \hat{H} \psi dx – E \int \psi^* \psi dx = 0\;\;\; \cdots(3)
\end{align}
\)
という式を得ます。この問題は
条件\(\int \psi^* \psi dx=1\)の下で、\(\int \psi^* \hat{H} \psi dx\)の極値を求めよ、
という問いになります。
上記問題を考えた時、(3)と同一の式になります。

すなわち、原理的にはどんな\(\psi(x)\)を与えても\(E(=E[\psi])\)の値を計算できるわけです。
その無数に考えられる\(\psi\)の中で、特に(3)の解が停留する\(\psi\)、すなわち
\(
\displaystyle \delta \int \psi^* (\hat{H}-E) \psi dx = 0
\)
となる\(\psi\)を探せ、という問題になります。

↑がもしも停留する\(\psi\)が見付かれば、それはシュレーディンガー方程式と同一になるということは以下の通りに示せます。

\(
\displaystyle \delta \int \psi^* (\hat{H}-E) \psi dx =0 \;\;\;\cdots (4)
\)
厳密な\(\psi\)の変分に対して停留、すなわち\(\psi\to \psi+\delta \psi\)、もしくは複素平面で\(\frac{\pi}{2}\)ずれた方向\(\psi\to \psi+i\delta \psi\)に対して0でなければなりません。
(4)に対して2つの方向に対して変分を取ると、
\(
\begin{align}
\int (\delta \psi)^* (\hat{H}-E)\psi dx &+ \int \psi^* (\hat{H}-E)(\delta\psi) dx =0 \\
-i\int (\delta \psi)^* (\hat{H}-E)\psi dx &+ i\int \psi^* (\hat{H}-E)(\delta\psi) dx =0
\end{align}
\)
と2式が導けます。2式は変分をとった方向にかかわらず同時に満たされていなければなりません。同時に成立するためには
\(
(\hat{H}-E)\psi=0,\;\; \text{or}\;\; \psi^*(\hat{H}^{\dagger}-E^*)=0
\)
を満たしている必要があります。これは時間依存しないシュレーディンガー方程式と一致します。

以上をまとめますと、

厳密な\(\psi\)の変分はシュレーディンガー方程式に一致する

近似的に\(\psi\)を変分するとシュレーディンガー方程式の近似的な解に帰着されるだろう

近似的に\(\psi\)を構成し、変分問題を解く

代表的な方法として、
手法1)直接法…\(\psi\)の具体的な形を変分パラメータと共に与える
手法2)間接法…変分パラメータとしてある正規直交基底関数系を採用する
の2つがありますが、数値計算で使われる間接法を説明します。

間接法


間接法による変分原理は、

  • 固有値の上限を与える(変分原理によるエネルギー固有値は必ず厳密なエネルギー固有値よりも大きくなる)
  • 基底関数の数を増やすほど精度は上がる(完全系の場合)

という特性があります。

さて、今解きたい問題は
条件
\(
\displaystyle \int \psi^* \psi dx =1
\)

の下で
\(
\displaystyle \int \psi^* \hat{H} \psi dx
\)

の極値を与える\(\psi\)を探すこと、です。

すなわち
\(
\displaystyle \int \psi^* \hat{H} \psi dx – E \int \psi^* \psi dx =0 \ \ …(5)
\)

を考えて変分します。求めたい波動関数は
\(
\displaystyle \psi(x)=\sum_n c_n \varphi_n(x) \ \ …(6)
\)

と展開できる、と考えます。ここで、\(c_n\)は変分パラメータ、\(\varphi_n(x)\)は区間[a,b]で定義されている正規直交関数系です。正規直交関数系\(\{\varphi_n(x)\}\)は以下の性質を持ちます。
\(
\displaystyle \int_a^b \varphi_n(x) \varphi_m(x)dx =\delta_{nm}
\)

この制約が無くても変分は可能ですが、以降の式が簡略になるので、正規直交関数だけを考えます。

(6)を(5)に代入すると、
\(
\begin{align}
\sum_{n}\sum_{m} c_n^* c_m \left[H_{nm}-E \delta_{nm}\right]=0 \\
\end{align}\ \ …(7)
\)

なります。ここで、
\(
\begin{align}
H_{nm}&=\int_a^b \varphi_n(x) \hat{H} \varphi_m(x)dx
\end{align}
\)

とおきました。(7)式を満たす係数\(\{c_n\}\)の内、自明な解(\(c_n=0\ (for\ all\ n)\))はすぐ見つかります。しかし、この解は\(\psi(x)=0\)であり、条件\(
\displaystyle \int \psi^* \psi dx =1\)を満たしません。なので解としては不適切です。

では非自明な解を見つけましょう。(7)式の解を探すため、極値を探します。今、係数の組を変分パラメータとみなしているので、(変分パラメータで微分したもの)=0を計算します。すなわち
\(
\displaystyle \frac{\partial}{\partial c_n^*} (for\ all\ n)
\)

を(7)式に作用させます。すると、
\(
\begin{align}
\sum_{n}\sum_{m} \left[H_{nm}-E\delta_{nm}\right]c_m=0
\end{align}
\)

となります。\(n,m\)を行列の行、列ととらえれば、

\(
\begin{align}
\left(
\begin{array}{cccc}
H_{11} & H_{12} & \cdots & H_{1N}\\
H_{21} & H_{22} & \cdots & H_{2N}\\
\ddots & \ddots & \ldots & \ddots \\
H_{N1} & H_{N2} & \cdots & H_{NN}
\end{array}
\right)
\left(
\begin{array}{c}
c_{1} \\
c_{2} \\
\ddots \\
c_{N}
\end{array}
\right)
=
E
\left(
\begin{array}{c}
c_{1} \\
c_{2} \\
\ddots \\
c_{N}
\end{array}
\right)
\end{align}
\)

となるわけです。ということは、左辺の正方行列を対角化すれば右辺のE,すなわち固有値が求まる、という事になります。

行列要素の計算


行列を対角化するためには行列要素\(H_{nm}\)が分かっていないとなりません。
これからこの行列要素を計算していきます。

ハミルトニアンが
\(
\displaystyle \hat{H}=-\frac{1}{2}\frac{d^2 }{d x^2} + V(x)
\)

として表現されている場合、行列要素は
\(
\begin{align}
H_{nm}&=K_{nm}+V_{nm} \\
&\hspace{3em} K_{nm}=-\frac{1}{2}\int_a^b \varphi_n(x)\frac{d^2}{dx^2} \varphi_m(x)dx \\
&\hspace{3em} V_{nm}=\int_a^b \varphi_n(x)V(x) \varphi_m(x)dx
\end{align}
\)
と書き表され、これらを数値計算で求めます。

\(K_{nm}\)の計算


まず、微分を含んだ
\(
\displaystyle K_{nm}=-\frac{1}{2}\int_a^b \varphi_n(x)\frac{d^2}{dx^2} \varphi_m(x)dx
\)

の項を考えます。微分は数値計算で行いません。
というのも正規直交関数系はこちらで用意する知っている関数系なので、求めることが出来てしまうからです。
多くの正規直交関数系には漸化式と微分方程式が存在します。一番簡単な基底としてsin関数系を考えてみましょう。
区間\([a,b]\)で正規直交化されたsin関数系は
\(
\displaystyle \varphi_n(x)=\sqrt{\frac{2}{b-a}}\sin\left(\frac{n\pi}{b-a}(x-a)\right),\;\;\;(n=1,2,\cdots,N)
\)

と書き表せます。sin関数なので、2階微分はもちろん
\(
\displaystyle \frac{d^2}{dx^2}\varphi_n(x)=-\left(\frac{n\pi}{b-a}\right)^2\varphi_n(x)
\)

を満たします。これを積分の中に代入して、
\(
\begin{align}
K_{nm}&=-\frac{1}{2}\int_a^b \varphi_n(x)\frac{d^2}{dx^2} \varphi_m(x)dx \\
&=\frac{1}{2}\left(\frac{n\pi}{b-a}\right)^2\int_a^b \varphi_n(x) \varphi_m(x)dx \\
&=\frac{1}{2}\left(\frac{n\pi}{b-a}\right)^2\delta_nm
\end{align}
\)

となります。積分を計算することなく求める事が出来ました。

より一般には一回微分の存在しない形、もしくは漸化式によって手計算で微分を含む積分を片付けてしまったり、部分積分を用いて何とか計算します。
正規直交系であれば、うまく整理することによって誤差が全く入らず、かつ高速に計算できます。

\(V_{nm}\)の計算


さて、残る問題は
\(
\displaystyle V_{nm}=\int_a^b \varphi_n(x)V(x) \varphi_m(x)dx
\)

の計算です。基底関数の次数が低い場合、この計算は何も問題がありません。
調和ポテンシャル\(V(x)=x^2\)の時、\(n=1,m=2\)の時の被積分関数は

と非常になだらかですが、\(n=5,m=99\)の時の被積分関数は

となります。なので細かい分点をとることになり、非常に積分が難しいです。

離散フーリエ変換によって求められそうな気もしますが、周期境界条件が満たされていないことは明らかなのであきらめましょう。
有限区間で早く振動する関数の積分は難しいです。二重指数関数型数値積分も振動する関数は苦手です。

こういった場合の計算方法として、
ゼロ点の位置が分かっている場合、そこで区間を分割して積分する
という直観的な方法があります[2]。
ゼロ点の位置が分かっていれば、その間の被積分関数の振る舞いはシンプルな形になるので滑らかな関数として考えられます。

また、端点での値がゼロだと分かっている場合、Gauss-Lobatto求積法が非常に有用です。

Gauss-Lobatto求積法による積分の計算


Gauss-Lobatto求積法はGauss-Legendre求積法と同じく被積分関数を多項式で近似します[1]。
Gauss-Legendre求積法と異なる点は端点で分点を持つということです。端点を明示的に書けば、積分は
\(
\displaystyle \int_{-1}^1 f(x)dx\sim \frac{2}{N(N-1)}(f(-1)+f(1))+\sum_{i=2}^{N-1}\omega_i f(x_i)
\)

と書けます。ここで、\(\omega_i,x_i\)はGauss-Lobatto求積法における重み、求積点です。
Gauss-Lobatto求積法の重みは
\(
\displaystyle \omega_i=\frac{2}{N(N-1)[P_{N-1}(x_i)]^2},\ \ x_i\ne \pm 1
\)

であり、求積点は
\(
\displaystyle \frac{P_{N-1}(x)+P_{N}(x)}{x-1}
\)

のゼロ点です。\(P_{n}(x)\)は\(n\)次Legendre多項式を表します。

Gauss-Lobatto求積法は
\(f(\pm 1)\)が既知の値(特に0)を持つ場合
に便利なことが分かります。
関数\(f(\pm 1)\)がゼロであれば、その点を計算する必要性がなくなりますし、見かけ上の特異点であってもこちらで値を指定すれば参照されることはなく、計算できないという状況に陥りません。

ゼロ点に挟まれる区間の積分


ゼロ点に挟まれる区間の積分は、結局以下の式になります[2]。
\(
\displaystyle \varphi_j(x)=\sqrt{\frac{2}{b-a}}\sin\left(\frac{j\pi}{b-a}\right)
\)

とすると、積分は
\(
\begin{align}
&\int_a^b \varphi_i(x) V(x) \varphi_j(x) dx \\
&=\sum_{n=0}^{j-1} \int_{x_n}^{x_{n+1}} \varphi_i(x) V(x) \varphi_j(x)dx \\
&\approx \sum_{n=0}^{j-1}\sum_{m=2}^{M-1}\frac{h}{2} \omega_m^{(GL)}\varphi_i(X_n^{(m)}) V(X_n^{(m)}) \varphi_j(X_n^{(m)}) \\
&=\frac{h}{2}\sum_{m=2}^{M-1}\omega_m^{(GL)}\varphi_j(X_0^{(m)})S_j^{(m)} \\
\end{align}
\)
ここで、
\(
\begin{align}
&x_n=nh+a,\ h=\frac{b-a}{j},\ (n=0,1,\cdots,j) \\
&X_n^{(m)}=\frac{h}{2}x_m^{(GL)}+\frac{x_{n+1}+x_n}{2} \\
&X_{n+1}^{(m)}=X_n^{(m)}+h \\
&S_j^{(m)}=\sum_{n=0}^{j-1}(-1)^n \varphi_i(X_n^{(m)})V(X_n^{(m)})
\end{align}
\)
と置きました。
\(x_m^{(GL)},\ \omega_m^{(GL)}\)は区間\([-1,1]\)でのM点Gauss-Lobatto求積法の求積点と重みです。
計算を可能な限り減らすため、\(\varphi_j(x)\)の周期性を用いています。

数値計算プログラム


行列要素を組み立てた後は適当なルーチンを用いて対角化します。
今回の場合、基底関数をsin関数に選んだため、強制的に固定端条件が課されます。
なので、この基底関数系を用いるといかなる問題も端点\(a,b\)に無限に高いポテンシャルの壁がある、と読み替えることが出来ます。
よって行列は対称行列、エルミートです。故に固有値は必ず実数です。

この問題を解く場合、実対称行列を対角化するルーチン、lapackのdsyevdを用いるのが良いでしょう。
対角化のプログラムに関する詳細は対角化へどうぞ。

以下のプログラムは行列要素を上記手法で作り、対角化し、得られた波動関数を規格化し、出力するものです。

スポンサーリンク

コンパイルと実行は

$ ifort -mkl main.f90
$ ./a.out

などとしてください。lapackを用いますので、それと一緒にコンパイルしてください。

ラグランジュの未定乗数について


ラグランジュの未定乗数は特に物理的な意味を持たない魔法の変数、として置きますが、今回の場合は未定乗数がエネルギー固有値というちゃんとした意味を持つのです。

正規直交系を用いない場合


\(
\begin{align}
\sum_{n}\sum_{m} c_n^* c_m \left[H_{nm}-E I_{nm}\right]=0 \\
\sum_{n}\sum_{m} c_n^* c_m \left[J_{nm} H_{nm}-E \right]=0
\end{align}\ \ …(7)
\)

となります。ここで、
\(
\begin{align}
H_{nm}&=\int_a^b \varphi_n(x) \hat{H} \varphi_m(x)dx\\
J_{nm}&=I_{nm}^{-1}\\
I_{nm}&=\int_a^b \varphi_n(x) \varphi_m(x)dx
\end{align}
\)
と置きました。行列の形では、
\(
\begin{align}
\left(
\begin{array}{cccc}
J_{11}H_{11} & J_{12}H_{12} & \cdots & J_{1N}H_{1N}\\
J_{21}H_{21} & J_{22}H_{22} & \cdots & J_{2N}H_{2N}\\
\ddots & \ddots & \ldots & \ddots \\
J_{N1}H_{N1} & J_{N2}H_{N2} & \cdots & J_{NN}H_{NN}
\end{array}
\right)
\left(
\begin{array}{c}
c_{1} \\
c_{2} \\
\ddots \\
c_{N}
\end{array}
\right)
=
E
\left(
\begin{array}{cccc}
1 & 0 & \cdots & 0\\
0 & 1 & \cdots & 0\\
\ddots & \ddots & \ldots & \ddots \\
0 & 0 & \cdots & 1
\end{array}
\right)
\left(
\begin{array}{c}
c_{1} \\
c_{2} \\
\ddots \\
c_{N}
\end{array}
\right)
\end{align}
\)
と書けます。

\(
\begin{align}
\sum_{n}\sum_{m} c_n^* c_m \left[H_{nm}-E I_{nm}\right]=0 \\
\sum_{n}\sum_{m} c_n^* c_m \left[I_{nm}^{-1} H_{nm}-E \right]=0
\end{align}
\)

研究分野での数値計算法について


実際の研究の場では対角化はしますが、上記方法によって行列要素を作ることはしません。遅いからです。
特別な基底関数を用いると行列要素\(H_{nm}\)が1~2項、しかも非常に簡単な式に置き換えることが出来ます。精度も速度も上記方法以上です。
この計算方法に関してはこれ以上は述べることはしませんが、こういった方法も存在する、とだけ述べておきます。

参考文献


[1]P. J. Davis, P. Rabinowitz著, 森正武訳, 「計算機による数値積分法」, Gauss-Lobatto求積法 p.88
[2]P. J. Davis, P. Rabinowitz著, 森正武訳, 「計算機による数値積分法」, ゼロ点に挟まれる区間の積分 p.131
[3]ランダウ=リフシッツ,「量子力学1」第3刷, p81-84

水素原子の解析解1/3 -相対座標と変数分離

ボーアモデルにおける水素原子の解析解、すなわち電子1個(電荷-e)と陽子1個(電荷+e)の系でスピンを無視する場合を考えます。

戦略1)
この問題をそのまま解こうとすると電子で3次元、陽子で3次元なので、3+3=6次元の問題になります。
私達は今、クーロンポテンシャルが電子陽子間でしか依存しないことを知っています。
なので、古典力学で知られている通り、相対座標と重心座標を使い分離できそうだ、と考えられます。
やってみましょう。

2体問題における変数分離


2体系を考えます。(参考文献[1])
2つの物体(質量\(m_1, m_2\))には2体間の相対距離だけで決まる、時間に依存しないポテンシャル\(V({\bf r_1}-{\bf r_2})\)で相互作用しています。
よって、この系のハミルトニアン\(H\)は
\(
\displaystyle H=\frac{{\bf p_1}^2}{2m_1}+\frac{{\bf p_2}^2}{2m_2}+V({\bf r_1}-{\bf r_2})
\)

となるわけです。それに対応するシュレディンガー方程式は、
\(
\displaystyle i\hbar\frac{\partial}{\partial t}\Psi({\bf r_1},{\bf r_2},t)=\left[ -\frac{\hbar^2}{2m_1}\nabla^2_1-\frac{\hbar^2}{2m_2}\nabla^2_2 +V({\bf r_1}-{\bf r_2})\right]\Psi({\bf r_1},{\bf r_2},t)
\)

となります。ここで相対座標系
\(
\displaystyle {\bf r}={\bf r_1}-{\bf r_2},\;\;\; {\bf R}=\frac{m_1{\bf r_1}+m_2{\bf r_2}}{m_1+m_2}
\)

をとります。この座標系ではシュレディンガー方程式は以下のように書き直されます。
\(
\displaystyle i\hbar\frac{\partial}{\partial t}\Psi({\bf R},{\bf r},t)=\left[ -\frac{\hbar^2}{2M}\nabla^2_{\bf R}-\frac{\hbar^2}{2\mu}\nabla^2_{\bf r} +V({\bf r})\right]\Psi({\bf R},{\bf r},t)
\)

ここで\(M=m_1+m_2\)は系の全質量、\(\mu=\frac{m_1m_2}{m_1+m_2}\)は換算質量です。
今、ハミルトニアンが\({\bf R}と{\bf r}\)のそれぞれの和で書けています。
このことは変数分離が可能であることを示しています。なので波動関数を
\(
\Psi({\bf R},{\bf r},t)=\Phi({\bf R})e^{-iE_{CM}t/\hbar}\cdot\psi({\bf r})e^{-iEt/\hbar}
\)

として分離し、シュレディンガー方程式に代入、整理すれば\(\Phi({\bf R})\)と\(\psi({\bf r})\)に関する方程式
\(
\displaystyle -\frac{\hbar^2}{2M}\nabla^2_{\bf R}\Phi({\bf R})=E_{CM}\Phi({\bf R}) \;\;\;\; \ldots(A)
\)

\(
\displaystyle \left[-\frac{\hbar^2}{2\mu}\nabla^2_{\bf r} +V({\bf r})\right]\psi({\bf r})=E\psi({\bf r})\;\;\;\; \ldots(B)
\)

が導けます。
式\((A)\)はエネルギー\(E_{CM}\)の質量\(M\)の自由粒子に対するシュレディンガー方程式とみなせ、
式\((B)\)はポテンシャル\(V({\bf r})\)中を質量\(\mu\)を持った粒子に対するシュレディンガー方程式であるとみなせます。
系の全エネルギー\(E_{tot}\)は
\(
E_{tot}=E_{CM}+E
\)

です。
適切な座標変換により、2体問題2つの1体問題に置き換えた。そういう事を今やったのです。

——
水素原子に適応しましょう[2]。
クーロンポテンシャル\(V(r)\)は
\(
\displaystyle V(r)=-\frac{e^2}{4\pi\varepsilon_0 r}
\)

と書けます。ここで、\(r\)は電子陽子間距離を表します。
電子の質量を\(m\), 陽子の質量を\(M\)と記述します。

重心と共に動く座標系で考えましょう。
今、重心は止まっています(重心の運動量\({\bf P}=0\))。
なので、解くべき問題は1体の時間依存しないシュレディンガー方程式だけで、
\(
\displaystyle \left[-\frac{\hbar^2}{2\mu}\nabla^2 -\frac{e^2}{4\pi\varepsilon_0 r}\right]\psi({\bf r})=E\psi({\bf r}) \;\;\;\;\; \ldots(1)
\)

を解けばよい、ということになりました。\(\mu=\frac{mM}{m+M}\)はこの系の換算質量を表します。

戦略2)
6次元の問題から重心の自由度を3つ減らして3次元の問題に落とせたわけですが、3次元の問題を解くのはまだまだ大変です。相互作用ポテンシャルが中心力ポテンシャルであることから、3次元極座標で考えればrに関する項は変数分離できて、少なくともrに関する次元が更に1つ落とせそうです。あわよくば、残りの2次元も変数分離出来たら1次元の問題が3つになって嬉しいのですが、果たしてそううまくいくのでしょうか・・・

3次元極座標を用いた変数分離


ここでは3次元極座標、特に球面座標を使って変数分離を試みましょう(実は球面座標以外でも放物座標系と呼ばれる座標を使っても変数分離できます)

球面座標は、空間を\(r,\theta,\phi\)を用いて指定する座標です。デカルト座標との間には、

\(
\begin{eqnarray}
\left\{
\begin{aligned}
r &=\sqrt{x^2+y^2+z^2} \\
\theta&= \arccos{\frac{z}{r}}\\
\phi&=\tan^{-1}\frac{y}{x}
\end{aligned}
\right.
\;\;\;\;\;\;\;\;
\left\{
\begin{aligned}
x &= r\sin{\theta}\cos{\phi} \\
y &= r\sin{\theta}\sin{\phi} \\
z &= r\cos{\theta}
\end{aligned}
\right.
\end{eqnarray}
\)

の関係があり、体積要素\(d{\bf V}\)は
\(
d{\bf V}=r^2\sin{\theta}dr d\theta d\phi
\)

と表されます。

波動関数\(\psi({\bf r})\)を動径方向について分離します。
\(
\psi({\bf r})=R(r)Y(\theta,\phi)
\)

この座標系でのラプラシアン\(\nabla^2\)は
\(
\displaystyle \nabla^2=\frac{1}{r^2}\frac{\partial}{\partial r}\left(r^2\frac{\partial}{\partial r}\right)+\frac{1}{r^2\sin{\theta}}\frac{\partial}{\partial \theta}\left(\sin{\theta}\frac{\partial}{\partial \theta}\right)+\frac{1}{r^2\sin^2{\theta}}\frac{\partial^2}{\partial \phi^2}
\)

であるので、これらを式(1)に代入すると、
\(
\begin{align}
\displaystyle \left[-\frac{\hbar^2}{2\mu}\left\{\frac{1}{r^2}\frac{\partial}{\partial r}\left(r^2\frac{\partial}{\partial r}\right)+\frac{1}{r^2\sin{\theta}}\frac{\partial}{\partial \theta}\left(\sin{\theta}\frac{\partial}{\partial \theta}\right)+\frac{1}{r^2\sin^2{\theta}}\frac{\partial^2}{\partial \phi^2}\right\} -\frac{e^2}{4\pi\varepsilon_0 r}\right]R(r)Y(\theta,\phi)\\
=ER(r)Y(\theta,\phi)
\end{align}
\)

を得ます。


調べた限りではここからの解法は2つに分かれます。
①数学的に直接解く
②物理的意味を持つ角運動量演算子を考えながら解く
です。
どちらでも最後は同じになるので構わないと思います。
もしも深く知りたい場合は①を解いてから②を解くことをお勧めします。
長くなるのでページを移します。

①直接解く場合は後ほど

②物理的意味を踏まえながら解く場合は後ほど

参考文献

[1]Physics of Atoms and Molecules 2nd Edition pp.109~111 by B.H. Bransden C.J. Joachain(2003)
[2]Physics of Atoms and Molecules 2nd Edition pp.147~148 by B.H. Bransden C.J. Joachain(2003)
※参考した本はこれです。

exp(ikx)が進行波であることの証明

\(e^{i{\bf k\cdot r}}\)が\({\bf k}\)の方向に進行すること、\(e^{-i{\bf k\cdot r}}\)が\(-{\bf k}\)の方向に進行することの所以を証明します。

まとめ

波は、実部と虚部を並べた時、虚部がある方向に進行していきます。

証明


一次元自由粒子の時間発展を考えます。
自由粒子の時間依存シュレーディンガー方程式は
\(
\displaystyle i\hbar\frac{\partial}{\partial t}\Psi(x,t)=\left(-\frac{\hbar^2}{2m}\frac{\partial^2}{\partial x^2}\right)\Psi(x,t)
\)

と書けます。時間依存しないシュレーディンガー方程式を求めるために、
\(\Psi(x,t)=\psi(x)e^{-i\frac{E}{\hbar}t}\)
と置いて代入すると
\(
\displaystyle -\frac{\hbar^2}{2m}\frac{d^2}{dx^2}\psi(x)=E\psi(x)
\)

となります。
上の微分方程式の解は、
\(
\psi(x)=e^{\pm i{\sqrt{\frac{2mE}{\hbar^2}}x}}
\)

です。通常、波数\(k=\sqrt{\frac{2mE}{\hbar^2}}\)を定義して、
\(
\psi(x)=e^{\pm ikx}
\)

と書きます。連続状態を考えているので、エネルギー\(E\)は、
\(
\displaystyle E=\frac{\hbar^2k^2}{2m}
\)

です。

時間依存しないシュレーディンガー方程式が解けたので、時間依存するシュレーディンガー方程式の解\(\Psi(x,t)\)は、角周波数\(\omega=\frac{E}{\hbar}\)を用いて、
\(
\begin{align}
\Psi(x,t)&=e^{\pm ikx}e^{-i\frac{E}{\hbar}t} \\
&=e^{i(\pm kx-\omega t)}
\end{align}
\)

と記述されます。

\(e^{ikx}\)の解について


時間依存しない解が
\(e^{ikx}\)
に対する時間発展を考えます。
この時、時間依存する解は\(e^{i(kx-\omega t)}\)ですが、変形して、
\(
\displaystyle e^{ik(x-\frac{\omega}{k} t)}
\)

と書きます。
時刻\(t=0\)に対して、波動関数は\(e^{ikx}\)
時刻\(t=\Delta t\)に対して、波動関数は\(e^{ik(x-\frac{\omega}{k}\Delta t)}\)
です。よって、時刻\(t=0\to \Delta t\)の間に、関数が\(f(x)\to f(x-a)\)になった、つまり位置\(+a\)だけ平行移動された関数になった、とみることができます。
なので進行する波であると言えるのです。
また、その時の進行速度を考えますと、時間\(\Delta t\)の間に距離\(\frac{\omega}{k}\Delta t\)だけ進んだわけですから、速度\(v\)は、
\(
\displaystyle v=\frac{\frac{\omega}{k}\Delta t}{\Delta t}=\frac{\omega}{k}
\)

です。
グラフで時間経過を表せば、以下のようになります。
時刻t=0で\(e^{ikx}\)で表される波形(実線、実部:紫線、虚部:緑線)が時間経過して\(e^{ik(x-\frac{\omega}{k}\Delta t)}\)となると破線になります。
expikx_c
↑画像のgnuplotスクリプト

\(e^{-ikx}\)の解について


次に時間依存しない解が
\(e^{-ikx}\)
に対する時間発展を考えます。
この時、時間依存する解は\(e^{-i(kx+\omega t)}\)ですが、変形して、
\(
\displaystyle e^{-ik(x+\frac{\omega}{k} t)}
\)

と書きます。
時刻\(t=0\)に対して、波動関数は\(e^{-ikx}\)
時刻\(t=\Delta t\)に対して、波動関数は\(e^{-ik(x+\frac{\omega}{k}\Delta t)}\)
です。よって、時刻\(t=0\to \Delta t\)の間に、関数が\(f(x)\to f(x+a)\)になった、つまり位置\(-a\)だけ平行移動された関数になった、とみることができます。
なので後進する波であると言えるのです。
また、その時の進行速度を考えますと、時間\(\Delta t\)の間に距離\(-\frac{\omega}{k}\Delta t\)だけ進んだわけですから、速度\(v\)は、
\(
\displaystyle v=\frac{-\frac{\omega}{k}\Delta t}{\Delta t}=-\frac{\omega}{k}
\)

です。
グラフで時間経過を表せば、以下のようになります。
時刻t=0で\(e^{-ikx}\)で表される波形(実線、実部:紫線、虚部:緑線)が時間経過して\(e^{-ik(x+\frac{\omega}{k}\Delta t)}\)となると破線になります。
exp-ikx_c
↑画像のgnuplotスクリプト


波動関数は、
進行する波は虚部が先に来て後から実部(実部の位相が\(\pi/2\)遅れる)
後進する波は実部が先に来て後から虚部(虚部の位相が\(\pi/2\)遅れる)
のように見えます。
言い換えれば、波は虚部が先行している方向に進む、ということです。

3次元の場合


3次元でも同じことが言えます。
3次元自由粒子のハミルトニアンは
\(
\displaystyle \hat{H}=-\frac{\hbar^2}{2m}\left(\frac{\partial^2}{\partial x^2}+\frac{\partial^2}{\partial y^2}+\frac{\partial^2}{\partial z^2}\right)
\)

で変数分離可能なので、3次元の1つの方程式ではなく、1次元の3つの方程式に変換可能であり、x,y,z方向は独立に決められます。
なので、
\(
\displaystyle \Psi(x,y,z,t)=e^{ik_x(x-\frac{\omega_x}{k_x} t)+ik_y(y-\frac{\omega_y}{k_y} t)+ik_z(z-\frac{\omega_z}{k_z} t)}
\)

と表されます。
ここで\(\omega_x, \omega_y, \omega_z\)は、1次元の3つの方程式に分けた時のそれぞれのエネルギー固有値\(E_x,E_y,E_z\)に対応する角周波数を示しています(\(\omega_i=E_i/\hbar,\;\; (i=x,y,z)\))。(エネルギーは\(E=E_x+E_y+E_z\)です。)

上の場合は\({\bf k}=(k_x,k_y,k_z)\)方向に進行する波を表します。
通常は波数ベクトルとして書いて、
\(
\begin{align}
& \Psi({\bf r},t)= e^{i{\bf k \cdot r}-\omega t} \\
&\;\;\;\;\; {\small (\omega=\omega_x+\omega_y+\omega_z)}
\end{align}
\)

と書かれます。

1次元調和振動子の直接的解法

1次元調和振動子の良く知られた解法は3種類あります。

  1. 昇降演算子を利用し\(\hat{a}\psi(x)=0\)から求める方法
  2. 微分方程式を境界条件の下、直接解く方法
  3. ハイゼンベルグ方程式から解く方法

です。ここでは直接的に解く2番目の方法を紹介します。

対象にする微分方程式は一次元調和振動子の原子単位系での表現
\(
\displaystyle \left[-\frac{1}{2}\frac{d^2}{dx^2}+Ax^2\right]\psi(x)=E\psi(x)
\)

です。ここでAは正の実数です。
この微分方程式を境界条件\(\left. \psi(x)\right|_{x\to \pm\infty}=0\)の下で解きます。

まとめますと、エネルギー\(E_n\)と固有関数\(\psi_n(x)\)は、
\(
\begin{align}
E_n&=\sqrt{2A}\left(n+\frac{1}{2}\right)\\
\psi_n(x)&=\frac{\left(2A\right)^{1/8}}{\pi^{1/4}\sqrt{2^n n!}}H_n( {\scriptsize (2A)^{1/4}}x)\exp\left(-\sqrt{\frac{A}{2}} x^2\right)\;\;\;\; (n=0,1,2,\cdots)
\end{align}
\)

です。


導出


両辺を\(-2\)を掛けて式変形して
\(
\displaystyle \left[\frac{d^2}{dx^2}-2Ax^2+2E\right]\psi(x)=0
\)

常套手段ですが、この微分方程式を満たすべき関数\(\psi(x)\)の漸近形を考えます。
何とか解ける部分だけを解いていこうという方針です。

漸近での解


もしも\(x\to \infty\)だったならば、上式の\(2E\)の項は定数であるため、\(-2Ax^2\)と比べて無視できるようになるはずです。
ただし、微分の項\(\frac{d^2}{dx^2}\)は比較することはできないため無視はしないで置いておきます。
つまり、関数\(\psi(x)\)は\(x\to \infty\)の漸近で
\(
\displaystyle \left[\frac{d^2}{dx^2}-2Ax^2\right]\psi(x)=0
\)

の形の微分方程式を満足しなければなりません。
式変形して
\(
\displaystyle \frac{d^2}{dx^2}\psi(x)=2Ax^2\psi(x)
\)

関数を2回微分したら\(x^2\times\)(自分自身)が出てくる、という性質を持つ関数は
\(
\displaystyle \psi(x)= C \exp(Bx^2)
\)

という形で掛けそうです。実際にこの形の2階微分は
\(
\begin{align}
\frac{d\psi}{dx}&=C\cdot 2Bx \exp(Bx^2) \\
\frac{d^2\psi}{dx^2}&=C\cdot 2B(1+2Bx^2)\exp(Bx^2) \\
&\sim C\cdot 4B^2x^2\exp(Bx^2)
\end{align}
\)
です。今、2階微分で括弧内に定数”1”が出てきていますが、\(x\to \infty\)の漸近を考えているのでこの定数は\(2Bx^2\)と比較して無視されるはずです。よって最後の式変形をしています。

漸近での微分方程式に代入すると、
\(
\begin{align}
\frac{d^2}{dx^2}\psi(x)&=2Ax^2\psi(x) \\
C 4B^2x^2\exp(Bx^2)&=2Ax^2 C \exp(Bx^2) \\
2C(2B^2-A)x^2\exp(Bx^2)&=0
\end{align}
\)
上記の方程式は\(x\)の値に依らず常に満たされなければなりません。
よって、\(2C=0\)もしくは\(2B^2-A=0\)である必要があります。
\(C=0\)の解は波動関数が至る所でゼロである自明な解です。なので物理的には適しません。

※漸近でゼロであるだけで漸近じゃないところではokなんじゃないの?という疑問が出てきますが、\(C=0\)の場合、漸近で\(\psi(x)=0\)であり、その微分\(\frac{d\psi}{dx}=C 2Bx \exp(Bx^2)\)もゼロです。波動関数は全領域に渡って滑らかに続いていなければならないため、いつまでたっても波動関数はゼロです。ということは、存在確率がゼロ、有限にはならないつまらない解なのです。通常、\(x\to \infty\)の漸近で\(\psi(x)=0\)という条件は入れますが、その微分\(\frac{d\psi}{dx}\)もゼロという条件は入れません。

なので物理的に意味のある解は
\(
\begin{align}
2B^2-A=0 \\
\to B=\pm \sqrt{\frac{A}{2}}
\end{align}
\)
というBでなければなりません。よって波動関数\(\psi(x)\)は漸近で2つの解の線形結合として書かれ、
\(
\displaystyle \psi(x)= C_1 \exp\left(-\sqrt{\frac{A}{2}} x^2\right)+C_2 \exp\left(\sqrt{\frac{A}{2}}x^2\right)
\)

という形で記述できます。
しかし、境界条件\(\left. \psi(x)\right|_{x\to \pm\infty}=0\)を満たさなければならない場合、係数\(C_2=0\)である必要があります。よって波動関数の漸近形は
\(
\displaystyle \psi(x)= C\exp\left(-\sqrt{\frac{A}{2}} x^2\right)
\)

と求められます。

漸近ではないところの解


漸近では波動関数は
\(
\displaystyle \psi(x)= C\exp\left(-\sqrt{\frac{A}{2}} x^2\right)
\)

という形を満たさなければならないことが分かりました。残りは漸近ではない場所の解です。
定数変化法を使います。定数をxに依存する関数として全領域に渡る波動関数を
\(
\displaystyle \psi(x)= f(x)\exp\left(-\sqrt{\frac{A}{2}} x^2\right)
\)

と仮定します。この操作は\(\psi(x)\)に関する方程式から\(f(x)\)に関する微分方程式に焼直す操作です。
全領域に渡る波動関数の満たすべき微分方程式は
\(
\displaystyle \left[\frac{d^2}{dx^2}-2Ax^2+2E\right]\psi(x)=0
\)

なので、これに代入します。
計算過程で出てくる文字を減らすために、\(B=-\sqrt{\frac{A}{2}}\)として、\(B\)で記述していきます。
\(
\begin{align}
&\left[\frac{d^2}{dx^2}-2Ax^2+2E\right]\psi(x)= \\
&\;\;\;\;\;\;\;\left\{\frac{d^2f}{dx^2}+4Bx\frac{df}{dx}+(4B^2-2A)x^2+(2B+2E)f\right\}e^{Bx^2}=0
\end{align}
\)
と求める事ができ、\((4B^2-2A)=0\)を利用すると、\(f(x)\)の満たすべき微分方程式は、
\(
\displaystyle \frac{d^2f}{dx^2}+4Bx\frac{df}{dx}+(2B+2E)f=0
\)

と求まります。
この微分方程式の解はエルミートの微分方程式と呼ばれる形をしています。
エルミートの微分方程式は
\(
\displaystyle \frac{d^2 H(x)}{dx^2}-2x\frac{dH(x)}{dx}+2nH(x)=0
\)

の微分方程式を満たす直交多項式です[1]。
そのままエルミートの微分方程式を今回の問題に適応することはできません。なぜならば1階微分の係数が違っているからです。
なので変数変換により、\(x=\alpha y\)とおいて、一階微分の係数\(4B\)をちょうど\(-2\)にする定数\(\alpha\)を探しましょう。
\(
\begin{align}
\frac{d}{dx}&=\frac{dy}{dx}\frac{d}{dy}=\frac{1}{\alpha}\frac{d}{dy}\\
\frac{d^2}{dx^2}&=\frac{1}{\alpha^2}\frac{d^2}{dy^2}
\end{align}
\)
を代入して、
\(
\displaystyle \frac{1}{\alpha^2}\frac{d^2 f(x)}{dy^2}+4By\frac{df(x)}{dy}+(2B+2E)f(x)=0
\)

xとyが混在しているので気持ち悪いです。\(f(x)\to f(y)\)とおきましょう。そして最後に求められたyを用いたf(y)をxの式f(x)に焼直します。
両辺に\(\alpha^2\)を掛けると
\(
\displaystyle \frac{d^2 f(y)}{dy^2}+4B\alpha^2 y\frac{df(y)}{dy}+(2B+2E)\alpha^2 f(y)=0
\)

です。今、一階微分の係数を\(-2\)にする定数\(\alpha\)を探していました。
なので、
\(
\begin{align}
4B\alpha^2&=-2 \\
\alpha&=\pm\sqrt{-\frac{1}{2B}}=\pm\left(\frac{1}{2A}\right)^{1/4}
\end{align}
\)
となります。ここの\(\alpha\)のプラスマイナスはどちらでもいいです。
プラスに取ることにすれば、\(f(y)\)についての微分方程式は、
\(
\displaystyle \frac{d^2 f(y)}{dy^2}-2 y\frac{df(y)}{dy}-\left(\frac{B+E}{B}\right) f(y)=0
\)

であるため、エルミートの微分方程式の形になりました。
エルミートの微分方程式は\(f(y)\)に掛かる項、すなわち\(-\left(\frac{B+E}{B}\right)\)の部分が\(2n\), (\(n=0,1,2,\cdots\))である時だけ解を満たします。
よって、
\(
\begin{align}
-\frac{B+E}{B}&=2n\\
E=E_n&=-2B(n+\frac{1}{2})\\
&=\sqrt{2A}(n+\frac{1}{2})\;\;\; (n=0,1,2,\cdots)
\end{align}
\)
を満たすエネルギー\(E=E_n(n=0,1,2,\cdots)\)でなければ解は存在しません。
その時の解\(f(y)\)も\(n\)で番号付けされて、
\(
f(y)=H_n(y)\;\;\; (n=0,1,2,\cdots)
\)

と記述され、右辺はエルミート多項式と呼ばれます。
エルミート多項式\(H_n(y)\)は
\(
\begin{align}
H_0(y)&= 1\\
H_1(y)&= 2y\\
H_2(y)&= 4y^2-2\\
H_3(y)&= 8y^3-12y\\
H_4(y)&= 16y^4-48y^2+12\\
H_5(y)&= 32y^5-160y^3+120y\\
H_6(y)&= 64y^6-480y^4+720y^2-120\\
\ldots
\end{align}
\)
と言う関数です。三項漸化式
\(
H_n(y)=2yH_{n-1}(y)-2(n-1)H_{n-2}(y)
\)

を用いて\(n=2\)以上のエルミート多項式は記述されます。

以上より微分方程式を満たす解とエネルギーは、整数(\(n=0,1,2,\cdots\))で番号付けされた解を持つことが分かりました。
\(y\)を変数変換に従い\(x\)に戻します。
\(
\begin{align}
x&=\alpha y \\
&\alpha=\left(\frac{1}{2A}\right)^{1/4}
\end{align}
\)
なので、
\(H_n(y)=H_n(x/\alpha)\)
です。
漸近形での波動関数と合わせて、本当に求めたかった解\(\psi(x)\)は、
\(
\begin{align}
\displaystyle \psi(x) &= f(x)\exp\left(-\sqrt{\frac{A}{2}} x^2\right) \\
&=C H_n(x/\alpha)\exp\left(-\sqrt{\frac{A}{2}} x^2\right)
\end{align}
\)

であることが分かりました。定数\(C\)の自由度がありますが、これは物理的に合った条件を用いてなされます。

おさらい


解きたかった問題は
\(
\displaystyle \left[-\frac{1}{2}\frac{d^2}{dx^2}+Ax^2\right]\psi(x)=E\psi(x)
\)

を境界条件\(\left. \psi(x)\right|_{x\to \pm\infty}=0\)の下で解くことでした(A>0)。

この微分方程式を何とか解くために、まず\(x\to \pm\infty\)で解の満たすべき形を考えました。
その後、全領域で境界条件を満たす解を定数変化法によって求めました。
計算の結果、微分方程式の解は与えられた適当なエネルギー\(E\)に対していつも存在するわけではなく、特別な\(E=E_n\)で番号付けされた値でなければ解が存在しないことが分かりました。
\(n\)番目のエネルギー\(E_n\)と関数\(\psi_n(x)\)は、
\(
\begin{align}
E_n&=\sqrt{2A}\left(n+\frac{1}{2}\right) \\
&\psi_n(x)= C_n H_n(x/\alpha)\exp\left(-\sqrt{\frac{A}{2}} x^2\right) \\
&\;\;\;\;\alpha=\left(2A\right)^{-1/4}
\end{align}
\)

です。与えられた条件下では規格化までは追求していないのでここでお仕舞いですが、規格化を考えてみましょう。

規格化


波動関数\(\psi_n(x)\)を規格直交化(正規直交化とも言う)します。
規格化は、全空間での存在確率
\(
\displaystyle \int_{-\infty}^{\infty}|\psi(x)|^2dx
\)

を1にすること。
直交化は異なる状態間をかけ合わせて積分した時ゼロになること・・・すなわち内積をゼロにします。
\(
\displaystyle \int_{-\infty}^{\infty}\psi_m^*(x)\psi_n(x)dx=C\cdot\delta_{mn}
\)

ここで\(C\)は任意の定数で、\(\delta_{mn}\)はクロネッカーのデルタを表します。
両方を考慮すると、規格直交化とは、
\(
\displaystyle \int_{-\infty}^{\infty}\psi_m^*(x)\psi_n(x)dx=\delta_{mn}
\)

を満たすように係数を決定します。

実際に波動関数を代入すると、
\(
\begin{align}
\int_{-\infty}^{\infty}\psi_m^*(x)\psi_n(x)dx&= C_m^*C_n \int_{-\infty}^{\infty} H_m(x/\alpha)H_n(x/\alpha)\exp\left(-\sqrt{2A} x^2\right)dx
\end{align}
\)

変数変換を行います。
\(
\begin{align}
x&=\alpha y \\
&\alpha=\left(\frac{1}{2A}\right)^{1/4}
\end{align}
\)
より、
\(
\begin{align}
& C_m^*C_n \int_{-\infty}^{\infty} H_m(x/\alpha)H_n(x/\alpha)\exp\left(-\sqrt{2A} x^2\right)dx \\
&= C_m^*C_n \int_{-\infty}^{\infty} H_m(y)H_n(y)\exp\left(-\sqrt{2A}\alpha^2y^2\right)\cdot \alpha dy \\
&= C_m^*C_n\alpha \int_{-\infty}^{\infty} H_m(y)H_n(y)\exp\left(-y^2\right) dy \\
\end{align}
\)
ここで、エルミート多項式の特性を使います。エルミート多項式の関係式の中に、
\(
\displaystyle \int_{-\infty}^{\infty} H_m(x)H_n(x)\exp\left(-y^2\right) dx=2^n n!\sqrt{\pi}\delta_{mn}
\)

という関係式があります。
なので、これを使うと、
\(
\displaystyle \int_{-\infty}^{\infty}\psi_m^*(x)\psi_n(x)dx=C_m^*C_n\alpha 2^n n!\sqrt{\pi}\delta_{mn}
\)

と求められます。今、\(({\mbox 右辺})=\delta_{mn}\)にしたいので、係数について、
\(
\begin{align}
C_m^*C_n\alpha 2^n n!\sqrt{\pi}\delta_{mn}&=\delta_{mn} \\
(C_m^*C_n\alpha 2^n n!\sqrt{\pi}-1)\delta_{mn}&=0
\end{align}
\)
にするように\(C_n\)を決定すればよいことになります。
クロネッカーのデルタより、\(m=n\)以外の時は係数に依らずゼロになるため、自明です。
なので、\(m=n\)の時,
\(
\begin{align}
|C_n|^2\alpha 2^n n!\sqrt{\pi}&=1\\
\to |C_n|&=\sqrt{\frac{1}{\sqrt{\pi}2^n n!\alpha}} \\
C_n&=\frac{1}{\sqrt{\alpha}\pi^{1/4}\sqrt{2^n n!}}\cdot \exp(i\theta)
\end{align}
\)
と規格化定数が求められます。位相\(\exp(i\theta)\)の任意性がありますが、この位相に特に意味は無く、自由に選べます。よって、式が簡単になる\(\theta=0\)を取ることにします。
以上の結果から、規格直交化された波動関数は\(\psi_n(x)\)は、
\(
\begin{align}
\psi_n(x)&=\frac{1}{\sqrt{\alpha}\pi^{1/4}\sqrt{2^n n!}}H_n(x/\alpha)\exp\left(-\sqrt{\frac{A}{2}} x^2\right)\\
\alpha&=\left(2A\right)^{-1/4}
\end{align}
\)
となります。\(\alpha\)を使わず、\(A\)だけで表せば、
\(
\displaystyle \psi_n(x)=\frac{\left(2A\right)^{1/8}}{\pi^{1/4}\sqrt{2^n n!}}H_n( {\scriptsize (2A)^{1/4}}x)\exp\left(-\sqrt{\frac{A}{2}} x^2\right)
\)

です。

参考文献

[1]小野寺嘉孝著「物理のための応用数学」裳華房 p.71

バレル内部でのBB弾の方程式

バレル内部でのBB弾の運動方程式です。

目的は、

BB弾は何秒間バレル内部に存在しているのか?
バレルが長いと減速になりうるのか?

を知ることです。
ただし、今は時間が無いため運動方程式を提示する、にとどめます。

弾道計算(BB弾)の理論
弾道計算(BB弾)の結果
弾道計算の結果2, 比較と詳細データ
弾道計算(BB弾)のコード(fortran90)
バレル内部でのBB弾の方程式←今ここ


電動ガンの場合です。
バレル内部1
空気抵抗は初速0[m/s]から90[m/s]前後にまで加速されるわけですから、空気抵抗の粘性抵抗、慣性抵抗どちらの項も無視することはできないでしょう。

いくつか仮定があります。

  1. 気体の状態方程式\(PV=nRT\)が成立している(空気を理想気体として扱う, 補遺1)
  2. BB弾とバレルの間に空気の漏れは無い(補遺2)

ということです。
この条件下では、バレル内部の運動方程式は以下のように導くことができます。

ピストン-BB弾間の体積\(V\)はピストンの位置とBB弾の位置の2つに依存します。
なので、ピストンの位置による時間変化を\(x_0(t)\)、BB弾の位置を\(x(t)\)、バレルの断面積を\(S\)とすると
\(
V(t)=S\cdot (x(t)-x_0(t))
\)

と記述されます。
気体の状態方程式より、
\(
PV=nRT
\)

であり、右辺を定数\(c_0\)と置き、圧力\(P(t)\)をBB弾に掛かる力F(t)、すなわち\(P(t)=\frac{F(t)}{S}\)とすれば、
\(
\begin{align}
\frac{F(t)}{S}\cdot S\cdot (x(t)-x_0(t))=c_0 \\
F(t)=\frac{c_0}{x(t)-x_0(t)}
\end{align}
\)
と書けます。よって、空気抵抗と共に書けば、運動方程式は
\(
\displaystyle m\frac{d^2}{dt^2}x(t)=\left(\frac{c_0}{x(t)-x_0(t)}-P_0\cdot S\right)-6\pi \eta R \cdot v(t)-\frac{1}{2}C_d \rho \pi R^2 \cdot v^2(t)
\)

と導けます。ここで、\(P_0\)は大気圧、\(\eta\)は粘性率、\(R\)はBB弾の半径、\(C_d\)は抗力係数、\(\rho\)は空気の密度を意味します。
空気抵抗に関する詳細は球体の空気抵抗と係数をご覧ください。

あとはこの運動方程式を解いて、BB弾が出てくるバレル位置において適当な速度になるような物理的な意味を十分に持つ\(c_0, x_0(t)\)を入れて計算すればいいのです。

ここで1つ言えることは、外部と内部の圧力が一定になる最適なバレル長というものが存在する、ということです。
このバレル長を使えば、BB弾が射出されたときに、銃口で気圧差は無く、それによって生じる不確定な要素を最小にすることができるでしょう。

そのバレル長とは、式(4)の右辺が0になる時であり、計算によって求めることができます。
いくつになるかは実際求めてみないと想像できません。数cmであるのか数十mであるのか、はたまた30cmとかその位なのか・・・。

ちなみに実銃ではほぼ延々に火薬によって非常に多くのガスが生成され、加速されるためこのようなことは起こらないでしょう(あるかもしれませんが、とてつもなく長いバレルでなければ実現しないと思われます)。

補遺


補遺1
理想気体として扱える条件は、
低い気圧(分子の数が少なく、衝突等が無視できる)かつ高い温度(分子間力が分子の運動エネルギーに比べて無視できる)
です。
どうやら実在気体では10気圧以下、室温以上でこの条件は良く満たされ、理想気体とのずれは1%以内のようです[1]。
[1]を引用すると、

一般に,沸点の低い酸素・窒素・水素・ヘリウム等は,室温またはそれ以上の温度で10atm以下の圧力の場合,理想気体の値の1%以内で理想気体に近い性質を示す。

とあります。ピストンで空気が圧縮されたとき、BB弾とピストン間の圧力が10気圧以上にならなければ良い近似だと言えるでしょう。

補遺2
BB弾とバレル内側との間における空気の漏れはどの程度あるのかはわかりません。

バレルの内径:6.08mm (東京マルイのものはこれらしいですが、信頼できるページが見当たりませんでした。)
BB弾の直径:5.95mm [2]

完全に理想気体として、時間を考慮しないで扱うのであれば、少しでも空間が空いていればどれだけ早く押し込もうとも気圧は一定になろうとするため、BB弾に掛かる力はゼロです。
なので時間を考慮した方程式が必要でしょう。
連続の式
\(
\displaystyle \frac{\partial \rho}{\partial t}+\nabla \cdot \vec{j}=0
\)

を用いれば時間と共に漏れ出る空気量を知ることが出来そうです。
ただ、まずは簡単のため、空気が漏れない場合をまず考えましょう。

簡単に空気の漏れを表現するには、あたかもピストン-BB弾間の体積\(V\)が時間と共に増加していく、と考えることである程度の効果を入れることができると考えられます。

実際に式を解いてから、空気の漏れが起こる場合と起こらない場合でどれほどの差が出てくるのか?を考えることでこの効果を入れるべきか入れないべきかを考えるのがいいでしょう。

良く空気の漏れがーと言われているので本当かどうか理論的に検証するべきでしょう。

参考文献


[1] 実在気体 -第2節 気体の状態方程式
[2]PERFECT HIT -TOKYO MARUI

『集弾性アップへの道』 BB弾とバレル内部に隙間があることが写真で確認できます。

非線形Schrödinger方程式のソリトン解

非線形シュレディンガー方程式
\(
\displaystyle i\frac{\partial \Psi}{\partial t}=-\frac{1}{2}\frac{\partial^2 \Psi}{\partial x^2}-g|\Psi|^2\Psi
\)

にはある解析解が存在します。それがソリトン(soliton)解と呼ばれるもので,上式のソリトン解は
\(
\displaystyle \Psi(x,t)=\sqrt{\Omega}\ {\rm sech}\left\{\sqrt{\Omega}\left(x\sqrt{g}-\frac{gV}{2}t\right)\right\}\cdot \exp\left\{i\frac{V\sqrt{g}}{2}x-i\frac{g}{2}\left(\frac{V^2}{4}-\Omega\right)t\right\}
\)

です。(\(g>0\),\(\Omega\):ソリトンの振幅、\(V\):ソリトンの速度に関するパラメータ、ソリトン自体の速度は\(V\sqrt{g}/2\))


スポンサーリンク


ソリトンの歴史的背景


「非線形」とは重ね合わせの原理が成り立たない系です。

1844年、スコットランドのJ.Scott-Russellによって孤立した波(solitary wave)を観測した事が報告されました J.Scott-Russellによる報告”Report on Waves”(リンク先のSR44.pdf, 16.3MB))。
当時の認識では、波は波動方程式で記述され、その波の速度\(v\)は\(v=f\lambda\)の元、一定である。だからパルス状の波は異なる波長の波の重ね合わせで書けているはずで、時間と共に分散していくはず。なのになぜ時間が経過しても孤立した波が存在できるのか?という事で大きな論争となりました。

60年後の1895年、オランダのKortewegとde Vriesによって”浅い水の波”を記述する非線形偏微分方程式(KdV方程式)が提出され、この方程式の特解として孤立波が与えられました。
孤立波は、

  1. 空間的に局在した波が、その性質(速さや形)を変えずに伝搬する
  2. 孤立波は互いの衝突に対して安定であり、各々の個別性を保持する

という性質を持つ非線形波動と定義されます[1]。
2番目の、粒子のような性質を持つことから、solitary に接頭語-on をつけ、soliton(ソリトン)と名づけられました。

その後、1981年に佐藤幹夫がソリトンの統一理論(佐藤理論やKP理論)を発表しました。
これによりソリトン方程式(ソリトンを記述し,かつ厳密に解ける方程式)に決着が付きました。
ソリトン方程式は非線形なのに厳密に解ける、可積分系である。

ソリトン方程式を解く方法は([4]を引用しますが)

上で指摘したように,logistic方程式が解けるからくりとソリトン方程式が解けるからくりはよく似ています.違いは,logistic方程式が変数変換一発で線形常微分方程式になってしまったのに対し,ソリトン方程式の場合は変数変換で双線形形式になり,双線形形式の解として行列式が現れ,行列式の中身に簡単な線形方程式が現れるというところです.しかし,離散化で保存するべき構造は明らかです.まず,解の中身の線形方程式を離散化し,行列式の構造をそのまま使って双線形形式を作る.最後に変数変換して非線形のレベルに戻ればよい.

となるそうです。

また、ソリトン方程式の特徴である、無限個の対称性(無限個の保存量)は、Gardner変換という変換をすることで証明できるそうです[5]。
これ以上はこの分野の専門家ではないので話せません。

ちなみに津波もソリトンの一つとみなせます。

ソリトン解が生まれるイメージ


なぜソリトン解が生まれるのでしょうか。
今、孤立した波(空間的に凸)を考えます。この時、

エネルギー的に安定になろうとして密度を均一にするために広がろうとする効果

粒子間を結び付ける引力相互作用(例えば水面だったら水と水との分子間力等)のため集まろうとする効果

のつり合いによって、丁度均衡が保たれるとき、このソリトン解が生まれます。

・・・実は、ソリトン解には2種類あります。
それは明るい(Bright)ソリトン解と暗い(Dark)ソリトン解です。
今まで話していたのは全て明るいソリトン解です。
暗いソリトン解とはどういったものでしょう。
暗いソリトン解とは、ある部分が空間的に凹んでいる、孤立した解です。

エネルギー的に安定になろうとして密度を均一にするためにその凹みを埋めようとする効果

粒子間の斥力相互作用のために粒子間を避けようとする効果

のつり合いによって、丁度均衡が保たれるとき、この暗いソリトン解が生まれます。
暗いソリトン解が生まれるのは斥力相互作用の時で、斥力相互作用を持つ系というのは、調べた限りでは量子力学のボーズ・アインシュタイン凝縮体で、暗いソリトンは渦ソリトンという形で現れるそうです。これ以上の具体例は分かりませんでした。もしも具体例を知っているという方は教えていただければ幸いです。
暗いソリトンの解析解は参考文献[1]の本に紹介されているので、それをご参考にしてください。

非線形シュレディンガー方程式におけるソリトン解


では本題の、非線形シュレディンガー方程式における(明るい)ソリトン解を考えましょう。
下の形の非線形シュレディンガー方程式を考えます。
\(
\displaystyle i\frac{\partial \Psi}{\partial t}=-\frac{1}{2}\frac{\partial^2 \Psi}{\partial x^2}-g|\Psi|^2\Psi
\)

ここで、\(\Psi=\Psi(x,t)\)で、\(g\)はの値で相互作用の強さ(この場合、引力相互作用)を表します。

この非線形シュレディンガー方程式のソリトン解\(\Psi(x,t)\)は、
\(
\displaystyle \Psi(x,t)=\sqrt{\Omega}\ {\rm sech}\left\{\sqrt{\Omega}\left(x\sqrt{g}-\frac{gV}{2}t\right)\right\}\cdot \exp\left\{i\frac{V\sqrt{g}}{2}x-i\frac{g}{2}\left(\frac{V^2}{4}-\Omega\right)t\right\}
\)

であり、\(\Omega\)はソリトンの振幅の大きさ、\(V\)はソリトンの速度を決めるパラメータを表します。ソリトン自体の速度は\(V\sqrt{g}/2\)となります([2]を参考)。
また、\({\rm sech}(x)\)は双曲線関数の一種(双曲線正割と呼ばれる)であり、
\(
\displaystyle {\rm sech}(x)=\frac{1}{\rm cosh(x)}=\frac{2}{e^{x}+e^{-x}}
\)
を表します。

解析解のプロット

解析解をプロットします。fortranコードは下のほうに載せておきます。
\(g=2, V=1, \Omega=1\)とすると、以下の振る舞いが観測されます。
ここで紫はソリトン解の実部、緑は虚部、青は絶対値2乗を表します。
動画は1枚当たり原子単位系で0.1秒、合計で10秒間のシミュレーションです。
また、このソリトンの速度は\(V\sqrt{g}/2\sim 0.7071\)です。
soliton1

スポンサーリンク

fortranコード


非線形シュレディンガー方程式のソリトン解(解析解)を出力します。
gnuplot上で下のスクリプトを実行してください。
(ただし、gnuplot ver4.6以降に限る)

omega=1e0
V=1e0
g=2e0
x0=-5e0 # initial position

set xr[-10:10]
set yr[-1.5:1.5]
set samples 1000
set xl "x[a.u.]"

sech(x)=2e0/(exp(x)+exp(-x))
amp(x,t)=sqrt(omega)*sech(sqrt(omega)*((x-x0)*sqrt(g)-g*V*t*0.5e0))
phase(x,t)=V*sqrt(g)*x*0.5e0-g*t*0.5e0*(V*V*0.25e0-omega)
soliton(x,t)=amp(x,t)*exp({0e0,1e0}*phase(x,t))


#set term gif animate delay 10 optimize size 960,720
#set output 'movie.gif'
do for[i=0:100:2]{
   t=i*0.1e0
   plot abs(soliton(x,t))**2 lw 3 lc 1 lt 1 ti sprintf("|\psi|^2, t=%2.1f",t),\
        real(soliton(x,t)) lw 3 lc 2 lt 2 ti "Real",\
        imag(soliton(x,t)) lw 3 lc 3 lt 3 ti "imag"
}
#set out
#set terminal wxt enhanced

もっとソリトンについて知りたい方はまず参考文献[3]を読むことをお勧めします。
その後、[4]を読み、[5]を読み、[1]の本を読むのが良いと思われます。
[3],[4]は簡単な表現を用いてソリトンとその後の発展について記述されています。

参考文献


[1]和達三樹著 『非線形波動 (現代物理学叢書) 』岩波書店 (2000年) p.7
[2]和達三樹著 『非線形波動 (現代物理学叢書) 』岩波書店 (2000年) p.29

[3]ソリトンの数学 – Researchmap
[4]ソリトン ~ 不思議な波が運んできた,古くて新しい数学の物語 ~
[5]〔連載〕非線形波動―ソリトンを中心として―第5章 逆散乱法


↑画像のフォントはキユマヤ園様による数式フォント -びゅんびゅん→SSSです!


球体の空気抵抗と係数

まとめ


半径\(R\)の完全な球体に働く空気抵抗力\(F_d\)は
\(
\displaystyle F_d = 6\pi R \eta v + \frac{1}{2} C_d \rho \pi R^2 v^2
\)

です。方向も含めるのであれば
\(
\displaystyle \vec{F}_d = -6\pi R \eta |\vec{v}| \frac{\vec{v}}{|\vec{v}|} – \frac{1}{2} C_d \rho \pi R^2 |\vec{v}|^2 \frac{\vec{v}}{|\vec{v}|}
\)

となります。

ここで、

  • \(\eta\ \mathrm{[kg \cdot m^{-1}\cdot s^{-1}]}\) : 空気の粘性率(Viscosity) ※動粘度ではないです。
  • \(\rho\ \mathrm{[kg \cdot m^{-3}]}\) : 空気の密度
  • \(v\ \mathrm{[m\cdot s^{-1}]}\) : 球体の速度の大きさ
  • \(R\ \mathrm{[m]}\) : 球の半径

を表します。
\(C_d\)は抗力係数(Drag Coefficient)と呼ばれ、厳密にはレイノルズ数\(R_e\)の関数\(C_d=C_d(R_e)\)となります。
レイノルズ数\(R_e\)とは無次元量で、
\(\displaystyle R_e=\frac{v \rho L}{\eta}\)
という量です。ここで\(L\)は、物体の大きさを表し、半径\(R\)の球の場合は直径\(L=2R\)に相当する量です。

球の場合の抗力係数\(C_d(R_e)\)をフィッティングしたものは、論文[2]より、
\(
\displaystyle \small C_d=\frac{24}{R_e}+\frac{2.6\left(\frac{R_e}{5.0}\right)}{1+\left(\frac{R_e}{5.0}\right)^{1.52}}
+ \frac{0.411\left(\frac{R_e}{263000}\right)^{-7.94}}{1+\left(\frac{R_e}{263000}\right)^{-8.00}} + \left(\frac{R_e^{0.80}}{461000}\right)
\)

と表されます。
Cdグラフ

20℃で1気圧(101325[Pa])、湿度0%の空気中では、
\(
\begin{align}
\eta &= 18.2\times 10^{-6}\mathrm{[kg \cdot m^{-1}\cdot s^{-1}]} \\
\rho &= 1.205\mathrm{[kg \cdot m^{-3}]}
\end{align}
\)

という値になります。

温度、圧力、湿度を考慮する場合

温度\(T[\mathrm{K}]\)、圧力\(P[\mathrm{Pa}]\)の時,
空気の粘性率\(\eta\mathrm{[kg \cdot m^{-1}\cdot s^{-1}]}\)は、
\(\displaystyle \eta\mathrm{[kg \cdot m^{-1}\cdot s^{-1}]}=1.487 \times 10^{-6}\cdot \left(\frac{T\mathrm{[K]}^{1.5}}{T\mathrm{[K]}+117}\right)\)
であり、
空気の密度\(\rho\ \mathrm{[kg \cdot m^{-3}]}\)は、
\(\displaystyle \rho\ \mathrm{[kg \cdot m^{-3}]}=0.0034856447\cdot \frac{P\mathrm{[Pa]}}{T\mathrm{[K]}-0.670}\)
です。

この2つの式を[1]が正しいと思って確かめた温度の範囲-10℃~40℃で大体3桁くらい一致します。
この式は下に説明してある空気の場合を代入し、求めたものです。

2016/03/13 追)
湿度がある場合、空気の密度\(\rho_{w}\mathrm{[kg \cdot m^{-3}]}\)は湿度0%の時と比べて減少します。
温度\(t\mathrm{[{}^{\circ}C]}\)、湿度M(湿度60%の時、\(M=0.60\))、気圧\(P\mathrm{[Pa]}\)の場合、空気の密度\(\rho_{w}\mathrm{[kg \cdot m^{-3}]}\)は
\(
\displaystyle \rho_{w}\mathrm{[kg \cdot m^{-3}]}=
\left(0.0034856447\cdot\frac{P\mathrm{[Pa]}}{t\mathrm{[{}^{\circ}C]}+272.48}\right)\cdot
\left(1-M\cdot 0.378\cdot\frac{100\cdot 6.1078\times 10^{\frac{7.5\cdot t\mathrm{[{}^{\circ}C]}}{t\mathrm{[{}^{\circ}C]}+237.3}}}{P\mathrm{[Pa]}}\right)
\)

と記述されます。
ここで用いた関係式は
[5]に圧力pの水蒸気を含んだ空気の密度の関係式

[6]、[7]の飽和水蒸気圧を求めるTetensの式
を用いています。また、atmによる単位系では[8]に記述があります。


スポンサーリンク



スポンサーリンク

風の抵抗力の詳細


これと同じことは弾道計算(BB弾)の理論に記述してあります。
その中の風の抵抗力を部分を抜き出したものになります。

風(流体)の抵抗力の大きさ\(F_d\)は一般には次元解析により
\(
\displaystyle F_d=\frac{\eta^2}{\rho}\sum_{n=1}^{\infty}K_n \left(\frac{v \rho L}{\eta}\right)^n
\)

として与えられます[3]。ここで

  • \(\eta\ \mathrm{[kg \cdot m^{-1}\cdot s^{-1}]}\) : 流体の粘性率(Viscosity) ※動粘度ではないです。
  • \(\rho\ \mathrm{[kg \cdot m^{-3}]}\) : 流体の密度
  • \(v\ \mathrm{[m\cdot s^{-1}]}\) : 物体の速度の大きさ
  • \(L\ \mathrm{[m]}\) : 物体の大きさ(半径\(R\)の球の場合、\(L\)は直径\(L=2R\)に相当する量)
  • \(\frac{v \rho L}{\eta}\) : レイノルズ数(\(R_e\),Reynolds number)と呼ばれる無次元量

です。次元解析から具体的な定数の値\(K_n\)を求めることはできません。

速度の3乗以上に比例する項というのは見かけません。
おそらく、物体の動きを記述するために必要なパラメータがこの2つを考慮すれば含まれるため、\(n=3\)以上の項は冗長になるからでしょう。つまり、第一項には空気の密度\(\rho\)が含まれておらず、第二項まで入れることで\(\eta,\rho\)があらわに含まれるようになるからだと思います。

\(n=3\)以上の高次の効果は、抗力係数\(C_d\)に押し込めていると思います。

上の式の\(n=1,n=2\)をそのまま書き下すと
\(F_d=K_1\eta L v +K_2\rho L^2 v^2\)
となりますが、慣例として同じ次元を持つように式を変形させて
\(F_d=C_1 L\eta v +C_d S \left(\frac{1}{2}\rho v^2\right)\)
の形で記述されます。ここでSは速度に垂直な、物体の断面積です。

\(C_1\)は物体が完全な球で、小さく、その速度がゆっくりである場合(レイノルズ数が2以下くらい)、流体の運動を記述する方程式、
ストークス方程式(Stokes’ law、ナビエ・ストークス方程式のレイノルズ数が小さい極限の方程式)
から導くことができて\(C_1=3\pi\)であることが理論的に導けます。

\(C_d\)は抗力係数(Drag Coefficient)と呼ばれ、
これはレイノルズ数の関数となり、
\(C_d=C_d(R_e)\)
となります。
この抗力係数は非常に厄介らしいです。物体の形状に依存したりするらしく、球の場合でも、どうやら実験的に求められているようです。
Wikipediaの抗力係数のページに若干の記述があります。
球の場合の抗力係数の関数\(C_d(R_e)\)をフィッティングした関数は論文[2]より、
\(
\displaystyle \small C_d=\frac{24}{R_e}+\frac{2.6\left(\frac{R_e}{5.0}\right)}{1+\left(\frac{R_e}{5.0}\right)^{1.52}}
+ \frac{0.411\left(\frac{R_e}{263000}\right)^{-7.94}}{1+\left(\frac{R_e}{263000}\right)^{-8.00}} + \left(\frac{R_e^{0.80}}{461000}\right)
\)

として書けるそうです。実際にプロットしてみるとこんな感じです。
Cdグラフ
確かに抗力係数のグラフと同じようになります。

よって半径Rの球の場合、風の抵抗力\(F_d\)の大きさは結局、
\(
\begin{align}
F_d &= C_1 L\eta v +C_d S \left(\frac{1}{2}\rho v^2\right) \\
&= 3\pi \cdot 2R\cdot \eta v +\frac{1}{2} C_d \rho \pi R^2 v^2 \\
&= 6\pi R \eta v + \frac{1}{2} C_d \rho \pi R^2 v^2
\end{align}
\)
として与えられます。

各定数の値
名称 捕捉
空気の粘性率(Viscosity)\(\eta\) \(\small 18.2\times 10^{-6} \\ \small \mathrm{[kg \cdot m^{-1}\cdot s^{-1}]}\) 空気中20度において。また、粘性率は数十気圧から数気圧の間圧力には依存しません。
が、温度に依存し、ある温度\(T_1[K]\)での粘性率\(\eta_1\)が分かっているならば、
温度\(T_2[K]\)での粘性率\(\eta_2\)はサザランドの定数Cを用いて、
\(
\eta_2=\eta_1 \left(\frac{T_1[K]+C}{T_2[K]+C}\right)\left(\frac{T_2[K]}{T_1[K]}\right)^{3/2}
\)
と記述されます。空気の場合、C=117です[4]。
乾燥した空気の密度\(\rho\) \(\small 1.205\mathrm{[kg \cdot m^{-3}]}\) 乾燥した空気の密度は温度と圧力に依存します。温度\(t[℃]\)、圧力をH[torr]とすると
\(
\rho=\frac{1.293}{1+0.00367t[℃]}\cdot \frac{H[torr]}{760}
\)
です[5]。

参考文献


[1]水・空気の物性 密度 粘度 動粘度
[2]Faith A. Morrison, “Data Correlation for Drag Coefficient for Sphere,” Department of Chemical Engineering, Michigan Technological University, Houghton, MI
[3]伊東敏雄著『な~るほど!の力学』学術図書出版社 (1994) p.220
[4]国立天文台編『理科年表 平成21年 第82冊』丸善株式会社 (2009) p.378
[5]国立天文台編『理科年表 平成21年 第82冊』丸善株式会社 (2009) p.376
[6]相対湿度の月別平年値 -理科年表オフィシャルサイト
[7]飽和水蒸気量 -wikipedia
[8] 空気 -wikipedia

ちなみに、
\(1 [\mathrm{torr}] = 133.322368 [\mathrm{Pa}]\\ 1[\mathrm{Pa}]=1[\mathrm{N\cdot m^{-2}}]=1[\mathrm{kg \cdot m^{-1}\cdot s^{-2}}]\)
という関係です。

弾道計算(BB弾)の結果2、違う重さでの比較

本稿の主要な結果は、様々なパラメータでの、BB弾の軌道の詳細なデータです。

この結果に付随して、
0.25g0.8Jの軌道は
0.20g1.2Jの軌道と同じであることが分かりました。

上下振れ幅が最小になる軌道における、重さとエネルギー、ゼロイン位置を様々にとった時のそれぞれのデータです。
弾道計算を数値計算によって行い、結果を考察することが本稿の目的となります。
※本稿では規定のエネルギーを超える場合のデータがありますが、このデータは全て数値シミュレーションであるため、こういったエアガンを持っている訳ではありません。

弾道計算に関するその他ページ
弾道計算(BB弾)の理論
BB弾の回転量について(実験との比較)
弾道計算(BB弾)の結果
弾道計算の結果2, 比較と詳細データ←今ここ
弾道計算(BB弾)のコード(fortran90)
バレル内部でのBB弾の方程式
水中下でのBB弾の弾道計算


様々なパラメータの最適な軌道


ここで載せるデータは、以下の組み合わせです。

BB弾の重さ: 0.20g , 0.25g
エネルギー:0.60J 0.65J 0.70J 0.75J 0.80J 0.85J 0.90J
ゼロイン位置:25m 30m 35m 40m 45m 50m 55m 60m(0.25gのみ)
各々でどのような軌道を描くかの計算結果を載せます。
破線はホップ無しの軌道に対応しています。

↓これは低画質です。高画質版はこのリンク(4.3MB)を踏んでください。
弾道詳細データ_高画質2 - コピー_c


最適な軌道に合わせるためには下向きに初速を与えます。
その時、下向きに何度傾けて撃てばいいか、のデータはこちらです。

ゼロインと下向き角度


前提


まず、規定速度うんぬんは別にしまして、事実を述べます。

  1. BB弾は重いほど弾はまっすぐ飛ぶ
  2. 射出速度が早いほど弾はまっすぐ飛ぶ
  3. ホップを掛けるほど弾は上下に揺れる

これらは紛れもない事実です。
よって、BB弾を飛ばす最高の条件とは、
「射出速度を早くし、ホップは余りかけないで重いBB弾を使うこと」
となります。
だからこそ、射出速度を速めようとして規定速度の話になります。

最適な軌道とは?


BB弾を飛ばすうえで最適な軌道とはどんな軌道でしょう?
それはホップによる上下方向の揺れを最小限に抑える軌道です。
この最適な軌道とは、BB弾の重さとエネルギーとゼロイン位置を決めた時の理想的な軌道、ということです。

ゼロイン位置と上下の振れ幅とは何かは、下の画像をご覧ください。
0.25,0.8_compressed

重さとエネルギーを決めても、どの角度で射出すればいいか、回転数はいくつか、など他のパラメータが残ります。
それを決めるため、ホップによる上下方向の揺れを最小限に抑える軌道を定めます。

この軌道を ”最適な軌道” と呼ぶことにします。

重さの違いによる軌道の具体的な影響


今度は、ゼロイン位置を固定します。
そして、BB弾の重さとエネルギーを変化させたとき、どのような軌道をたどるか見てみましょう。
zeroin50_2_cc

この結果から分かることは2点あります。

1つ目は、エネルギーを上げていくと軌道の変化が小さくなる
2つ目は、重いBB弾の軌道ほどまっすぐ飛ぶ

ということです。

エネルギーを上げると軌道の変化が小さくなる

0.20gの場合でゼロインを50mに合わせ、エネルギーを0.1Jずつ増やしていったとき、上下振れ幅がどのくらい減少していくのか見てみましょう。

0.60J→0.70J … 14.2cm
0.70J→0.80J … 10.4cm
0.80J→0.90J … 8.0cm
0.90J→1.00J … 6.2cm
1.00J→1.10J … 5.0cm
1.10J→1.20J … 4.2cm
1.20J→1.30J … 3.6cm
1.30J→1.40J … 3.0cm

となります。エネルギーを変えるよりも重さを変えるべきです。
それだけで軌道は大きく変わります。

例えば、0.80Jで射出できるエアガンで、使うBB弾を0.25gで射出すると、この時、軌道は0.20g1.2Jの軌道にほぼ一致します。
m020J12m025J08_2_c

まぁ、軌道自体は一致するのですが、残念ながら到達時間は一致しません。
約0.1秒(まばたき程度)、0.20g1.2Jのほうが早くなります。
zeroin50_2_time_c

追記)
0.25gで0.90Jは、0.20gの1.35Jの軌道におおよそ相当します。
0.30gで0.80Jは、0.20gの1.60J, 0.25gの1.05Jの軌道におおよそ相当します。
0.30gで0.90Jは、0.20gの1.80J, 0.25gの1.20Jの軌道におおよそ相当します。

0.25gで0.8Jでうまく合わせられたら50m飛ばしても上下の振れ幅40cmです。
上に20cm、下に20cmしかずれません。
また、エネルギーを上げても空気抵抗の強さは速度の2乗に比例して強くなるため、それに見合うようなの飛距離の伸びは得られません。


補足


「流速チューン」というものがあるそうですね。
流速チューンの問題点と改善案について -週休5日

流速チューン比較テスト その1「初速変化」 -Gunsmithバトン
にあるように。

本稿ではこれらのことについては全く触れていません。
というのも、銃口から射出された直後の初速と回転数の2つだけが軌道に影響するためであり、発射するまでの過程に何があったか?などどうでもいいのです。
数値計算の観点から言いますと、通常の軌道よりも「まっすぐ遠くまで飛ぶ」場合、初速が上がっていること以外には考えられません
BB弾を飛んでいる最中で加速させる要因など無いのです。


※1
なぜBB弾を使う限り初速を稼ぐことが無駄なのか?
これは空気抵抗が関係します。
空気抵抗の大きさは弾道計算(BB弾)の理論で書いたように、BB弾の半径\(R\)と速度\(v\)にのみ依存します。
空気抵抗力を半径\(R\)と速度\(v\)の関数として\(F_d(R,v)\)、質量\(m\)とすると、その運動方程式は
\(
m\frac{d^2\vec{r}}{dt^2}=F_d(R,v)
\)
であり、両辺をmで割れば
\(
\frac{d^2\vec{r}}{dt^2}=F_d(R,v)/m
\)
となります。BB弾の半径や、温度、空気は変えられないので、可変なパラメータは質量と速度のみです。
そして右辺だけに注目すれば、質量が大きいほど空気抵抗力があたかも小さくなるのです。
すなわち、空気抵抗力を減らそうとするにはBB弾の重さを変えるほかない。
また、初速を変えても空気抵抗は速度の2乗で効いてくるため、皮肉なことに速度が早いほど空気抵抗力が強くなっていくのです。
よって、高いエネルギーでは初速を上げてもそれに見合っただけの結果は得ることはできないのです。