sikino のすべての投稿

sedコマンドによるファイルの書き換え

シェルスクリプトによってファイルを書き換えます。

  1. sedコマンドの概要
  2. 具体例一覧
    1. 特定の数値を変えたい場合
    2. 特殊文字(/等)を含む文字列の変更
    3. 置換する数値を変数での指定
    4. ある数値以降を全て変更したい
    5. ファイルの空白行の削除
    6. ファイルの行頭に数値を入れる
    7. ファイルの最後に改行を入れる
    8. 繰り返し変更したい場合

    9. do文で何回も変更したい
    10. 全く違う変数を繰り返し代入したい
    11. 整数の変数を規則的に何回も変更したい
    12. 小数点を含む変数を規則的に何回も変更したい
    13. あるディレクトリ以下にある全ファイルを変数にしたい

sedコマンドの概要


対象とするファイル名”input.ini”の中身に

&input
 a=7,
 b=13,
 c='./dat/',
&end

という文が書いてあるとします。
これのファイルの中身を端末からのコマンドで書き換えるためにはsedコマンドを利用します。

sedコマンドとは、置換用のコマンドです。
例えば、コマンドライン上で

sed -e "s/a=7/a=9/" ./input.ini

と入力すると”input.ini”の中身の”a=7″に一致する場所を”a=9″に置換して出力します。
ここで重要なのが、”a=7″に一致する場所だけなので、後ろにあるコンマは影響されません
実際に実行してみると

$ sed -e "s/a=7/a=9/" ./input.ini
&input
 a=9,
 b=13,
 c='./dat/',
&end

となってコマンドライン上に出力されます。
この時、書き換わった内容はコマンドラインのみに出力されるため、元のファイル自体に書き換えはされません

元のファイルの書き換えを行いたければ、出力先を適当なファイル、例えば”tmp”というファイルにし、それを”input.ini”にリネームすればいいです。
すなわち、

sed -e "s/a=7/a=9/" ./input.ini > tmp
mv tmp input.ini

とすればいいのです。実際に実行して中身を見ますと

$ sed -e "s/a=7/a=9/" ./input.ini > tmp
$ mv tmp input.ini
$ cat input.ini
&input
 a=9,
 b=13,
 c='./dat/',
&end

となります。

sedコマンドの概要はここで終了です。
以降は具体例とシェルスクリプトのコード(ファイル名”q.sh”),実行例を載せていきます。

特定の数値を変えたい場合


対象とするファイル”input.ini”

&input
 a=7,
 b=13,
 c='./dat/',
&end

“q.sh”の中身

#!/bin/sh
sed -e "s/a=7/a=9/" ./input.ini > tmp
mv tmp input.ini

実行例

$ sh q.sh
$ cat input.ini
&input
 a=9,
 b=13,
 c='./dat/',
&end

特殊文字(/等)を含む文字列の変更


sedコマンドを使うために/を使います。ならば/を置換したい場合どうすればいいでしょう。
この時はsedコマンドで使っている区切りを表す/を別の特殊文字(例えば%)に置き換えればよいです。

対象とするファイル”input.ini”

&input
 a=7,
 b=13,
 c='./dat/',
&end

“q.sh”の中身

#!/bin/sh
sed -e "s%c='./dat/'%c='./data/'%" input.ini > tmp
mv tmp input.ini

実行例

$ sh q.sh
$ cat input.ini
&input
 a=7,
 b=13,
 c='./data/',
&end

置換する数値の変数での指定


対象とするファイル”input.ini”

&input
 a=7,
 b=13,
 c='./dat/',
&end

“q.sh”の中身

#!/bin/sh
pr=5
sed -e "s/a=7/a=${pr}/" ./input.ini > tmp
mv tmp input.ini

実行例

$ sh q.sh
$ cat input.ini
&input
 a=5,
 b=13,
 c='./dat/',
&ends

ある数値以降を全て変更したい


コンマも含めて消したいとします。
対象とするファイル”input.ini”

&input
 a=7,
 b=13,
 c='./dat/',
&end

“q.sh”の中身

#!/bin/sh
sed -e "s/c=.*/c=12345/" ./input.ini > tmp
mv tmp input.ini

実行例

$ sh q.sh
$ cat input.ini
&input
 a=7,
 b=13,
 c=12345
&end

ファイルの空白行の削除


ファイルに含まれる空白行を全て削除します。
対象とするファイル”input.ini”

&input
 a=7,

 b=13,

 c='./dat/',


&end

“q.sh”の中身

#!/bin/sh
sed -e '/^ *$/d' input.ini > tmp
mv tmp input.ini

実行例

$ sh q.sh
$ cat input.ini
&input
 a=7,
 b=13,
 c='./dat/',
&end

ファイルの行頭に数値を入れる


対象とするファイル”input.ini”

&input
 a=7,
 b=13,
 c='./dat/',
&end

“q.sh”の中身

#!/bin/sh
i=1
while [ ${i} -le 1000 ]
do
    sed -e "${i}s/^/${i} /" input.ini > temp
    mv temp input.ini

    i=`expr $i + 1`
done

実行例

$ sh q.sh
$ cat input.ini
1 &input
2  a=7,
3  b=13,
4  c='./dat/',
5 &end
6

“q.sh”中、1000をファイルに応じて変更してください。

※また、行番号を入れるだけの場合、catコマンドを利用して、

cat -b input.ini > tmp
mv tmp input.ini

とするのが簡単だと思います。
オプション”-b”は空行を含めず行番号を付けるという意味です。空行を含めたい場合は”-n”でできます。

ファイルの最後に改行を入れる


対象とするファイル”input.ini”

&input
 a=7,
 b=13,
 c='./dat/',
&end

“q.sh”の中身

#!/bin/sh
sed '$a\\n' input.ini > tmp
mv tmp input.ini

実行例

$ sh q.sh
$ cat input.ini
&input
 a=7,
 b=13,
 c='./dat/',
&end
   

$

※どうしても2行分改行されてしまうようです。
sedコマンドで1行分の改行は、どうすればよいか分かりませんでした。

苦肉の策ですが、適当なファイル(例えばnn.txt)に改行だけを入力し、catコマンドでファイルを結合すると一応できます。
—nn.txt—

————
を作っておいて、
#!/bin/sh
cat input.ini nn.txt > tmp
mv tmp input.ini
とすれば”input.ini”の後ろに”nn.txt”が結合されるため、1行分の改行が可能です。

do文で何回も変更したい


対象とするファイル”input.ini”

&input
 a=7,
 b=13,
 c='./dat/',
&end

“q.sh”の中身

#!/bin/sh
for pr in 0 1 2
do
    sed -e "s/a=.*/a=${pr},/" input.ini > tmp
    mv tmp input.ini
   
    cat input.ini
done

実行例

$ sh q.sh
&input
 a=0,
 b=13,
 c='./dat/',
&end
&input
 a=1,
 b=13,
 c='./dat/',
&end
&input
 a=2,
 b=13,
 c='./dat/',
&end

全く違う変数を繰り返し代入したい


対象とするファイル”input.ini”

&input
 a=7,
 b=13,
 c='./dat/',
&end

“q.sh”の中身

#!/bin/sh
i=1
for j in 3d0 939d0
do
    echo "count -> ${i}"
    sed -e "s/a=.*/a=${j},/" ./input.ini > tmp
    mv tmp input.ini
   
    grep "a=" input.ini
   
    i=`expr $i + 1`
done

実行例

$ sh q.sh
count -> 1
 a=3d0,
count -> 2
 a=939d0,
$

シェルスクリプトの中の”grep”は”a=”に一致する行を書き出させるものです。

整数の変数を規則的に何回も変更したい


対象とするファイル”input.ini”

&input
 a=7,
 b=13,
 c='./dat/',
&end

“q.sh”の中身

#!/bin/sh
i=1
while [ ${i} -le 3 ]
do
    echo "count -> ${i}"
    sed -e "s/a=.*/a=${i},/" ./input.ini > tmp
    mv tmp input.ini
   
    grep "a=" input.ini
   
    i=`expr $i + 1`
done

実行例

$ sh q.sh
count -> 1
 a=1,
count -> 2
 a=2,
count -> 3
 a=3,
$

小数点を含む変数を規則的に何回も変更したい


小数点を含む演算は”expr”コマンドではできません。
シェルスクリプトで数値計算を行う”bc”コマンドを用いましょう。
ここでは、変数\(v\)に、
\(v=i\times0.7,\;\;(i=1,2,3)\)
を代入する演算を考えましょう。

対象とするファイル”input.ini”

&input
 a=3,
 b=13,
 c='./dat/',
&end

“q.sh”の中身

#!/bin/sh
i=1
while [ ${i} -le 3 ]
do
    echo "count -> ${i}"

    v=$(echo "scale=3; $i*0.70 " | bc | sed -e 's/^\./0./g')
    echo "$v"
   
    sed -e "s/a=.*/a=${v},/" ./input.ini > tmp
    mv tmp input.ini
   
    grep "a=" input.ini
   
    i=`expr $i + 1`
done

実行例

$ sh q.sh
count -> 1
0.70
 a=0.70,
count -> 2
1.40
 a=1.40,
count -> 3
2.10
 a=2.10,
$

あるディレクトリ以下にある全ファイル名を変数にしたい


状況

.
├── a.sh
├── dat
├── data
│   ├── asdfgh.d
│   ├── qwerty.d
│   └── zxcvbn.d
└── input

inputの中身

&input
  filename="./data/xxxxx.d",
  key = 0,
&end

a.shの中身

#!/bin/sh
i=1
for j in ./data/*
do
    sed -e "s%filename=.*%filename="${j}",%" input > tmp
    mv tmp ./dat/${j##*/}_ex
    echo ${j##*/}
    grep "filename=" ./dat/${j##*/}_ex
done

———————-

a.sh実行

$sh a.sh
asdfgh.d
  filename="./data/asdfgh.d",
qwerty.d
  filename="./data/qwerty.d",
zxcvbn.d
  filename="./data/zxcvbn.d",

a.sh実行後の構成

.
├── a.sh
├── dat
│   ├── asdfgh.d_ex
│   ├── qwerty.d_ex
│   └── zxcvbn.d_ex
├── data
│   ├── asdfgh.d
│   ├── qwerty.d
│   └── zxcvbn.d
└── input

例えばdat/asdfgh.d_exの中身

&input
  filename="./data/asdfgh.d",
  key = 0,
&end

参考先


フィルタを使用した文字列操作 1 -UNIX & Linux コマンド・シェルスクリプト リファレンス
sed の使い方 -サーバエンジニアの知恵袋
制御構文 -シェルスクリプト入門
How to pass results of bc to a variable – Ask Ubuntu
bcによる少数の演算 -SWEng TIPs

水素原子の解析解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)
※参考した本はこれです。

openMPによる並列計算

まとめ


fortran90では、

program main
  implicit none
  integer::i,omp_get_thread_num

  ! number of threads for parallel computation
  call omp_set_num_threads(2)

  !$omp parallel do
  do i=1,100
     write(6,*) omp_get_thread_num()
  enddo
  !$omp end parallel do  

  stop
end program main

とすれば2個のスレッドで並列計算され、コンパイル時に
intel® fortranでは、

ifort -openmp main.f90

gfortranでは、

gfortran -fopenmp main.f90

とすれば良いです。


fortran90でのお話です。

openMPは複数のCPUを使って並列計算をするインターフェースのことです
C/C++,fortranで対応しています。
何も指定しない場合、計算は基本的に1つのコアのみを使います。
Linuxを使っている場合、恐らくデフォルトで使えるようになっているはずです。

並列計算を行うまでの大まかな手順は2つ。

  1. プログラム内部に文章を書く
  2. コンパイル時にオプションを追加
    • intel®fortranコンパイラの場合
    • gfortranコンパイラの場合

です。

※今使っているパソコンの並列計算できるCPUの数の確認は、

cat /proc/cpuinfo | grep processor

です。proc/cpuinfoの場所は使っているOSによって違う可能性があります。
”linux mint cpu 確認”
等とgoogle等で検索するのが手っ取り早いでしょう。

プログラム内部に文章を書く


プログラム内でopenMPを使用するdoループを挟むように、

!$omp parallel do
do i=1,100
    call test(i)   
enddo
!$omp end parallel do

と書きます。基本はこれだけです。
並列計算に使えるスレッドが2個である場合、並列計算というのは

スレッド0

do i=1,50
    call test(i)   
enddo

スレッド1

do i=51,100
    call test(i)   
enddo

と割り振る操作なのです。

上記のようにdoループを分割するため、注意しなければならない点として、
doループに用いる変数に対して別々に計算しても計算結果が変わらない事
が挙げられます。
上記例では、i=0から始めて、i=1,2,3,・・・と順次計算する必要は無く、i=5,10,1,3,・・・みたいな乱雑な順番でも計算結果が変わらない状況で使わなければなりません。
特に、i=2を計算したい時にi=1の計算結果を使う場合は並列化することは出来ません。

実際には同期したりとかで、出来なくはないですが面倒くさく、僕はやったことが無いので書きません。

もしも並列計算したいスレッド数を制限したい場合、

!$omp parallel do

!$omp parallel do num_threads(2)

のように指定するか(この例ではその部分の並列計算を2つのスレッドで行う)、
並列計算が始まる前のプログラムの初めの方に

call omp_set_num_threads(2)

と書きましょう。これを書くと全ての並列計算箇所が2個のスレッドで行われることになります。

また、関数omp_get_thread_num()はスレッドにつけられた番号を返す事が出来ます。引数は要りません。
なので、これを利用すると何個のスレッドが使われているのか直に確認できます。

補足ですが、MKL(マス・カーネル・ライブラリ)のルーチン内で使われる並列計算のスレッド数は

call mkl_set_num_threads(2)

と書くことで制御できます。

コンパイル時のオプション


intel® fortranでは、

ifort -openmp main.f90

gfortranでは、

gfortran -fopenmp main.f90

としてコンパイル、実行しましょう。

実例


例として

program main
  implicit none
  integer::i,omp_get_thread_num

  ! number of threads for parallel computation
  call omp_set_num_threads(2)

  !$omp parallel do
  do i=1,4
     write(6,*) omp_get_thread_num()
  enddo
  !$omp end parallel do

  stop
end program main

を動かすと同出力されるか見てみましょう。
gfortranコンパイラで、使えるスレッドの数は2つです。
上記のコードを

gfortran -fopenmp main.f90

としてコンパイルし、実行すると

./a.out
           1
           0
           0
           1

という結果が得られます。
ここに出ている数字はスレッドの番号で、0,1と番号付けされたスレッドが確かに出力されている事が分かります。1,0,0,1の順番はどうでも良いです。実行するごとに変わるので、これは同時並行処理されているためでしょう。

call omp_set_num_threads(1)

として計算を行うと、

./a.out
           0
           0
           0
           0

と表示されることから、確かに使うスレッド数が制御されていることが分かります。

doループの並列処理を処理が終わった順に行う場合


doループ内で、あるときだけ計算が重くなる場合があります。
この場合、均等にdoループの変数を分けてしまうとCPU1つあたりの処理が大きく異なってしまいます。そこで、doループの変数を終わった順に処理していく指定方法があります。それがschedule構文です。以下のようにすると、i=1から計算が終わった順にiを増やしていきます。

program main
  implicit none
  integer::i,omp_get_thread_num,N

  ! number of threads for parallel computation
  call omp_set_num_threads(2)

  N=100
  !$omp parallel do schedule(dynamic,1)
  do i=1,N
     write(6,*) omp_get_thread_num()
  enddo
  !$omp end parallel do
 
  stop
end program main

とします。ただし、1回あたり計算時間がほとんどかからない場合、オーバーヘッド時間が無視できなくなるため、使うときは吟味してから使用するべきです。

また、並列化しても共通の変数を使うことができます。以下のようにすると、計算の進み具合がわかります。

program main
  implicit none
  integer::i,omp_get_thread_num,N,prog

  ! number of threads for parallel computation
  call omp_set_num_threads(2)

  N=100
  prog=0
  !$omp parallel do shared(prog) schedule(dynamic,1)
  do i=1,N
     write(6,*) omp_get_thread_num()
     prog=prog+1
     if(mod(prog,50).eq.0)write(6,'(A,i6,A,i6)')" ",prog," / ",N
  enddo
  !$omp end parallel do
 
  stop
end program main

とします。

経験的に、openmpを使用する際はただ一つのサブルーチンを用いて、並列個所を簡単にするのが良いです。

  !$omp parallel do
  do i=1,N
     a(i)=0
     do j=i,N
         a(i)=a(i)+1
     enddo
     a(i)=a(i)*i     
  enddo
  !$omp end parallel do

と、!$~!$につらつらと書くのではなく、

  !$omp parallel do
  do i=1,N
    call test(i,N,a(i))
  enddo
  !$omp end parallel do

として、サブルーチン1つにするだけの方がバグが出にくいです。

openMPの機能はこれだけではありません。
調べることをお勧めします。ここではこれ以上書きません。

参考文献


インテル® コンパイラーOpenMP*入門

正多角形とスピログラフの数式

正多角形とスピログラフの数式です。

  1. 正多角形(2次元)
  2. 正多角形(3次元)
  3. 星型多角形(2次元)
  4. 星型多角形(3次元)
  5. トロコイド(2次元)
  6. トロコイド(3次元)
  7. スピログラフ

2016/10/24
幾何学模様を作りたい方は幾何学模様(gnuplot)へどうぞ。

[adsense1]

正多角形の式(2次元)


正n角形の数式は[1]を参考に修正しました。
正n角形の数式は、
\(
\begin{equation}
\left\{
\begin{aligned}
x&\displaystyle =\cos t \cdot \frac{\cos \frac{a_n}{2}}{\cos\left[a_n\left\{\frac{t}{a_n}-{\rm floor}\left(\frac{t}{a_n}\right)\right\}-\frac{a_n}{2}\right]}\\
y&\displaystyle =\sin t \cdot \frac{\cos \frac{a_n}{2}}{\cos\left[a_n\left\{\frac{t}{a_n}-{\rm floor}\left(\frac{t}{a_n}\right)\right\}-\frac{a_n}{2}\right]}
\end{aligned}
\right.
\end{equation}
\)
です。ここで、\(a_n=\frac{2\pi}{n}\)と置きました。また、\({\rm floor}(x)\)は床関数(床関数と天井関数)を表します。
媒介変数\(t\)の範囲は(\(0\le t \lt 2\pi\))です。
グラフを書かせてみますとこうなります。

polygon3_10_c
n=3,4,5,6
n=7,8,9,10の順で、単位円を一緒に記述しています。
線の色の濃さは媒介変数の値に応じて変化させています。

gnuplotのスクリプトはこちら

多面体の式(3次元)


半径1の多面体は、平面のx,yの媒介変数による正n角形の関数をそれぞれpc(t),ps(t)と置くと、
\(
\begin{equation}
\left\{
\begin{aligned}
x& = {\rm ps}(u)\cdot {\rm pc}(v) \\
y& = {\rm ps}(u)\cdot {\rm ps}(v) \\
z& = {\rm pc}(u)
\end{aligned}
\right.
\end{equation}
\)
と書けます。
媒介変数\(u,v\)の範囲は(\(0\le u \lt \pi, \;\; 0\le v \lt 2\pi\))です。

プロットすると以下のようになります。
polygon3d_3_10_c
n=3,4,5,6
n=7,8,9,10の順です。

gnuplotスクリプトはこちら

星型多角形(2次元)


\(
\begin{equation}
\left\{
\begin{aligned}
x&\displaystyle =\cos t \cdot \frac{\cos a_n}{\cos\left[a_n -\frac{1}{n}\cos^{-1}\{\cos(nt)\}\right]}\\
y&\displaystyle =\sin t \cdot \frac{\cos a_n}{\cos\left[a_n -\frac{1}{n}\cos^{-1}\{\cos(nt)\}\right]}
\end{aligned}
\right.
\end{equation}
\)
です(\(n\ge 5\))。ここで、\(a_n=\frac{2\pi}{n}\)と置きました。
グラフを書かせてみるとこのようになります。
starpolygonp_c

gnuplotスクリプトはこちら。
※\(\cos^{-1}(\cos(x))\)は不安定です。
場合により書けなくなることがあります。ここでは上記関数の代わりに、正弦波ではない周期関数の数式
で記述した関数\(s(a,x)\)を用いています。

星型多角形(3次元)


3次元の場合は多面体の時と同じように、平面のx,yの媒介変数による星型多角形の関数をそれぞれstpc(t),stps(t)と置くと、
\(
\begin{equation}
\left\{
\begin{aligned}
x& = {\rm stps}(u)\cdot {\rm stpc}(v) \\
y& = {\rm stps}(u)\cdot {\rm stps}(v) \\
z& = {\rm stpc}(u)
\end{aligned}
\right.
\end{equation}
\)
と書けます。
媒介変数\(u,v\)の範囲は(\(0\le u \lt \pi, \;\; 0\le v \lt 2\pi\))です。
グラフはこのようになります。
starpolygon3dp_c

gnuplotスクリプトはこちら。

トロコイド(2次元)


トロコイドと呼ばれる曲線があります[2]。
これをgnuplotで書いてみましょう。数式は[2]を参考にして以下の通りです。
\(
\begin{equation}
\left\{
\begin{aligned}
x& = \left[(n-1)\cos t +h\cos((n-1)t)\right]/n \\
y& = \left[(n-1)\sin t +h\sin((n-1)t)\right]/n
\end{aligned}
\right.
\end{equation}
\)
\(h=1\)の時、半径が1の円に内接するように決めています。
また、\(n\)が整数の時、内接する円が外側の円をちょうど一周回った時に元の位置に戻る条件です。
別に整数でなくても構いません。その時はもっと複雑なものになることでしょう。

書いてみますと曲線はこうなります。
自由度\(h\)を1にしてグラフを書いてみましょう。
trocoid2d_3_10_c

自由度hを少し変えてみましょう。
trocoid2dh_3_10_c

トロコイド(3次元)


3次元の場合は多面体の時と同じように、平面のx,yの媒介変数によるトロコイドの関数をそれぞれhcyc(t),hcys(t)と置くと、
\(
\begin{equation}
\left\{
\begin{aligned}
x& = {\rm hcys}(u)\cdot {\rm hcyc}(v) \\
y& = {\rm hcys}(u)\cdot {\rm hcys}(v) \\
z& = {\rm hcyc}(u)
\end{aligned}
\right.
\end{equation}
\)
と書けます。
媒介変数\(u,v\)の範囲は(\(0\le u \lt \pi, \;\; 0\le v \lt 2\pi\))です。

trocoid3d_3_10_c
n=3,4,5,6
n=7,8,9,10の順です。

gnuplotスクリプトはこちら

[adsense2]

スピログラフ(2次元)


スピログラフは遊んだ方も多いと思います。これです。
1024px-Spirograph
©Alexei Kouprianov/2007/CC-BY 2.5
スピログラフはトロコイドの特別な時です。スピログラフをもっと一般化したのがトロコイドです[2]。

図示すれば、スピログラフと言うのは
speq_compressed
という二つの円からなる軌道です。黄色の点線の長さが等しい、という束縛条件が課せられます。
数式は上のトロコイドと同じで
\(
\begin{equation}
\left\{
\begin{aligned}
x& = \left[(n-1)\cos t +h\cos((n-1)t)\right]/n \\
y& = \left[(n-1)\sin t +h\sin((n-1)t)\right]/n
\end{aligned}
\right.
\end{equation}
\)
と書けます。

もう少し考えを一般化して、周りを取り囲んでいる形状が円ではない時を考えましょう。その時の数式を考えます。
考え方はトロコイドと一緒ですが、もう少し拡張しましょう。

今、2つの媒介変数による軌道を持っていたとします。1つの媒介変数の軌道上を動かせばスピログラフの考えと同じになります。図で示すと以下のようになります。
spirograph_any_c
スピログラフは両方円であり、滑りが無い(cを決定する条件)場合に対応するわけです。ただ、円である必要性はないので、多角形であろうが、何だろうが数式上では書けるわけです。
興味深い事に、大きな軌道に沿って小さな軌道を動かした、と見ようが小さな軌道に沿って大きな軌道を動かしたとみても数式上では一緒なのです。面白いですね。

いくつか例を挙げます。
五角形上を動かす場合はこうです。
parameter1

数式として,
\(
\begin{equation}
\left\{
\begin{aligned}
x& = f_{1x}(t)+r\cdot f_{2x}(s\cdot t) \\
y& = f_{1y}(t)+r\cdot f_{2y}(s\cdot t)
\end{aligned}
\right.
\end{equation}
\)
と書き表すことにします。これを元にgnuplotのスクリプトを書けば以下のようになります。

色々に変化させた時の色々なときのグラフです。
spr2_any_compressed

rを綺麗な数字にしなければもっと複雑な図が得られます。
また、この考え方は3つの媒介変数の軌道、4つの媒介変数の軌道…として容易に拡張できます。
もっと複雑な模様が欲しければ、もう少し増やしてみましょう。

例えば3つの場合はこうなります。
parameter3_spi
このスクリプトを置いておきます。

この考えを拡張していくと”フーリエ級数の図示化”と言うのができるわけです。

実際にやっている方がいるので紹介しておきますね。

参考文献


[1]Is there an equation to describe regular polygons? -MATHEMATICS
[2]Hypotrochoid -Wolfram Mathworld

[3]Is there an equation to describe regular polygons? -Mathematics
2021/11/29追記)滑らかに星と円を遷移する場合、このページが丁寧です。

アニメーションに便利な関数

アニメーションに便利な関数、手法です。

点の表示


plot sprintf("< echo ''%f %f''", x0,y0)

とすると点\((x_0,y_0)\)が描写されます。

窓関数


窓関数です。関数\(W(x)\)は、
\(
\displaystyle W(x)=\frac{r^{2^n}}{(x-x_c)^{2^n}+r^{2^n}}
\)
書かれます。この関数は\(r\)が窓の横幅を示し、\(n\)が窓の端の様子を決めるパラメータです。窓の広がりはおよそ\(2r\)です。

fwi(n,r,i,ic)=(real(r)**(2**n))/((real(i-ic))**(2**n)+real(r)**(2**n))

加速、減速


extanh
\(
\begin{align}
y(x)&=C_0\tanh(a(x-x_c))+C_1 \\
& C_0 = \frac{f_0-f_1}{t_0-t_1} \\
& C_1 = \frac{-t_1f_0+t_0f_1}{t_0-t_1} \\
& t_0 = \tanh(a(x_0-x_c)) \\
& t_1 = \tanh(a(x_1-x_c)) \\
& x_c = \frac{x_0+x_1}{2}
\end{align}
\)
の関数は点\((x_0,y_0)\)から点\((x_1,y_1)\)に\(\tanh(ax)\)で与えられる関数で書きます。
関数の中心は\((x_c,y_c)=((x_0+x_1)/2,(y_0+y_1)/2)\)になります。
画像の2つの赤い点は\((x_0,y_0)\)と\((x_1,y_1)\)を表しています。
\(a\)は変化の度合いを調節するパラメータです。

この関数は非常に使えます。ラミエルを作っている最中にこの関数を考えましたが、想像以上に役立ちます。

fth(a,x,x0,f0,x1,f1)=(f0-f1)*tanh(a*(real(x)-(x0+x1)*0.5e0))/(tanh(a*(real(x0)-(x0+x1)*0.5e0))-tanh(a*(real(x1)-(x0+x1)*0.5e0)))+(-f0*tanh(a*(real(x1)-(x0+x1)*0.5e0))+f1*tanh(a*(real(x0)-(x0+x1)*0.5e0)))/(tanh(a*(real(x0)-(x0+x1)*0.5e0))-tanh(a*(real(x1)-(x0+x1)*0.5e0)))
 # x0,f0 --> x1,f1

この関数を用いた応用です。
extime_tanh

これはある時間内で高さを反転させるものです。コードはこちら。

set xr[-6:6]
set yr[0.8:2.2]
set size ratio 0.5 1
   
a=1e0
at=0.4e0
x0=-5e0
f0=1e0
x1=5e0
f1=2e0

#set term gif enhanced animate optimize delay 8 size 800,500
#set output 'extime_tanh.gif'

i0=0
i1=20
do for[i=i0:i1]{
    ft0=fth(at,i,i0,f0,i1,f1)
    ft1=fth(at,i,i0,f1,i1,f0)
    plot fth(a,x,x0,ft0,x1,ft1) w l lt 2 lc 3 lw 5 ti sprintf("a=%3.2f, at=%3.2f, (x_0,f_0)=(%3.2f,%3.2f), (x_1,f_1)=(%3.2f,%3.2f)",a,at,x0,f0,x1,f1), sprintf("< echo ''%f %f''", x0,ft0) pt 7 ps 2 lc 7, sprintf("< echo ''%f %f''", x1,ft1) pt 7 ps 2 lc 7

  }

#set out
#set terminal wxt enhanced

Wolfram alphaの便利な使い方

オンラインで使える数学用の検索エンジン
Wolfram alpha
を使うといろいろな事が出来ます。それらを紹介します。

積分


integral from 0 to infty sin(x)/x

基本中の基本、積分です。複雑な積分でも定積分でも不定積分でも無限大でも受け付けてくれます。
ただし、計算してくれるかは分かりません。
例えば\(\rm{sinc}\)関数\(\left(\rm{sinc}(x)=\frac{sin(x)}{x}\right)\)のゼロから無限大までの積分は、
20160202-222431_c
のように出力してくれます。\(\pi\)で出力されるのは凄いです。

\(\rm{sinc}\)関数のように有名な関数であれば、

integral sin(x)/x

としたほうがたくさんの情報が得られます。

微分方程式の解


y''(x)+y(x) = sin(x)

と入力すると微分方程式を解いてくれます。また、

y''(x)+y(x) = sin(x), y(0)=1, y'(0)=0

な風にすると、初期条件を与えたうえで解いてくれます。
20160202-221952_c
自分の与えた数式通りに認識しているか、確かめましょう。

行列の対角化


Eigenvectors[{-2,2,4},{-2,4,2},{-2,1,4}]

行列の対角化を行ってくれます。上を入力すると行列
\(
\begin{eqnarray}
\left(
\begin{array}{ccc}
-2 & 2 & 4 \\
-2 & 4 & 2 \\
-2 & 1 & 4
\end{array}
\right)
\end{eqnarray}
\)を対角化し、その固有値と固有ベクトルを出力してくれます。

級数展開


関数の級数展開を出力してくれます。
例えば、

Series[BesselJ(1,x),{x,0,5}]

と入力するとベッセル関数\(J_1(x)\)を\(x\)について、\(x=0\)回りで、\(x^5\)まで展開して出力してくれます。
20160202-223943_c

また、

Series[cos(x),{x,1,10}]

とすると、\(\cos(x)\)を\(x\)について、\(x=1\)回りで、\(x^{10}\)まで展開して出力してくれます。

また、

Series[exp(ix)/x,{x,0,4}]

としても展開してくれます。この場合はローラン(Laurent)展開であり、特異点周りの展開となっています。
上の場合、出力として、
20160202-221030_c
という結果が得られます。
もう一つ、ガンマ関数\(\Gamma(x)\)の\(x=0\)まわりでのローラン展開です。

Series[gamma(x),{x,0,3}]

と入力するとガンマ関数\(\Gamma(x)\)を\(x\)について、\(x=0\)回りで、\(x^3\)まで展開して出力してくれます。
20160202-223636_c
すごい・・・

キャラクターの曲線


様々なキャラクターの曲線のグラフを出力してくれます。
その曲線はWolfram alpha named parametric curvesにまとめられていましたので紹介します。
もしくは、↓をクリックすることで展開されます。

上の他には、wolframで

fictional character curves

animal curves

person curves

と入力して,moreを押していってください。

例を挙げます。

tachikoma like curve

“tachikoma”は、攻殻機動隊に出てくる戦車です。こんな感じです。
tachikoma_wolfram_c

——————–

charmander like curve

“charmander”は、ヒトカゲの事のようです。
charmander_wolfram_c
———————
もしも、上のグラフの数式をコピペしたい場合は、各グラフの下に
plaintext_wolfram
という欄の、”plain text”というところをクリックすれば数式を取得できます。

vmware playerの容量を減らす

VMware playerの容量が大きくなってきた時に大きく効果のあったことです。

ここでは、
ホストOS:windows8.1
上に
vmware playerにインストールされてるゲストOS:Linux Mmint 17.1 rebecca
がある、という状況で行った手順をまとめます。

ゲストOS(VMwareを立ち上げたlinux内)のコマンドラインで以下のことを行ってください。VMware toolsを使います。インストールされてない方はLinuxMint導入時にやることに書いたので、VMware toolsのインストールを行ってください。

sudo vmware-toolbox-cmd disk list
sudo vmware-toolbox-cmd disk shrink /

で後は1~2時間待てばいいです。
僕はこれによって、windows上のフォルダ”Virtual Machines”の占める容量を43GB→22GBにまで減らせました。


  1. 仮想マシンのシャットダウンと起動(サスペンドではなく、シャットダウンです。)
  2. ファイル整理
    (特にwindowsから仮想OSへドラッグ&ドロップする際に生成される

    ~/.cache/vmware/drag_and_drop/

    の削除)

  3.     sudo vmware-toolbox-cmd disk list
        sudo vmware-toolbox-cmd disk shrink /

    の実

参考文献


VMWareのゲストイメージ圧縮方法

三角波、のこぎり波、矩形波、その他の数式

正弦波ではない周期関数、三角波、矩形波等々の数式による表現です。
if文は使っていません。
床関数\(floor(x)\)を主に用いています。

これらが唯一の作り方ではなく、最善の作り方でもありません。

非連続になる関数の場合の特定な関数の場合、中間の値を取るようにしています(フーリエ級数の収束より)。


階段関数

——-
f1
z_f1_x_c
——-
f2
z_f2_x_c
——-
fl
z_fl_x_c


のこぎり波

s1
z_s1_x_c
s2
z_s2_x_c
saw
z_saw_x_c


矩形波

k1,k2,k,ka,kb,kk
z_k1_abx_c

z_k2_abx_c

z_k_abx_c

z_ka_abx_c

z_kb_abx_c

z_kk_abx_c


その他

sl,s,spk
tr,se,de,mo
well1,well2

sl
z_sl_abx_c

s
z_s_ax_c
三角波は、\({\rm acos}(\cos(x))\)とも書けます。

spk
z_spk_abx_c

tr
z_tr_abx_c

se
z_se_abx_c

de
z_de_abx_c

mo
z_mo_abx_c

well1
zwell1_abx_c

well2
z_well2_abx_c

well

追記
\(\displaystyle
f(a,b,x) = sgn(abs(2*(x-b)/a)-1e0)
\)
でokです。aは井戸の幅、bは井戸の中心を示します。


まとめ
periodic_functions_c

全てをまとめたgnuplotのコードは以下のものです。

set size ratio -1
set yr[-2:2]
set xr[-5:5]
set samples 101
set grid
set xtics 1
set ytics 1
set mxtics 2
set mytics 2

f1(x)=floor(x)
f2(x)=-floor(-x)-1e0
fl(x)=0.5e0*(floor(x)-floor(-x)-1e0)

s1(x)=x-floor(x)
s2(x)=x+floor(-x)+1e0
saw(x)=x-0.5e0*(floor(x)-floor(-x)-1e0)

k1(a,b,x)=f1((x+a)/(a+b))-f1(x/(a+b))
k2(a,b,x)=f2((x+a)/(a+b))-f2(x/(a+b))
k(a,b,x)=fl((x+a)/(a+b))-fl(x/(a+b))
kk(a,b,x)=2e0*(k(a,b,x)-0.5e0)

ka(a,b,x)=0.5e0*(sgn(k(a,b,x)-0.25e0))+0.5e0
kb(a,b,x)=0.5e0*(sgn(k(a,b,x)-0.75e0))+0.5e0

sl(a,b,x)=kk(a+b,a+b,x+0.5e0*b+a)*k(a,b,x)

s(a,x)=2e0*abs(s2(x/(a*2))-0.5e0)*kk(a*2e0,a*2e0,x+a)
tr(a,b,x)=k(a,b,x)*s(0.5e0*(a+b),x-0.5e0*b)*(1e0+b/a)+sl(b,a,x-b)

spk(a,b,x)=-abs(kk(a,b,x))+1e0

se(a,b,x)=(x-b-(a+b)*fl(x/(a+b)))/a*k(a,b,x)+0.5e0*(k(a,b,x)*(1e0-b/a)+1e0)*spk(a+b,a+b,x)

de(a,b,x)=-se(a,b,x)*2e0*(k2(a+b,a+b,x-0.5e0*b)-0.5e0)
mo(a,b,x)=de(a*0.5e0,b+a*0.5e0,x+a*0.5e0)+de(a*0.5e0,b+a*0.5e0,-(x+a*0.5e0))

well1(a,b,x)=0.5e0*(sgn((abs(2e0*(x-b)/a)-1e0)+ka(1e0,1e0,abs(2e0*(x-b)/a))-0.5e0))+0.5e0
well2(a,b,x)=0.5e0*(sgn((abs(2e0*(x-b)/a)-1e0)+kb(1e0,1e0,abs(2e0*(x-b)/a))-0.5e0))+0.5e0

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