sikino のすべての投稿

水素原子の原点に電子を見出すか?

問題


水素原子の電子の基底状態において、原点に電子を見出すことができるでしょうか?

問題設定


非相対論の範囲で、水素原子の電子を考えます。
電子の満たすシュレーディンガー方程式は、プランク定数\(\hbar\),電子の電荷\(e\),電子の質量\(m\)
をそれぞれ\(\hbar=1, e=1, m=1\)とする原子単位系で

と書けます。ここで、\(\nabla\)はデカルト座標系\((x,y,z)\)における微分演算子\(\nabla=\left(\frac{\partial }{\partial x}, \frac{\partial }{\partial y}, \frac{\partial }{\partial z} \right)\)であり、\(\psi(\mathbf{r})\)は波動関数を表します。
規格化は

で行います。特に、これから中心力に対する問題を考えていくので、球面座標系で変数分離を考えます。球面座標系において動径方向\(r(=\sqrt{x^2+y^2+z^2})\)と角度方向\(\theta(=\cos^{-1}(z/r)), \varphi((=\tan^{-1}(y/x)))\)の\((r,\theta, \varphi)\)座標における波動関数は

と書きます。ここで規格化は

とします。過程は飛ばしますが、シュレーディンガー方程式の固有値\(E=-1/2\)に属する基底状態の波動関数は、量子数\((n, l, m)=(1,0,0)\)を持つ状態として指定され、\(R(r)\)は

と書けます。\(\theta, \varphi\)方向の関数は\(Y_{0,0}=({4\pi})^{-1/2}\)ですので、波動関数は

となります。

さて、存在確率密度はどう考えればよいでしょう?
まず規格化から、

と計算されるわけですから、原点からの距離\(r\sim r+dr\)の範囲に電子を見出す存在確率密度\(f(r)dr\)の係数\(f(r)\)は

と書けます。

一方、存在確率密度は波動関数の絶対値二乗ですので、

ですから、位置\((x,y,z)\sim (x+\Delta x,y+\Delta y,z+\Delta z)\)の範囲に見出す確率\(g(x,y,z)dxdydz\)の係数\(g(x,y,z)\)は

と書けます。
さて、ここで疑問が生じます。原点\((x,y,z)=(0,0,0)\)において、電子を見出す可能性はあるのでしょうか?

つまり、動径方向の関数については原点における\(f(0)\)は式(8)より

であり、\(g(0,0,0)\)は式(10)より

と書けるのでどちらが正しいのか?本当に観測を行ったとき、原点で電子を見出すことはあるのか?

という問題です。つまり、グラフにすると以下のようになります。\(g(x,y,z)\)については、\(x=y=0\)としてグラフにしています。

これで注目するのは、原点における値がゼロか有限か?で大きく異なっている点です。

座標系が違うだけで物理が変わることはありませんので、これは解釈の問題です。
現実には原点には原子核があるので原点における電子を観測することはできないので、頭の中で考えます。
結論としては、どちらも正しくて解釈が異なる。そして、原点で電子を見出すことがあっても良いです。
その理由を述べていきます。

考察


まず、数値実験を行い事実を確認します。
デカルト座標系\((x,y,z)\)における波動関数\(\psi(\mathbf{r})=\frac{1}{\sqrt{\pi}}e^{-\sqrt{x^2+y^2+z^2}}\)は紛れもない事実ですので、フォン・ノイマンの棄却法に従って、\((x,y,z)\)空間で一様な変数を作り出し、この確率密度に従った乱数を表示させます。すると以下のような図を得ます。

例えば\(x=0, y=0\)に固定してz軸上の波動関数は、

となりますので、\(z=0\)においてゼロではない有限の値をとるため確率密度が存在しそうなので、原点に電子を見出してもよさそうです。

続いて動径方向の分布を考えましょう。棄却法によって採用された\(n\)番目の点\((x_n,y_n,z_n)\)(図に表示されている点)の原点からの距離\(r_n=\sqrt{x_n^2+y_n^2+z_n^2}\)を調べてみます。そして、\(r\sim r+\Delta r\)の範囲にある点数を数え、観測された個数と原点からの距離の関数を考えます。これは動径分布関数の\(f(r)\)に等しいはずです。すると、以下のような図を得ます。

確かに、確率密度の係数にかかる分布\(f(r)\)になっていることが分かります(最大値を1にするように規格化しています)。

さて、以上の数値実験からやはりどちらも数式通りであり、正しいことが分かりました。
すると解釈の問題でしょう。どのように解釈していけばよいかを考えていきます。

まず動径方向の確率密度\(f(r)dr=r^2R^2(r)dr\)ですが、これは正確には

と書いたほうが良いでしょう。意味は、半径が\(r\sim r+dr\)の間を占める球殻の体積は\(4\pi r^2 dr\)であり、その体積の中に粒子の見出しやすさ\(\frac{R^2(r)}{4\pi}\)が掛かっている、という意味です。

半径が\(r\to 0\)の極限において、式を解釈すると、

球殻の体積がゼロのとき、その中に粒子を見出す確率はいくつか?

という問いになります。見出すことができる体積はゼロなので、体積ゼロの領域に電子を見出すことはできません。なので、

となります。すなわち、\(r^2 dr\)の部分がゼロになるだけで、\(R^2(r)\)の部分は値を持っていても良い、ということになります。少し詳しく言えば、球殻の体積がゼロの時に確率があると確率が無限大になってしまうので、少なくとも\(r\to 0\)で\(R(r)\propto r^{n},~(n>-1)\)であることが課されます。

以上の話が正しければ、原点では体積がゼロなので、電子を見出す確率\(r^2R^2(r)\)は\(r=0\)でゼロ、すなわち電子を見出さない、という結論になりそうです。

これまでの考察を行っても矛盾が生じたままで、何一つ解決していません。
「原点において電子を見出す体積がゼロだから、電子を見出せない」という結論が誤りなのか、もう少し疑ってみます。

体積がゼロでも電子を見出すことがあり得ることを仮定した場合、これはどういう意味を持つでしょうか。

2つの解釈を説明します。

  • 有限の大きさによって、本当の原点に電子を見出すことはない
  • 例えば将来的に超高精度な観測機器ができた場合、原点に電子を見出す場合にしても、現実ではどう頑張ってもプランク長(10^{-35}\mathrm{m})以上の分解能は無いはずです(電子の大きさもありますしね。しかもプランク長以下の長さは物理的な意味がないと言われています。一応、プランク長よりも短くなる長さがあるようですが、無限小よりかは大きいでしょう)。なので、原点に電子を見出した!と思っても、実は原点からほんのちょっとずれているのかもしれません。
    だから、本当の原”点”に電子は見出さないのかもしれません。

  • 点で電子を見出すか、密度を考えるかで区別しなければならない
  • 点自体に体積は無いので、確率密度がゼロだからと言って原点に電子を見出しても良い

二番目の点自体に体積は無いので、確率密度がゼロだからと言って原点に電子を見出しても良い、が正しいと思いますが…すみません。調べたりしたのですが、なかなか目当ての問題や解答がなく、探すのを諦めました。
存在確率密度なので、(体積中の電子を見出す確率)/(見出す体積)ですが、原点であればこの分母がゼロです。
この場合に、原点に電子を見出したからと言って、それは体積でも何でもない0次元の情報ですから、確率「密度」の情報に対応させることができません。なので見出しても問題がないと考えます。

一方、原点以外で波動関数の節となる確率密度がゼロになる点においては体積が有限であり、その点では本当に見出すことは本当の意味でありません。

結論


点自体に体積は無いので、確率密度がゼロだからと言って原点に電子を見出しても良い。

…で合っていると思いますが、明確な参考文献を見つけられませんでした。
メモとして、これまでの考察を書きました。

補足)カスプについて


一つ、本当の原点に電子を見出してしまった場合、ポテンシャルエネルギーが発散しないのか心配になります。
すなわち、ポテンシャルエネルギー項が\(-\frac{1}{r}\)であり、無限大になるので非物理的な気がします。

だからと言って、ポテンシャルエネルギーが発散するから非物理的である、という結論は間違いです。
この原点における発散はカスプ(Cusp)と呼ばれ、運動エネルギーの項と相殺するため、消えてしまう、もしくは消えなければなりません(カスプ条件)。カスプについて説明するために、もう一度シュレーディンガー方程式に立ち返って考えてみます。

左辺のハミルトニアンのポテンシャル項\(-\frac{1}{r}\)は\(r=0\)で発散しています。しかし束縛状態では、波動関数の二乗は電子の存在確率密度ですので、全空間で積分したら有限にならなければなりません。つまり右辺\(E\psi(\mathbf{r})\)は発散してはなりません。
以上から、\(\psi(\mathbf{r})\)は連続でなければならず、\(r=0\)で

を計算したら\(\frac{1}{r}\psi(\mathbf{r})\)の項が出てきてポテンシャルエネルギーの項を打ち消さなければなりません。
波動関数に課されるこの条件は、カスプ条件と呼ばれます。

陰関数の等高線を数値計算で求める

まとめ


関数

を求める問題を考えます。
\(f(x_0, y_0)=0\)を満たす\((x_0, y_0)\)が与えられている場合、\((x_1\equiv x_0+\Delta x, y_1\equiv y_0+\Delta y)\)を求めることを考えていきます。
まず、

の変換を行い、初期値\((x_0, y_0)\)を\((x’_0, y’_0)\)に変換します。
ここで、\(\theta\)は定数で

です。\((x’_0, y’_0)\)を新たな初期値として微分方程式

を解き、\(x_0′\)から微小距離\(h\)だけ発展させた\(x_1’\equiv x_0′ + h\)とその値\(y_1’\equiv y_0′ + \Delta y\)を求めます。
その後、逆変換

を経ることで\(f(x_1,y_1)=0\)を満たす\((x_1,y_1)\)の値を得ることができます。

導出


陰関数

を満たす\((x,y)\)を求めます。ある点\((x_0, y_0)\)において\(f(x_0, y_0)=0\)を満たすことが分かっている場合、
陰関数の全微分より、\((x_0, y_0)\)の近傍で

が満たされています。\(f(x_0,y_0)=0\)であり、\(\Delta x\approx \Delta y \to 0\)ならば、\(O(\Delta^2)\)の項を無視して

が成立する変化のみが許されます。両辺を\(\Delta x\)で割って、\(\Delta x\to 0\)の極限をとれば、

を得ます。すなわち、式(9)を何らかの手段、例えば数値計算で解けば等高線を得ることができます。

発散への対処


さて数値計算によって微分方程式(9)を解いて等高線を描く場合、\(y\)が\(x\)の一価関数になることはほとんどありません。
そのため、式(9)の右辺の分母がゼロになることがあり、微分が発散するため計算を行うことができなくなることがあります。言い換えれば、発散しない範囲しか計算できなくなるということです。
これを回避するためには様々な方法があると思いますが、例えば変数変換を適宜行いながら積分を実行していくことでこの問題を回避することができます。

発散を回避するためのアイデアは以下の通りです。
\(y\)を\(x\)の関数と見たとき、微分が発散するということは、分母がゼロになることを意味しています。
しかし、その点において90°回転させて\(y\)を\(x\)の関数としてみれば微分はゼロに等しいため、問題なく計算を進めることができます。
図示すれば以下の通りになります。

左図は座標を変化せずにそのままの関数ですが、右図のようにある点を計算したいときにその点近くの傾きと等しいように回転した座標から見るようにします。こうすることで発散を回避することができます。

これから変数変換について考えていきます。新たな変数\(x’, y’\)を用いて\(x=x(x’, y’),~y=y(x’, y’)\)(またはその逆関数)と書けているならば、

と書いても良いです。すなわち式(9)を解く代わりに\(x’\)と\(y’\)に関する方程式

を解いても良いのです。具体的に変数変換


を考えます。これは座標を\(\theta\)回転させた座標系への変換です。もちろん、\(\theta=0\)のとき、\(x=x’, y=y’\)が成立し、恒等変換となります。微分方程式は

であり、\((x’, y’)\)の空間で解いていきます。

さて、回転角\(\theta\)は任意ですが、式(15)の右辺が発散しないような\(\theta\)を選ぶ必要があります。
また、数値計算を行う上では、傾きがほとんど0に近いほうが有利です。そのため、方針としては現在の傾きを計算し、その傾き分だけ逆方向に座標回転させることでこれを実現させます。
すなわち、

分だけ回転させます。現在の傾きと座標軸の傾きを一致させることで傾きがほとんどゼロの座標軸に回転させるのです。

プログラム

以下にプログラムを示します。
必要な情報は\(f(x,y)\)の\(x, y\)方向の偏微分\(\frac{\partial f}{\partial x}, \frac{\partial f}{\partial y}\)の解析解です。関数名”fx”, ”fy”にそれぞれを入力してください。
計算は刻み幅制御の陽的ルンゲクッタ法を用いています。

具体例


具体的に関数\(f(x,y)=x^3-2xy+y^3\)の等高線を考えていきましょう。
関数の値がそれぞれ、\(0.01,1,2,5\)に等しいときの等高線を考えます。
具体的な初期値は、
\(\displaystyle
f(-2,1.3657)\approx 0.01,~f(-2,1.4646)\approx 1,~f(-2,1.5568)\approx 2,~f(-2,1.7977)\approx 5
\)
とします。
具体的に解いたのが左の図です。

”正解”としたgnuplotの等高線表示の結果と同じ結果が得られました。座標回転がうまく働いているようですね。微分の発散を抑えるだけでなく、数値計算量も抑えられているようです。

刻み幅\(h\)が正である場合と負である場合で進む方向が変わります。
また計算の途中で分岐がある場合、どうなるかは保証できません。

また、閉じた等高線が複数ある場合、その全てを計算することは本稿のプログラムだけではできません。初期値を求める別のプログラムを作成する必要があります。
更に、等高線がある値に等しい量を計算することはできません。初期値の値に等しい等高線を引きます。

追記


調べ物をしていたら

照井章, 佐々木建昭著『27.代数関数の陰関数描画について』
http://www.kurims.kyoto-u.ac.jp/~kyodo/kokyuroku/contents/pdf/0941-27.pdf
という素晴らしい資料がありました。

この文献では2つの方法が紹介されており、
1.\(f(x,y_0)=0\)を解いて\(x=x_1,x_2,\cdots,x_n\)を解く
2.陰関数定理から求める
という方法です。資料では、1.の方法を採用していましたが、本稿の方法は2.の方法でした。
1.はもれなく等高線を見つけることができますので、考えられる結果を得たい場合によい方法です。
一方、2.は全てのゼロ点を見つけることができませんが、計算が早くできます。別の方法などで注目したい等高線が分かっている場合に良い方法でしょう。

効率的に解くためには、場合によって使い分けることが一番です。

エアガンの集弾限界

理想的なエアガンとBB弾、それに風などが無い理想的な環境があったとします。

バレルとBB弾に直径差があることによって生じるばらつきは、

 30m先 | 直径3cm
 50m先 | 直径10cm

となります。この直径以下にばらつきを抑えることは構造的に不可能です(60cmのバレルで0.9Jで射出する場合)。
つまり、このばらつきよりも大きくばらつく場合、バレルとBB弾の直径差以外にばらつく原因が存在します。

バレルとBB弾に直径差があることによって生じるこのばらつきを軽減するには、
 ・バレルの直径を小さくする
 ・バレルを長くする
 ・バレル内部の素材とBB弾の反発が起こりにくいものにする
にすることで軽減されますが、他の要因が場合はこの限りではありません。


ここでの理想的とは、

・球の直径のばらつきは無い
・いつも同じ初速で撃ち出せる
・無風

とします。この場合でも同じ点に集弾することはありません。
なぜなら、バレルの直径とBB弾の直径が異なる為です。この差によってどの位集弾性が悪くなるのか、見積もりましょう。

これは構造的な問題であり、ばらつきの原因の中で取り除く事ができない1つの原因です。

ばらつく原因として、以下の三つが考えられますが、まず本稿では理想状態のばらつきを考えます。

ばらつきの種類 理想状態のばらつき(本稿) 製品誤差によるばらつき セッティングによるばらつき
BB弾重さ あり なし
BB弾の大きさ なし あり なし
回転のばらつき なし なし あり
手振れ なし なし あり

ばらつきが生じる原因


何故ばらつきが生じるかを考えましょう。
全くの理想であれば、銃口から同じ初速、角度、回転量で射出されたBB弾にばらつきが生じることはありません。
しかし、現実にはバレルの大きさとBB弾の大きさに差があります。これによって銃口から射出する時に進行方向と垂直な平面に速度を持つことになります。
銃口から出てきたBB弾の速度\(\mathbf{v}\)を以下のように表現します。

ここで、\(v_{z}, \mathbf{e}_{z}\)はバレルの方向に沿う速度、単位ベクトルで、
\(v_{x,y}, \mathbf{e}_{x,y}\)はバレルの方向に垂直な面への速度、単位ベクトルです。
すなわち、壁面に垂直な方向に対する速度ベクトルの大きさ\(v_\perp\)は、\(v_\perp=\sqrt{v_x^2+v_y^2}\)と表されます。

ばらつきが無くまっすぐ飛ぶのであれば\(v_{\perp}=0\)であり、そうでなければ\(v_{\perp}\ne 0\)です。

BB弾のばらつきが生じる原因は、射出方向に垂直な方向(横方向)に有限の速度が生じている、と仮定します。
この横方向の速度が生じる原因の一つは、ピストンでBB弾を空気で押す際に圧力が一定ではないとか、回転を掛けるためのゴムで生じる、などいろいろ考えられます。
ここでは、パッキンを通過して、もうこれ以上横揺れを増やすような原因が生じないのだ、と仮定します。壁の反射によって変化はあるものの、速度は減衰する一方であるとします。

それでは、定式化をしていきましょう。

定式化


横方向の速度はそれほど大きくないだろう、と予想するので空気抵抗は考えません。この条件のもと考えていきます。

バレルの直径\(d_s\)とBB弾の直径\(d_B\)の差があることによって、バレルに垂直な方向にBB弾が自由に進める距離\(l\)は

と書けます。\(n\)回目の衝突から\(n+1\)回目の衝突までにかかる時間\(t_n\)は、その間のBB弾の垂直方向の速度\(v_n\)に依存して、

と書けます。よって、\(N\)回衝突するのにかかる時間\(t\)はそれぞれの衝突までに掛かる時間を足し合わせればよいので、

です。バレルの内壁との衝突が起こることによってBB弾の速度の変化するとします。するとBB弾とバレルとの反発係数\(r\)を用いれば、

と表現されます。初速度を\(v_0\)と書いてこれを代入すれば、

と書くことが出来ます。ただし、\(r\ne 1\)です。\(r=1\)の場合は\(\displaystyle
t=\frac{l}{v_0}(N+1)
\)

です。\(r\ne 1\)の場合に\(N\)について変形すれば、

と書くことが出来ます。
初速\(v_0\)、反発係数\(r\), 銃口に到達するまでに掛かる時間\(t\)が分かれば、銃口から出た時のバレルに垂直な方向の速度\(v_\perp\)は式(5b)より、

と書けます。具体的に妥当な値を入れて射出時の、進行方向に垂直な方向の速度を計算してみましょう。

具体的な衝突回数の見積もり


60cmのバレル長を持ち、0.20gで0.9J程度の球が射出される場合、セットされた位置から射出までにかかる時間は約\(t=0.015\mathrm{[s]}\)と分かっていますのでこれを利用します(バレル内部の計算より)。
バレルBB弾の直径差は

とします。\(r\)はどうにかして求めることにして、\(v_0\)はとりあえず変数としましょう。

反発係数の見積もり


\(r\)を見積もります。
反発係数の見積もりは、実験して大雑把な値を見積もります。
\(r\)を見積もるために、家にある固い材質のものとBB弾とを衝突させて、高さを計測しました。以下の実測結果を得ました。

  • BB弾-フローリング 30cmから落として15cm上がる
  • BB弾-アクリル板 30cmから落として10cm上がる
  • BB弾-アルミ 10cmから落として6cmまで上がる

高さ\(h_0\)から静かに落として、反発した後に極大値をとるときの高さ\(h_1\)が分かっているとき、反発係数\(r\)は

と求められることが分かっているので、反発係数は

  • BB弾-フローリング  \(r\approx 0.71\)
  • BB弾-アクリル板  \(r\approx 0.57\)
  • BB弾-アルミ  \(r\approx 0.77\)

となります。バレル内部はどちらかといえば金属に近いと思いますので、\(r\approx 0.75\)と仮定して計算を進めていきます。

横軸に進行方向に垂直な方向の初速度…つまり、ホップを掛けるためのゴムなどで、進行方向ではない方向の速度のことを意味します…、縦軸に射出時の垂直な方向の速度をプロットしました。横方向の速度が仮に50m/sになっていても、10m/sになっていようともあまり変わらないことが分かります。なので、10m/sと仮定しましょう。
ばらつきの上限を与えるにはよい指標です(※1)。

仮に\(r=0.75, v_0=10\mathrm{[m/s]}\)とすると、\(N\approx 20.6\)と計算されます。つまり、20.6回バレルに衝突してから飛び出していくわけです。実際には衝突回数は整数しかありえないので、小数点以下を切り捨てて20回衝突が起こって飛び出していきます。
ただし計算を行う上では衝突回数を切り捨てると不連続性が発生しますので、小数点の衝突回数を認めることにします。

この場合、銃口から飛び出す際の、進行方向に垂直な方向の速度\(v_\perp\)は

を得ます。この垂直方向の速度は上下左右に振れる可能性がありますが、最も触れる場合は横方向でしょう。

30, 50m先の広がり具合を考えます。
着弾までの時間\(t_{30m}, t_{50m}\)はそれぞれ\(0.5, 2 [s]\)分かっているので(弾道計算(BB弾)の結果)、着弾時の直径\(d_{30m}, d_{50m}\)は

と求められます。ここで垂直方向の速度が非常に遅いので、空気抵抗は無視して考えています。BB弾の重さは関係ありません。軽いものほど空気抵抗を受けて減衰するので、ぶれは小さくなります。

空気抵抗を考慮すると、これよりは小さくなるということです。

実測定との比較


さて、30mチャレンジというものがあります30mチャレンジ公式ランキング 2016年度
この記録によると、30m先で直径15cm程度の範囲に収まるようです。
すなわち、バレル-BB弾の直径差による集弾性の悪化よりも大きな影響を及ぼす要因がある、ということです。

恐らくは回転量が一定ではないとか、BB弾の重さにばらつきがあったり、銃身の振動があったりという要因のほうが大きいということでしょう。

逆に言えば、スタンダードな大きさのバレルとBB弾を用いる限り、30m先で3cm以内に収めることは不可能、ということです。

注釈

※1
\(v_0\to \infty\)で射出時の速度\(v_\perp\)は発散しますが対数の発散です。非常にゆっくり発散するため、\(v_0=10\)であろうが\(v_0=50\)であろうがほとんど変化しません。どこかでカットオフ(これ以上はとらない上限)となるような値を取れればよいかもしれません。
豆知識ですが、対数の発散をほとんど無視する目的としてカットオフを設けるやり方は、量子力学の繰り込み論が有名ですね。だからと言って、エアガンの計算で量子力学が現れる!などは言わないでください。全く関係ありません

波動関数の規格化とは?

非相対論的シュレーディンガー方程式の解を規格化するお話

Q氏とR博士の会話1


Q「波動関数の規格化は、なぜ行われるのでしょうか?」

R「異なる量子数に属する状態間を、分かりやすく比較したいからだ。」

Q「なぜ規格化をすると比較ができるのでしょうか?」

R「例えば、有限区間内で定義される2つの状態があるとき、片方の状態の波動関数の値が1、もう片方が0.01という、どちらも単なる定数であることが分かったとしよう。」
R「その場合、規格化を行わないと、”値1を持つ状態のほうが多いから、もう片方は無視できる”、と直感的にしてしまうかもしれない。」
R「その点、規格化をすればどちらも同じ値にできるから、優劣がなく比較をすることが直感的にできて、都合がよい。」

Q「なるほど。わかりやすさ、ですか。」
Q「しかし規格化するといっても、何かを等しくして比較するって意味ですよね。何を等しくして比較するのですか?」

R「個数で規格化するのだ。」
R「1つの電子について解いているとき、波動関数が示す電子数を”1”であるようにするのだ。」

Q「個数で規格化するのですか。では、2個の電子を表す1つの波動関数があったらどうです?規格化は”1”ですか、”2”ですか?」

R「どちらも正解だ。1でもいいし、2でもいい。」
R「2個の電子が相関して離れない現象ならば、それは”1つのペア”を”1”にするのが良いだろう。だから、その意味ならば1が良い。」
R「しかし、電子”1つ”が重要ならば、電子1つの波動関数を”1”にするべきだろう。それが2つあるので合計数は2だ。」
R「定義が書いてあれば、それでよい。」

Q「見たい現象によって分けるのですか。規格化と言いながら、いかなる問題に対する規格は無いのですね。」

Q氏とR博士の会話2


Q「波動関数の規格化は”個数”で行われるんですよね。」

R「そうだ。」

Q「波動関数の”個数”は全存在確率密度を1にするようにすることで実現できます。」

R「そうだ。」

Q「つまり、波動関数の絶対値の二乗の、全空間に渡る積分が収束しなければならない。」

R「その通り。」

Q「ならば、自由粒子のような連続状態の規格化は不可能なのですか?」
Q「無限遠まで振動していますし、絶対値の2乗は定数です。」
Q「全位置空間に渡る積分は無限になり、規格化定数が零になります。」

R「可能だ。だが、君の疑問はもっともであり、計算自体は正しい。」
R「しかし、全存在確率密度による規格化には隠れた前提がある。」
R「それは”束縛状態であるならば、全存在確率密度で規格化すると都合がよい”ということだ。」
R「連続状態はその前提に入らない。なので全存在確率密度による規格化が計算不能でも、何も問題はない。」

Q「では連続状態の規格化とは、なんなのでしょう?」

R「計算が楽だとか、都合がよい…例えば、計算したら有限の値になる何か、があれば何で規格化しても良い。」
R「束縛状態において、全存在確率密度はたまたま採用されたのだ。」
R「本来、規格化は数学的な手順であり、”分かりやすい表現”の定数倍という表現をしたいがために生まれた。」
R「何で規格化するのか?が変わったところで、定数倍の違いしかないから、物理は変わらない。」
R「自由粒子のような連続状態の規格化は、通常、デルタ関数で行われる。」
R「しかし、例えば”指数関数の係数を1にするように決めた”と主張して、物理を作り上げても何も問題はない。」
R「この例は任意性があるので、褒められない規格化だがね。」

Q「規格化が連続状態と束縛状態で違ってもよい、ということですか。」
Q「しかし、それでは連続状態と束縛状態で規格化した関数が、その境界で不連続になってしまうかもしれませんが、問題ないのですか。」

R「問題ない。連続状態と束縛状態を別々に規格化すると約束するならば、不連続でも問題ない。」
R「一方で、束縛状態でも連続状態でも統一的な規格化の方法を採用すれば、不連続性は起こらず、すっきりする。」
R「そのような規格化方法として、例えば波動関数を単なる2乗で規格化する方法がある。」
R「この場合は波動関数を複素平面に解析接続しなければならなくなるが。」

マイナス×マイナス=プラスはなぜ?

負は、正の世界と負の世界を入れ替える因子である。

Q氏とR博士の会話1


Q「なぜ負に負を掛けると正になるのですか?」

R「ゆっくり確認していこう。まず、正の世界に正の数を掛けたら正になる。」

Q「納得します。」

R「では、正の世界に負の数を掛けたらどうなるかな?」

Q「負でしょう。」

R「そうだ。正の世界から見たら、負の数を掛けることで世界が反転するのだ。」
R「では負の世界から、世界を反転させる負の数を掛けたらどうなるかな?」

Q「負の反転は正ですから、正です。」

R「その通り。なので負に負を掛けると正なのだ。」

Q氏とR博士の会話2


Q「正の数に、正を掛けたら正になるのですよね。」

R「そうだ。」

Q「負の数に、負を掛けたら正になるのですよね。」

R「そうだ。」

Q「自分の世界の数同士をかけたらどちらも正はおかしくないですか?」
Q「正×正=正に、負×負=負になるのが、道理でしょう?」

R「その通りだ。君の疑問は全く正しい。だが、人類の都合で決められてしまったのだ。」
R「正の数しか無かった頃、人類は”正×正=正”を基準にした。その後に”数”を負に拡張したのだ。」
R「負の数に広げても”正×正=正”を保ちつつ、”正×負=負”または”負×正=負”になるように拡張したいと考えた。」
R「そう制約したから、”負×負=正”となってしまったのだ。負は正の世界と負の世界を入れ替える因子である、と。」

Q「なるほど。もし仮に”負×負=負”を基準に選んでいたら、正の数が正の世界と負の世界を入れ替える因子となるわけですか。」

R「その通りだ。負×負=負”を基準とするならば、”正×負=正”または”負×正=正”となり、”正×正=負”になる。こういった世界を作ることも可能なのだ。」
R「どちらの考え方も間違いではない。なぜ負×負=正なのかの理由を突き詰めるならば、『正×正=正に選んだから』だ。」

Q「ならば、無理やり”正×正=正”と”負×負=負”を満たすように強引に数学を組み立てたら、”正×負”はどうなるのでしょう?」

R「未定だ。定義できない。正×負の単位を数えるための、新たな数の拡張が必要になるだろうね。」

秋間補間

秋間補間と呼ばれる補間方法があります。

こんな補間結果となります(akimaとmodified akimaの線が秋間補間の結果です)。

秋間補間の特徴は、

  • 与えられた点を必ず通る
  • オーバーシュートが発生しにくい
  • 補間したい点から遠い補間点の影響を受けない
  • 導関数の連続性は保証されない
  • 元のデータ点を大きく逸脱するような補間は苦手(線形的な補間に近い)
  • 極値におけるデータ点など、元の関数の特徴的な値が分かっているときに良い補間結果を与える

ということです。オーバーシュートとは、元のデータが急峻に変化した場合、その前後の補間結果が振動してしまい大きく外れてしまうことを意味します([5]にオーバーシュートの様子を見ることができます)。
ただし、3次スプライン補間のような微分の連続性が保証されません。

それでは計算方法を追っていきます。

導出


秋間補間は1970年に秋間浩によって提案された補間方法です[1]。地図の等高線を表示するために生成されたようです。その後修正された秋間補間が紹介されました[2]。

具体的な計算方法を見ていきます。
補間は補間したい点を中心に前後のデータ点を用いるので、両端では特別な扱いが必要です。それを区別するために領域を3つに分けて考えます。具体的には下の図のような感じです。

領域IIにおける補間


まずは領域IIを考えます。簡単に考えるために、以下の図のように補間に用いるための6点\((x_i, y_i), i=0,5\)があった時、中間の領域である\(x_2\le x\le x_3\)の中の値を補間して求めたいと考えます。

補間は3次関数で行います。すなわち\(x_2\le x\le x_3\)において

の形で補間を行います。係数\(a_0, a_1, a_2, a_3\)は補間に用いる6点のデータ点から計算できて、

の通り求められます。ここで、\(q_1, q_2, m_1,\cdots, m_5\)は


の通り求められます。\(w_1, w_2\)は使用する秋間補間の種類によって変わります。通常の秋間補間では、

の通り計算し[1,3,4]、修正秋間補間においては

で計算します[2,3,4]。

領域I, IIIにおける補間


領域I, IIIについては以下の関係を満たすように、本来は存在しない補間用のデータ点を予測します。
領域I:最も左端のデータ点を\(x_3\)、その右隣のデータ点を\(x_4, x_5\)と置きます。すると、\(x_3\)より左のデータ点\(x_2, x_1\)を関係式(7)より予測します。
領域III:最も右端のデータ点を\(x_3\)、その左隣のデータ点を\(x_2, x_1\)と置きます。すると、\(x_3\)より右のデータ点\(x_4, x_5\)を関係式(7)より予測します。

ここで、\(m_{ij}\)は

を意味します。すると、秋間補間で用いる(本来は存在しない)端のデータ点は、

領域I)

領域III)

と求められます。

プログラム例


上の式をプログラムすれば、このようになるかと思います。
注意してプログラムを作らないとゼロ割が発生します。\(x_i\ne x_j\)を仮定していれば、式(5),(6)の\(w_1, w_2\)を計算する際にこれが発生します。なので、この回避を取り入れてプログラムを作成しています。

計算例


参考文献[4]に秋間補間と修正秋間補間の振る舞いの例があります。それを本稿のプログラムを用いて計算すると同じ振る舞いを示すことが分かります。

それでは、3次スプライン補間と秋間補間、修正された秋間補間を比べて、秋間補間にどのような特徴があるかみましょう。

オーバーシュートの抑圧


本当の答えを定義することは難しいですが、変化する点(\(x=2.8\sim 3.2\))近傍で直線とすれば傾きは5になります。一階導関数はどの補間方法も連続(右から極限と左から極限が一致の意味で)ですが、秋間補間は二階導関数は連続ではありません。3次スプラインは二階導関数まで連続であるという仮定の下での補間なので、一応連続的な変化をしています。

秋間補間では3次スプライン補間よりはオーバーシュートが抑えられているものの、近傍 \((x=2.5, 3.5)\)において外れた値となっています。
修正された秋間補間では、その点が改善されオーバーシュートが発生していないことが分かります。

振動関数の補間

振動関数を補間する例として、元の関数が\(y=\sin(x)\)だったとします。

sin関数

この場合は元の関数が十分滑らかであり、高階微分も連続であることが分かっています。与えられたデータ点も誤差がないので、3次スプラインが有効でしょう。
秋間補間では本来極値をとるはずの\(x=n\pi/2,~(nは整数)\)あたりで値があまり大きくならないことから、かなり局所的な補間である特徴が分かります。すなわち、大きくオーバーシュートしない反面、予測は得意ではないということです。これを回避するためには、元の関数の極値の値など特徴的なデータが必要であることが分かります。

二次元平面の補間は本稿には載せていませんが、[4,6]にどのような違いが表れるか図があるので、興味がある方は参照してみてください。
確かに、本来の目的である地図の等高線作成には大きな効果が期待できそうなことが分かります。おそらく、取得するデータ点も標高など、地形の特徴的な点におけるデータを集めることが主となるので、秋間補間との相性が良いのでしょう。

参考文献


[1]Akima H. A New Method of Interpolation and Smooth Curve Fitting Based on Local
Procedures. Journal of the ACM, 17(4), 589 ― 602. (1970).
doi:10.1145/321607.321609

[2]Akima H. A method of univariate interpolation that has the accuracy of a third-degree
polynomial. ACM Transactions on Mathematical Software, 17(3), 341 ― 366. (1991),
doi:10.1145/114697.116810

[3]
makima 修正Akima区分的3次エルミート内挿 -MathWorks
interp1 1次データ内挿(テーブルルックアップ) -MathWorks

[4]Cleve Moler, Makima Piecewise Cubic Interpolation

[5]スプラインによる補間

[6]格子点に内挿して等値線をひく

emacsのメニューバー、ツールバーを消去する

emacsの操作に慣れてくると、emacs画面の上部にあるパネルやツールバーを使わなくなります。
大事な画面のリソースが減ってしまうので、これを消すことを考えます。

emacsの設定ファイル(~/.emacs)にメニューバーとツールバーの消去する設定

(menu-bar-mode -1)
(tool-bar-mode 0)

を記述することで消せます。書いてemacsを再度開くと、以下の通りとなります。

環境:Linux-mint 20
Emacs: GNU Emacs 26.3 (build 2, x86_64-pc-linux-gnu, GTK+ Version 3.24.14) of 2020-03-26, modified by Debian

メニューバーとツールバー有り

メニューバーとツールバー無し

参考

メニューバーの削除
Remove the Emacs menu bar

ツールバーの削除
in GNU Emacs OSX, how to hide title bar? stack overflow

Lenovo thinkpad X1とeGPUの接続時のフリーズ

現在、Lenovo社から発売されているThinkpadにeGPUを接続したとき、30分か1時間ほど経過するとフリーズする、という状態が頻発しました。

解決法はeGPU Freezing up every 30 minutes.
にあるように、

デスクトップを右クリックして”NVidia コントロールパネル”を表示
 1. 『3D設定』→『3D設定の管理』を選択
 2. 『グローバル設定』タブで『高パフォーマンスNVIDIAプロセッサ』を選択
 3. 『電源管理モード』を『パフォーマンス最大化を優先』に変更

とします。そうしたところ、現在安定しています。

環境と状況


環境
 品名: ThinkPad X1 Carbon (2019)
 OS:  Windows 10 Pro 64-bit
 eGPU Box: Razer Core X
 GPU: GeForce GTX 1660 SUPER VENTUS XS OC
 接続: ThinkPad Thunderbolt 3 Workstation ドック 2を利用したThunderbolt接続
 出力: eGPU BoxのGPUから外部ディスプレイに出力

状況と解決

eGPUを接続して30分~2時間は問題なくGPUが利用できるのですが、その後前兆もなくフリーズします。このフリーズはマウス、キーボードの入力を全く受け付けなくなる酷いフリーズです。電源ボタンを長押しして強制終了する以外にやれることがありません。
強制終了後に起動したのちに動く時間を利用して、再起動、リセット、GPUドライバの更新などを試しましたが状況が全く変わりませんでした。

いろいろ調べたところ、eGPU Freezing up every 30 minutes.に行きつき、上のような設定を行ったところ、今のところ安定して動いています。

論文に使える図をgnuplotで作る

本稿の目的は、論文に用いることができるレベルの図をgnuplotで制作することです。
gnuplotはVersion 5.2 patchlevel 8を利用しています。

方針


本稿の方針としては、実際の論文からとってきた図を再現することです。
元データはないので、あくまで枠線や、フォントなどを合わせることを目的とします。

対象とする図は以下の[1]の図3と[2]の図3にします。
これらを選んだ理由は特にないですが、これらを再現することを目標にします。

[1]の図3を再現


ではまず、Phys. Rev. A 102, 033316 (2020)の図3を再現しましょう。

[1]の図3

[1]の図3をgnuplotで再現した図

gnuplotで用いた全コードは以下に折りたたんでおきます。

データは適当に作成しました。データのダウンロードは以下からどうぞ。
https://slpr.sakura.ne.jp/qp/supplement_data/PRA102_033107_Fig3likedata.tar.gz

細かな説明をしていきます。

  • 線の色の指定
    線の色はlinetypeで決定していきます。

    set linetype  1 lc rgb '#0072bd' lw 1 # blue

    の意味は、gnuplotでプロットした時の1番目の線を、カラーコード’#0072bd’, 線の太さ”1″とすることを指定しています。プロット時に”lc 1″としてもこの1番目の色に変更できるようになります。linetype 6まで同じです。

  • 複数画面の表示位置指定
    続いて、

    # left top figure(1) position
    lx1=0.14; ly1=0.87
    rx1=0.47; ry1=0.57

    ~~~

    set lmargin screen lx1
    set tmargin screen ly1
    set rmargin screen rx1
    set bmargin screen ry1

    の部分について説明します。個別画面の左上の座標\((lx1, ly1)\)と右下の座標\((rx1, ry1)\)を決定します。座標の原点は全体の左下\((0,0)\), 全体の右上\((1,1)\)としたときに決められます。実際にその値を適用しているのがlmargin, tmarginなどです。詳しい説明は、gnuplotで複数の図を載せるを参照してください。

  • eps出力のためのターミナル指定
  • if(@ARG1 == 1){
      set terminal postscript eps size 5,5 dashed enhanced color font 'Roman,20'
      set output "PRA102_033107_Fig3.eps"
    }

    の部分は、出力ターミナルの設定を行っています。条件分岐は引数によってターミナルに出力するのか、epsとして画像出力するのかを指定するためです。使いやすさのためです。重要なのは、 size 5,5 の箇所です。ここを正方形にしておかないとxyの比が変わってしまうので扱いづらくなります。
    eps画像の切り抜き(Boundary boxの変更)はgsviewなどでできるので、画面四方の空白はいくらあってもよいです。そのサイズをgnuplot上で合わせるよりも、gnuplotで思った通りの出力をすること優先しましょう。
    フォントも重要ですね。経験上、2×2の画像であればsize5,5の時フォントサイズ20よいと思います。論文の本文はRomanフォントが使われるため、図もRomanにしておいたほうがよいでしょう。

  • 個別画面の細かい指定
    1. タイトルの指定
      set title "{/Gothic-bold=24 Linear scale}" offset 0,0

      タイトルを指定します。図のLinear scaleと書かれた太字です。[1]の図を見ると、どうもゴシック体の太字で書いているようです。なので、フォントをGothic-boldとして、文字の大きさを24に変更します。

    2. 凡例の指定
    3. set key reverse samplen 1 width 0

      凡例はkeyで指定します。reverseは凡例の中の”文字-線の種類”の表示順を”線の種類-文字”に変更します。samplenは、凡例内の線の長さを変更する(デフォルトで4)、wigthは線の種類と文字の間隔を指定するパラメータとして考えればよいでしょう。具体的な値は実際に出力させて何度も気に入るまでトライアンドエラーしています。

    4. x,y軸の設定
    5. unset log x
      set xr[-2:52]
      set xtics 10 nomirror
      set mxtics 1
      set format x "%2.0f"
      set xl "{/Roman-Italic {/=20 tJ}}/{/=14 @^{/Symbol=14 -}{/Roman-Italic {/=20 h}}}"

      軸設定をまとめて説明しますが、特に特殊文字の出力が難しいと思うのでここに注目します。プランク定数は特殊文字であり、postscriptターミナルでは出力することができないと思います。なので、文字hとハイフンを組み合わせて出力しています。私が思いついたわけではなくググりました。[4]に詳細が書いてありまして、それから持ってきました。文字サイズが違うので適当に合わせました。
      y軸も同じですので、省略します

    6. グラフのプロット
    7. plot "hJ5_fit.d" u 1:2 w l lw 2 lc 5 dashtype 2 ti "{/Gothic-Condensed linear fit}",\
           "hJ5.d" u 1:2 w l lw 2 lc 1 dashtype 1 ti "{/Roman-Italic M}=24   ",\
           "hJ5.d" u 1:2 w l lw 2 lc 2 dashtype 2 ti "{/Roman-Italic M}=32   ",\
           "hJ5.d" u 1:2 w l lw 2 lc 3 dashtype 4 ti "{/Roman-Italic M}=64   ",\
           "hJ5.d" u 1:2 w l lw 2 lc 4 dashtype 3 ti "{/Roman-Italic M}=128 ",\
           "hJ5.d" u 1:2 w l lw 2 lc 5 dashtype 5 ti "{/Roman-Italic M}=256 "

      最後にグラフをプロットします。dashtypeは実線、破線を決定するパラメータです。私の環境では実線が1, 破線は2以降になっていますが、環境によって違うかもしれません。各自合わせましょう。
      また、linear fitという文字は、なんというかギュッとまとまっています。しかもRoman体ではなくゴシック体です。これを指定するには、

      ti "{/Gothic-Condensed linear fit}"

      と入力することで実現可能です。

以上でおおよその説明が終わりです。

重要なのはトライアンドエラーです。ターミナルによって表示が異なるので、postscriptで出力すると決めたら確認はpostscriptで行いましょう。実際に最終出力図をqtターミナルの表示上では、このようになります。

ベクタ形式であるEPS画像のサイズは102.5kBです。論文に投稿する場合でも許容できるサイズではないでしょうか。

ラスタライズ形式であるPNGをimagemagickを用いて

convert -density 400x400 -background white -flatten -alpha off PRA102_033316_Fig3.eps PRA102_033316_Fig3.png

で変換すると、91.8kBです。この程度の差であれば、ベクタ形式のほうが良いでしょう。

[2]の図3を再現


続いてPhys. Rev. A 102, 033107 (2020)の図3を再現しましょう。
大きな違いはカラーマップを利用した図であることです。

[2]の図3

[2]の図3をgnuplotで再現した図

データは適当に作成しました。データのダウンロードは以下からどうぞ。
https://slpr.sakura.ne.jp/qp/supplement_data/PRA102_033107_Fig3likedata.tar.gz

※ギリシャ文字Φは、

{/Symbol-Oblique f}

とすることでギリシャ文字のイタリックを用いることができます。

[1]と同じ部分がかなりあるので、違う部分だけを説明していきます。

  • カラーマップの設定
    見た目で[3]の’hot’と呼ばれるカラーセットだとわかったので、それを用います。

    set palette rgb 21,22,23 negative # colormap setting

    ただし、カラーマップの色が反対になっているので、それを指定するnegativeを追加しています。

  • カラーボックスの設定
    unset colorbox

    ~~~~

    set cbtics 0.5 scale 0 # scale 0, to delete tics line in colorbox
    set cbr[-2:0]
    set colorbox user origin 0.91,0.35 size 0.02,0.25 border # change size of colorbox

    カラーボックスに注目すると、左図(a)はカラーボックスがないですが、右図(b)にはあります。なので、カラーボックスを出力しないunset colorboxを指定しています。
    また、(b)のカラーボックスはカラーの値を示すところにticsが入ってないため、cbticsのscale(ticsの線の大きさ)を0にしています。
    さらに見ると、カラーボックスの大きさが少し小さいです。なので、カラーボックスの大きさ位置を指定するために、set colorbox user origin …を指定しています。

  • テキストボックスの表示
    set label tc rgb "black" at graph 0.02, graph 0.95 front "(a)"

    グラフ上のテキスト(a), (b)は、テキストボックスを追加して画面上に表示します。上のような指定で表示できます。

ベクタ形式であるEPS画像のサイズは752kBです。論文に投稿する画像としては、許容できるといえばできるサイズですが、若干大きいです。

ラスタライズ形式であるPNGをimagemagickを用いて

convert -density 400x400 -background white -flatten -alpha off PRA102_033107_Fig3.eps PRA102_033107_Fig3.png

で変換すると296kBです。2倍少し容量が小さくなるので論文全体の画像数が多いのならば、ラスタライズ形式の方が良いかもしれません。共著者と相談するとよいでしょう。

以上で説明を終わります。

参考文献


[1]S. Goto and I. Danshita, “Measurement-induced transitions of the entanglement scaling law in ultracold gases with controllable dissipation”, Phys. Rev. A 102, 033316 (2020)

[2]Z. Chen and F. He, “Interference of nuclear wave packets carrying different angular momenta in the dissociation of H2+ in strong circularly polarized laser pulses”, Phys. Rev. A 102, 033107 (2020)

[3]color mapについて
http://www.gnuplotting.org/?s=color

[4]6.4 What if I need h-bar (Planck’s constant)? -gnuplot FAQ
http://www.gnuplot.info/faq/#x1-610006.4

そのほか有用なページ

gnuplotでoption詰め込んで論文用figureを作ってみた[1] -qiita, @ucc_white
https://qiita.com/ucc_white/items/544b85e009a58d046d3c

Gnuplotのいろは -東京工業大学 科学技術創成研究院 佐藤千明研究室
http://www.csato.pi.titech.ac.jp/From2014/software/how-to-gnuplot.html

連成振動

弦を考えます。弦は重さ\(M\), 長さ\(L\), ばね定数\(K\)を持ちます。

この弦を波動方程式に頼らないで定式化しましょう。
波動方程式は、2次元の波動のうち、1次元分の方向を\(y=y(x, t)\)として表現するため、x軸にとって多価関数となる弦を表現することが出来ません。物性物理学等で格子の振動を考える場合は、結晶を構成する原子が平衡位置から殆ど動かないのでこういった問題は起こりません(起こったとしても無視される位少ないので、大多数の振る舞いを考える場合には問題になりません)。
よって多次元の問題を考えるには波動方程式を元にして考えるのは不適切です。そのため、弦を単なる連成振動として考えて、完全に2次元の問題として扱いたいと思います。

図のような連成振動のモデルを考えます。

質量\(m\)を持つ\(i\)番目の質点が、位置ベクトル\(\mathbf{x}_i\)で指定され、隣り合う質点との間が、ばね定数\(k\)、自然長\(l\)を持つ理想的なばねで繋がれているとします\((i=1,2,\cdots, N)\)。この時、運動方程式を立てると
\begin{eqnarray}
m\frac{d^2 \mathbf{x}_i}{dt^2}&=&
+k(|\mathbf{x}_{i+1}-\mathbf{x}_{i}|-l) \frac{\mathbf{x}_{i+1}-\mathbf{x}_{i}}{|\mathbf{x}_{i+1}-\mathbf{x}_{i}|}
-k(|\mathbf{x}_{i}-\mathbf{x}_{i-1}|-l)\frac{\mathbf{x}_i-\mathbf{x}_{i-1}}{|\mathbf{x}_i-\mathbf{x}_{i-1}|}
\end{eqnarray}
となります。全く同じですが、見通しをよくするため式変形すれば、
\begin{eqnarray}
m\frac{d^2 \mathbf{x}_i}{dt^2}&=& k\Bigl[A_i \mathbf{x}_{i+1}-\bigl(A_i+A_{i-1}\bigr)\mathbf{x}_{i}+A_{i-1} \mathbf{x}_{i-1}\Bigr]
\end{eqnarray}
と表せます。ここで、
\begin{equation}
A_i = \frac{|\mathbf{x}_{i+1}-\mathbf{x}_{i}|-l}{|\mathbf{x}_{i+1}-\mathbf{x}_{i}|},~~(i=1,\cdots, N, A_N = 0)
\end{equation}
と置きました。

さて、こうしたモデルを考えた時に、連成振動として考えた\(N\)分割した1つ1つの小さいばね(”ミクロなばね”と名づけます)とその弦を1つにみなした時の大きなばね(”マクロなばね”と名づけます)の間の関係式を考えましょう。
仮定として、等間隔に分割し、弦の密度は一様であるとします。

この条件の下で、
マクロなばねの重さ\(M\), 自然長\(L\), ばね定数\(K\)

ミクロなばねの重さ\(m\), 自然長\(l\), ばね定数\(k\)
の間の関係式を導出します。

結論として、マクロとミクロなばねの間に成り立つ関係は
\begin{eqnarray}
L&=&N l \\
M&=&N m \\
K &=& k/N
\end{eqnarray}
となります。

自然長

自然長\(L\)の弦を\(N\)分割するので、
\begin{equation}
L=N l
\end{equation}
の関係があります。

重さ

重さ\(M\)の弦を\(N\)分割するので、
\begin{equation}
M=N m
\end{equation}
と書けます。

ばね定数

最後にばね定数の関係を考えましょう。この関係は、弾性エネルギーを考えることで導くことができます。
マクロなばねとして全体を見た場合に長さ\(D\)だけ縮んだ時に蓄えられる弾性エネルギーと
同じだけ縮んだ時にミクロな\(N\)個のばねに蓄えられる弾性エネルギーが同じであってほしいと要請します。すなわち、ミクロなばねの縮みの合計がマクロなばねの縮み\(D\)に等しいとすると、
\begin{equation}
D=N d
\end{equation}
と書けます。
実は、大胆な過程をしないとマクロなばねとミクロなばねの関係を求めることが出来ません。
なぜなら、マクロなばねはいつでも線形応答をするかのように現在考えていますが、ミクロなばねでは場所ごとに違っても良いという様に定式化しています。なのでもともと無理な話なのです。
ではどうやって妥当な関係式を導出すればよいかといえば、ばねが非常にゆっくり運動が行われるときを仮定するのです。ゆっくり動くときの振る舞いは全てのばねが同じ動きをするとして差し支えないでしょう。そうして定式化します。
よって、
\begin{eqnarray}
\frac{1}{2}K D^2 &=& N\cdot \frac{1}{2} k d^2 \nonumber \\
\to ~~ K &=& k/N
\end{eqnarray}
という関係式が導けます。

あとは数値計算を行うだけです。陽的ルンゲクッタ法で解いたプログラムを以下に示します。

https://slpr.sakura.ne.jp/qp/supplement_data/wave_particles.tar.gz

実行は、

$ gfortran rkf45.f90 main.f90
$ ./a.out

でokだと思います。初期状態や、重さや自然長はinputファイルを見てください。

実行後、動画をgnuplot上で出力したければ、

$ gnuplot
$ call "movie.plt" 100

としてください。100は表示する時間ステップを表します。この引数の値の最大はinputファイルの中のNtまでです。

数値計算結果


実際に数値計算を行う条件として、自然長や重力があると波動方程式と若干異なる振る舞いになってしまうので、対応を見る上ではなくしています。すなわち、自然長\(L=0\)、重力\(g=0\)として計算します。

三角パルス

三角形の波がパルス状に存在するとします。
そして初期状態として、初速度がゼロであることを仮定しましょう。

この場合、波動方程式の考察から、
同じ関数であらわされる右向きの波\(f(x-vt)\)と左向きの波\(f(x+vt)\)の重ね合わせとして書かれることが分かっています。すなわち、解\(y=y(x,t)\)は
\begin{equation}
y(x,t)=f(x-vt)+f(x+vt)
\end{equation}
と書かれます。図で表せば以下のように時間発展していくことが期待されます。

これが本稿の連成振動モデルで実際にそうなっているのかを確かめるのは良い問いかけでしょう。
実際に連成振動を数値的に解くと以下の通りになります。

ゆっくりと時間発展を見ていくと

となります。明らかに右向きに進む波と左向きに進む波として分離します。分離した結果から、初期状態がそれらの波の重ね合わせで表現されていたことが分かるでしょう。

矩形パルス

矩形パルスを考えてみます。実際に数値計算をするとこんな感じです。

結局右向きと左向きの重ね合わせとしてあらわされることも確認できますが、決して単純なそれらの足し合わせで書かれているわけではないことがわかります。この振る舞いは、波動方程式に対応するように連成振動モデルを作っていないことから由来します。波動方程式としてあらわした二階微分の偏微分方程式の場合は多価関数をとることができません。なので縦に弦が伸びる場合、間に質点を置くことができないので質点の位置がいきなり飛ぶことになります。しかし、連成振動モデルではその制約をして初期状態を準備しているわけではありませんのでその違いが出て来ます。図で表せば、以下のようになります。

両者の違いがあるので、波動方程式の結論である
\begin{equation}
y(x,t)=f(x-vt)+f(x+vt)
\end{equation}
が成り立っていない、と解釈することができます。

少しずつ時間発展をみるとこの通りになります。

三角形の波

端まで三角形のパルスであった場合にどのように振る舞うかを見てみましょう。特に面白い振る舞いがあるわけではなく、上のパルスの、右向きと左向きの波の重ね合わせで説明ができるので、単にシミュレーション結果を載せるだけにとどめます。実際に計算してみると以下のようになります。

星形

初期状態が星型である場合も見てみましょう。振る舞いも同じです。

注意点としては、実際にこの振る舞いは起こりえないことに注意しておきましょう。なぜなら、本当の弦で行った場合、重なる場所があると弦の衝突が起こります。しかし、この連成振動モデルでは計算上、衝突は起こりません。

偏りがある三角形のパルス

偏りがある場合も波動方程式では記述することができません。こんな感じの初期状態です。

これを実際にシミュレートすると以下の通りになります。

3次元の場合

3次元の場合も載せておきます。ふるまい自体はそんなに変わるわけではないので、シミュレーション結果を載せるだけにとどめておきます。