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

位相速度、群速度とは!?調べてみました!!!!

近頃SNSやインターネットでよく見る言葉、「位相速度」「群速度」

元ネタや原理を分かった気になって使っている方、多いのではないでしょうか。
せっかくなので調べてみました!!!

  1. そもそも「位相速度」「群速度」って何?
  2. 学会騒然。ある波数近傍だけ値を持つ場合にテーラー展開!?
  3. 位相速度はどうやって!?驚きの仮定
  4. 簡単だと思いました?注意です!
  5. いろいろなところに群速度に関する議論が!
  6. 最後に
  7. Note

そもそも「位相速度」「群速度」って何?


なにはともあれwikipedia!!
関連するキーワードをササッと知るには強い味方です。
「位相速度」「群速度」を調べるとこのように説明されていました。

群速度(ぐんそくど、英: group velocity)とは、複数の波を重ね合わせた時にその全体(波束)が移動する速度のことである。

群速度 -wikipedia

位相速度(いそうそくど、英語: phase velocity)は、位相、すなわち波の山や谷の特定の位置が移動する速度のことである。

位相速度 -wikipedia

どうやら、波束波の速度を表現する際に都合が良い言葉なようです。

波束とは何なのでしょう…?また変な言葉が出てきました。wikipediaの出番です!!

波束(はそく、英: wave packet, wave train)は、局所的に存在する波うち/波動であり、移動する1個の波動の塊のようにふるまう。

波束 -wikipedia

うーん、つまり一回だけ波打って、ある程度まとまった波みたいなものですかね。
水面に水が一滴だけ落ちたときに発生する波みたいなものでしょう。

具体的にどういう意味なのでしょうか。気になりますよね!

皆さんも群速度、位相速度の導出が気になったら「tweet」、お願いします!

学会騒然。ある波数近傍だけ値を持つ場合にテーラー展開!?


角周波数が波数の関数として書けている場合、分散関係があるといいます。
つまり、角周波数\(\omega\)が波数\(k\)に依存して\(\omega=\omega(k)\)と書けている場合、波形のほとんどは

とように書けます!!
本来、位置\(x\)と時間\(t\)で書ける関数は、\(k\)と\(w\)の二重積分で表現されなければなりませんが、分散関係のおかげで\(\omega\)が強制的に決まってしまうので、一重積分で良くなるのですね!
うーん、とても嬉しい!!!!

…でも、あまりに一般的過ぎて分かりませんね。
しかし!!
応用で重要なのは、ゆっくり振動する包絡線早く振動する波掛け算で表されていることが多いのです。

つまり、これから考える波が持つ波数は、\(k= k_0\)の周りしか値を持たないのです。ここで、早く振動する波の波数を\(k_0\)と置きました。

\(k_0=5\)の場合、これから考えていく\(f(k)\)はこんな感じの特徴を持っています!ババン!

こんな特徴を持つのならば、\(k=k_0\)以外ではどうせ\(f(k)\to 0\)となるので、考えてもしょうがないっていう近似が使えそうですよね。

そ こ で !

波数\(k=k_0\)周りでいろいろテーラー展開して使用していきます。つまり、\(k=k_0\)の近傍さえ合っていれば、\(k_0\)以外でどんな振る舞いになろうとも\(f(k)\)がゼロになるので大丈夫でしょう!ということです。

早速、分散関係についてテーラー展開を行うと

となります。式を簡単に書くために\(\omega'(k_0)\equiv \frac{d\omega}{dk}\bigr|_{k=k_0}\)と置きました。

計算を進めると、

より、

のように書けます。つまり、式(4)が言っていることは、
時刻\(t\)の波形とは、初期状態\(t = 0\)の波形 \(f(x,0)\)を\(\omega'(k_0)t\)だけ平行移動させた波に、時間だけに依存する位相に関する変化分 \(e^{−i[\omega(k_0)−k_0 \omega'(k_0)]t}\)が掛け合わされたもの、と表すことができる。と言っています。

以上から、初期状態の波が時刻\(t\)に至るまでの速度は\(\omega'(k_0)\)で決まることが分かり、これは
群速度
と呼ばれます!!

位相速度はどうやって!?驚きの仮定


さて、ここまでで群速度が求まりましたが、位相速度が求められていません。
そこで、もう少し具体的な波の形を定義しましょう。
波が、

と書けている場合を想定します。ここで、振幅\(A(x,t)\)は実数値関数で波の包絡線を表しています。また\(A(x,t)\)の変化は \(k_0\) に比べて非常にゆっくり変化しているものとします。
つまり、

今まで考えてきた仮定が全部使える

んですね!

式(5)に代入しますと

となります。つまり、振幅の変化はいつも初期状態の単なる平行移動として表現することができ、位相の変化も簡単に表せられます。

そこで、振幅が平行移動する速度を群速度\(v_g\)、位相が平行移動していく速度を位相速度 \(v_p\) として定義すると、 式(6c) よりそれぞれ

と定義してしまうのが良さそうです!

簡単だと思いました?注意です!


注意しなければならないのは、群速度と位相速度の概念は、

振幅と位相を記述する部分がはっきりと区別できるような状況でなければ、定義できない

ということです!振幅部が早く振動していたら、これまでの議論が使えなくなることに注意しましょう。
その場合、群速度と位相速度の概念が変わってしまうので、意味をなさなくなってしまいます。

いろいろなところに群速度に関する議論が!


様々なところにリンクがあります。やはり皆、群速度の概念を理解するのに苦労しているようですね。

群速度の謎 -平林 浩一, (C) 2006

「yam@広島大」物理化学Monographシリーズ ”3. 物体の速度と物質波の速度 -E=hνの本質の理解-”

最後に


いかがでしたか?気に入ったなら是非、記事の下にある「like!」ボタンや「tweet」して皆さんに共有して下さいね!
それでは皆さん、良い物理ライフを!

Note


本記事の大枠の構成はいかがでしたか? -ニコニコ大百科を参照しています。

再学習用の電磁気学

電磁気学を可能な限り短くまとめました。電磁気学を一度学んだ方の再学習用を想定しています。
歴史を追うような考えをしていません。なので、証明には「なるからなる」を多く取り入れています。
経緯を知りたい方は各種参考書をご覧ください。

https://slpr.sakura.ne.jp/qp/wp-content/uploads/2021/02/electromagnetism.pdf

electromagnetism

間違いがありましたら、コメントなどいただけると助かります。

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

問題


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

問題設定


非相対論の範囲で、水素原子の電子を考えます。
電子の満たすシュレーディンガー方程式は、プランク定数\(\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})\)の項が出てきてポテンシャルエネルギーの項を打ち消さなければなりません。
波動関数に課されるこの条件は、カスプ条件と呼ばれます。

エアガンの集弾限界

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

この場合、

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

以下のばらつきに抑えることは構造的に不可能です(60cmのバレルで0.9Jで射出する場合)。

バレルと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「この場合は波動関数を複素平面に解析接続しなければならなくなるが。」

連成振動

弦を考えます。弦は重さ\(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次元の場合も載せておきます。ふるまい自体はそんなに変わるわけではないので、シミュレーション結果を載せるだけにとどめておきます。

熱力学

  1. 準静的な過程
  2. 平衡状態と最適なエネルギー
    1. 熱のやりとり(\(d’Q=TdS\))と体積変化(\(dV\))を見たい時
    2. 熱のやりとり(\(d’Q=TdS\))と圧力変化(\(dp\))を見たい時
    3. 温度変化(\(dT\))と体積変化(\(dV\))を見たい時
    4. 温度変化(\(dT\))と圧力変化(\(dp\))を見たい時
  3. 熱力学変数の各種関係式
  4. 熱容量
    1. 定積熱容量
    2. 定圧熱容量
    3. 定積熱容量と定圧熱容量の関係
    4. \(dp\)の消去
    5. \(dV\)の消去
  5. エントロピーについて
  6. さいごに
  7. 注記
  8. 参考文献

準静的な過程


気体を状態変化させることを考えます。
状態変化とは、圧力や体積、温度、粒子数などが変化することを指します。

気体の持つエネルギーの微小変化\(dE\)は、外部からの気体へ与えられたエネルギー\(d’Q, d’W, d’G\)を用いて
\begin{align}
dE = d’Q+d’W+d’G
\end{align}
と書くことが出来ます。この式はエネルギー保存を表しており、左辺が気体の持つエネルギーの変化分、右辺が気体に加えられたエネルギーを意味しています。
ここで、
\(d’Q\)は気体が受け取った熱量(熱的なエネルギーの授受)、
\(d’W\)は気体が受け取った仕事(力学的なエネルギーの授受)、
\(d’G\)は気体が受け取った、粒子数を変化させるのに必要なエネルギー(粒子数変化に伴うエネルギーの授受)
を意味します。ここで、これらの量が表す主語は気体です。気体が熱を受け取ったならば\(d’Q\gt 0\)であり、気体が仕事を受け取ったならば\(d’W\gt 0\)です。
また、プライムの記号は気体へのエネルギーの与え方(経路)を指定していることを意味し、微小量だけれども経路に依存する量を示すための記号です。
もう一度言いますが、この式は単なるエネルギー保存則を表しているだけです。

\(d’G\)とはどんなエネルギーなのか考えてみましょう。
例えば箱の中に入っている気体の粒子が、原子間力や分子間力などにより斥力相互作用をもつ単一の粒子で構成されていると考えましょう(※1)。
 この場合、新しい粒子を入れるためにはエネルギーが必要です。この過程が断熱的で体積変化が無いと想定すると、熱エネルギーでも力学的なエネルギーでもない別のエネルギーによって粒子が押し込められたということですよね。このエネルギーを化学ポテンシャルと呼びます。
 なので、粒子間相互作用が全くなければ、このエネルギーは生じません。

さて、状態変化を物理的、数学的に取り扱い易くするため、状態変化に制限を加えます。
それは準静的な変化と呼ばれる過程です。準静的な変化とは簡単に言えば可逆過程のみを考えるということです。
例えば気体を膨張する操作において、ピストンの摩擦力などによって気体以外の場所にエネルギーが渡ることが無い、ということです。
別の言い方をすれば、状態変化の途中では常に平衡状態が保たれているということを意味します。しかし常に平衡状態を保ちながら変化する、という言葉は矛盾しています。なぜなら、平衡状態とは変化が無い状態の事であり、それが変化するとはおかしい事です。平衡状態を保ちながら変化するということは、無限大の時間を掛けて無限小の変化を引き起こすということです。本当に無限ではなくて、気体にとって無限大の時間をかけて無限小の変化を起こすという意味です。

準静的な過程を考えた時、与えた熱などのエネルギーは気体の持つ適当な状態量の変化として現れます。先人方の研究によって、微小な熱量、微小な力学的エネルギー、微小な粒子数変化に伴うエネルギーは、
\begin{align}
d’Q&= TdS \
d’W&= -pdV \
d’G&= \sum_i \mu_i dN_i
\end{align}
と書くことが出来ることが分かりました。
ここで、
\(T\mathrm{[K]}\)は温度、\(S\mathrm{[J/K]}\)はエントロピー、\(P\mathrm{[N/m^2]}\)は圧力、\(V\mathrm{[m^3]}\)は体積、\(\mu\mathrm{[J/個数]}\)は化学ポテンシャル、\(N\mathrm{[個]}\)は粒子数を表します。
また、添え字\(i\)は気体を構成する原子・分子種の違いを表します。気体が酸素\(O_2\), 窒素\(N_2\)の2種から出来ているのであれば、\(i=1,2\)となります。

もう少し先でエントロピーの説明をしますが、本稿ではエントロピーの定義は単に、温度に共役な量と定義してそれ以上は余り述べません。

ここで、気体のエネルギーとは何をするためのエネルギーなのか触れておきましょう。
気体の状態量があれば、それに対応したエネルギーがあります。気体を特徴づける状態量から次元解析によりエネルギーの次元を作ろうと思えば、
(温度\(T\mathrm{[K]}\))×(エントロピー\(S\mathrm{[J/K]}\))はエネルギーの次元、
(圧力\(P\mathrm{[N/m^2]}\))×(体積\(V\mathrm{[m^3]}\))はエネルギーの次元、
(化学ポテンシャル\(\mu\mathrm{[J/個数]}\))×(粒子数\(\mathrm{[J]}\))はエネルギーの次元を持つことが分かります。
それぞれのエネルギーに対する微小変化に対しては、線形結合で書くことが出来ます。
すなわち、気体のなんらかのエネルギーの微小変化\(dE\mathrm{[J]}\)は、
\begin{align}
dE=c_1 d(TS) + c_2 d(PV) + c_3 d(\mu N)
\end{align}
と書かれることが予想できるでしょう。ここで\(c_1, c_2, c_3\)は各状態量の定義次第で決まる定数です。

上式が意味することは、3個の自由度を決めることで気体のエネルギーを定義できる、ということを意味しています。独立な変数として取れると言っても良いでしょう。
3個の自由度とは、それぞれのエネルギーに対応する自由度で、

\(T, S\)の中から1つ、
\(P, V\)の中から1つ、
\(\mu, N\)の中から1つ、

という意味です。つまり気体のエネルギーとして定義可能な組み合わせは8種類存在するということです。
8個のエネルギーを表現する形式のうち、どれが正しいとかそういうことは無くて、どれが便利な表現か?という基準で選びます。
問題によって、最適なエネルギーの表現があるということです。
もし上で表現されていないエネルギーの種類、例えば電気的なエネルギーが加わればそれに応じで種類も増えます。
ちなみに\(T, S\)や\(P, V\)や\(\mu, N\)の関係は、片方の変数はもう片方の共役な変数と呼ばれます。エネルギーの次元を持つ量を作る二つの変数の組み合わせを共役な変数と呼びます(※2)。例えば、\(T\)は\(S\)の共役な変数である、と言います。

平衡状態と最適なエネルギー


さて、3個の自由度を決めることで気体の持つ何らかのエネルギーを定義できると述べました。
つまり気体の自由度はエネルギーの表現が可能な、共役な変数の組の数できまり、上記の気体を表現したい場合ならば気体が持つ自由度が3であることを言っています。

続いて気体の平衡状態とは何を指す言葉なのか考えましょう。
平衡状態とは気体に変化を与えた後、十分長い時間放置した時に気体の状態量に変化が現れなくなった時のことを指します。平衡状態は3個の自由度を適切に決めた時に初めて定義できます。最適なエネルギーは問題設定から自ずと分かるでしょう。少し計算をしてから平衡状態を表すのに最適なエネルギーを述べていきます。

方針は、気体のエネルギーの代表例である内部エネルギー\(U\)について述べます。その後、\(U\)を基準にして、その他のエネルギー表記をルジャンドル変換を利用して導出していきます。

  • 熱のやりとり(\(d’Q=TdS\))と体積変化(\(dV\))を見たい時

    まず、最も基本的なのが気体の内部エネルギー\(U\)と呼ばれている量です。
    \(U\)は独立変数として\(S, V\)を選んだ場合の気体のエネルギーを指します。すなわち、
    \begin{align}
    dU=TdS -p dV
    \end{align}
    と定義されます。このように書いた時、\(U, T, p\)は\(S, V\)の関数として\(U(S, V), T(S, V), p(S, V)\)と書くことを示すことを注意してください。
    内部エネルギーの表記が有用なのは、
    熱のやりとりと気体の体積変化を見たい時
    です。熱のやり取りは\(d’Q=TdS\)、体積変化は\(dV\)で表されており、その変化に注目したい時便利な表記です。

  • 熱のやりとり(\(d’Q=TdS\))と圧力変化(\(dp\))を見たい時

    独立変数として\(S, p\)を選んだ場合の気体のエネルギーを考察しましょう。
    今、内部エネルギー\(U\)の微小変化は
    \begin{align}
    dU=TdS -p dV
    \end{align}
    です。目的は右辺を\(dS,dp\)で表すことです。共役な変数を取り換えるために、ルジャンドル変換を利用します。その方法は共役な変数の組み合わせを足したり引いたりすることです。便宜的にあるエネルギー\(H\)を考えて、それを
    \begin{align}
    H=U +pV
    \end{align}
    と定義します。このように新たな変数を考えるのがルジャンドル変換です。\(H\)の微小量を考えると
    \begin{align}
    dH&=d(U +pV) \\
    &=dU + pdV+Vdp
    \end{align}
    です。\(dU=TdS -p dV\)を代入すれば、
    \begin{align}
    dH=TdS+Vdp \\
    \end{align}
    を得ます。右辺には熱のやり取り(\(d’Q=TdS\))と圧力変化(\(dp\))で表されており、その変化に注目したい時便利な表記です。\(H\)はエンタルピーと呼ばれています。

  • 温度変化(\(dT\))と体積変化(\(dV\))を見たい時

    独立変数として\(T,V\)を選んだ場合の気体のエネルギーを考察しましょう。
    今、内部エネルギー\(U\)の微小変化は
    \begin{align}
    dU=TdS -p dV
    \end{align}
    です。目的は右辺を\(dT,dV\)で表すことです。ルジャンドル変換を利用して、
    便宜的にあるエネルギー\(F\)を考えて、それを
    \begin{align}
    F=U -TS
    \end{align}
    と定義します。\(F\)の微小量を考えると
    \begin{align}
    dF&=d(U -TS) \\
    &=dU – SdT-TdS
    \end{align}
    です。\(dU=TdS -p dV\)を代入すれば、
    \begin{align}
    dF=-SdT-pdV \\
    \end{align}
    を得ます。右辺には温度変化(\(dT\))と体積変化(\(dV\))で表されており、その変化に注目したい時便利な表記です。\(F\)はヘルムホルツの自由エネルギーと呼ばれています。

  • 温度変化(\(dT\))と圧力変化(\(dp\))を見たい時

    独立変数として\(T,p\)を選んだ場合の気体のエネルギーを考察しましょう。
    今、内部エネルギー\(U\)の微小変化は
    \begin{align}
    dU=TdS -p dV
    \end{align}
    です。目的は右辺を\(dT,dp\)で表すことです。ルジャンドル変換を利用して、
    便宜的にあるエネルギー\(G\)を考えて、それを
    \begin{align}
    G=U-TS+pV
    \end{align}
    と定義します。\(G\)の微小量を考えると
    \begin{align}
    dG&=d(U -TS+pV) \\
    &=dU – SdT-TdS +pdV+Vdp
    \end{align}
    です。\(dU=TdS -p dV\)を代入すれば、
    \begin{align}
    dG=-SdT+Vdp \\
    \end{align}
    を得ます。右辺には温度変化(\(dT\))と圧力変化(\(dV\))で表されており、その変化に注目したい時便利な表記です。\(G\)はギブスの自由エネルギーと呼ばれています。

以上の話でいろいろなエネルギー名が出てきましたが、それぞれの名称を覚える必要はありません。重要なのは考え方であり、迷ったら確かめれば良いのです。どうせ\(U\)に\(pV, TS\)を足したり引いたりすれば良いのですから。

一般的に書けば、熱力学における便利なエネルギー\(E\)は、
\begin{align}
E=U+c_1 TS + c_2 pV
\end{align}
の形で書けるわけです。微小量を考えれば、
\begin{align}
dE&=dU + c_1 TdS + c_1 SdT + c_2 pdV + c_2 Vdp \
&=TdS-pdV+ c_1 TdS + c_1 SdT + c_2 pdV + c_2 Vdp \
&=c_1 SdT + (c_1 +1)TdS+ (c_2-1)pdV + c_2 Vdp
\end{align}
と書けるわけです。状況に応じて、2つの微小量を消すようにそれぞれの係数\(c_1, c_2\)を選べばよい訳です。
そうすれば最適なエネルギー表記を見つけることが出来ます。

熱力学変数の各種関係式


状態数の間には関係式が成立します。
内部エネルギー\(U\)を用いると、
\begin{align}
dU=TdS -p dV
\end{align}
と書けます。\(U, T, p\)は\(S, V\)の関数として\(U(S, V), T(S, V), p(S, V)\)で書けると述べました。
今、\(U(S, V)\)だけに注目します。\(U(S, V)\)を単なる二変数関数として見れば、
\begin{align}
dU=\left(\frac{\partial U}{\partial S}\right)_V dS+ \left(\frac{\partial U}{\partial V}\right)_S dV
\end{align}
と書くことが出来ます。
ここで、\(\left(\frac{\partial U}{\partial S}\right)_V \)は、\(V\)を固定して\(U\)をSで偏微分する、という意味です。
わざわざ偏微分に()と添え字を付けている理由は、取り得る変数の値が複数あるからです。意味は、
\(\left(\frac{\partial U}{\partial S}\right)_V \)は、\(V\)を固定して\(U\)をSで偏微分する、という意味です。
例えば\(\frac{\partial U}{\partial S} \)とだけ書いた時、もう片方の変数が何なのか指定が無いので何でもよい、ということになってしまいますが、\(U=U(S,V)\)と\(U=U(S,p)\)で偏微分の値が変わります。その誤解を招かないために括弧と添え字で表現しています。

話を戻すと、今\(S, V\)を独立変数として書くとき、
\begin{align}
dU&=TdS -p dV \
dU&=\left(\frac{\partial U}{\partial S}\right)_V dS+ \left(\frac{\partial U}{\partial V}\right)_S dV
\end{align}
と二通りの表現方法が今あるわけです。となると、係数同士を比較すると、
\begin{align}
\left(\frac{\partial U}{\partial S}\right)_V=T,~ \left(\frac{\partial U}{\partial V}\right)_S=-p \
\end{align}
という関係式が成立していることが分かります。

また、二変数関数\(f(x,y)\)の二階偏微分には
\begin{align}
\frac{\partial}{\partial y}\frac{\partial f}{\partial x} = \frac{\partial}{\partial x}\frac{\partial f}{\partial y}
\end{align}
という関係式があります。従って、
\begin{align}
\left(\frac{\partial }{\partial V}\left(\frac{\partial U}{\partial S}\right)_V\right)_S=
\left(\frac{\partial }{\partial S}\left(\frac{\partial U}{\partial V}\right)_S\right)_V
\end{align}
が成立しているので、全微分と熱力学第一法則の比較から得た関係式を代入して
\begin{align}
\left(\frac{\partial T}{\partial V}\right)_S=-\left(\frac{\partial p}{\partial S} \right)_V
\end{align}
が成立していることが分かります。このように二階偏微分を利用して比較して得た関係式はマクスウェルの関係式と呼ばれています。

マクスウェルの関係式を見て分かる通り、これはエネルギーを含まないので純粋な状態量同士の関係です。
同様の関係式は他のエネルギーを考えた結果でも分かることで、\(S, T\)と\(p, V\)の組み合わせから合計4つの便利なエネルギー表式(内部エネルギー、エンタルピー、ヘルムホルツの自由エネルギー、ギブスの自由エネルギー)から4つの関係式を導くことが出来ます。

熱容量


気体に熱量\(Q\)を加えた時、気体の温度\(T\)が何度上がるか考えてみましょう。
見るべき変化量は、気体のエネルギー\(TdS\)であり、その温度変化を考える事です。知りたい量を数式で表せば、
\begin{align}
\frac{d’Q}{dT}=\frac{TdS}{dT}
\end{align}
という値です。この量をエネルギーの関数として表すことが目標です。
まずはこの問題を解くにあたり\(S, T\)と\(p, V\)のどの組み合わせが良いか選びます。

今、\(d’Q\)を陽に表したいので、最適なエネルギーを選ぶ基準として\(d’Q=TdS\)が含まれるエネルギーの形式を選びます。すなわち、\(S\)を変数の1つとして選びます。残りは\(p, V\)の内から選びます。この自由度があることから、\(\frac{d’Q}{dT}\)は一意には決まりません。どちらかの変数を選んだ時の熱容量が定義できるだけです(※3)。

定積熱容量


まずは体積が変わらない時の熱容量、すなわち定積熱容量\(C_V\)を計算しましょう。数式で求めたい量を書けば、
\begin{align}
C_V=\left(\frac{TdS}{dT}\right)_V
\end{align}
という量です。右辺を見ると、

  • \(S\)の\(T\)による変化量を考えたい
  • \(V\)は変化しない

という事に注目したいので、\(S=S(T,V)\)と書けていれば都合が良さそう、ということが分かります。
すなわち、\(S(T,V)\)について
\begin{align}
dS(T,V)=\left(\frac{\partial S(T,V)}{\partial T}\right)_V dT + \left(\frac{\partial S(T,V)}{\partial V}\right)_T dV
\end{align}
の関係式を用いることが出来ます。この式の意味は、単なる全微分という意味だけではなく、変数\(S\)を\(T,V\)に変換する式とみることが出来ます。

左辺の\(dS(T,V)\)を気体の持つエネルギーの表現に代入したいと考えましょう。すなわち\(S, p\)の変数を持つ内部エネルギーが良いです(※4)。
\begin{align}
dU(S,V)=T(S,V)dS -p(S,V)dV
\end{align}
より、\(dS\)に代入することで\(S,V\)から\(T,V\)に変数変換して
\begin{align}
dU(S(T,V),V)&=T(S(T,V),V)\left[\left(\frac{\partial S(T,V)}{\partial T}\right)_V dT + \left(\frac{\partial S(T,V)}{\partial V}\right)_T dV \right] -p(S(T,V),V)dV \
&=T(S(T,V),V)\left(\frac{\partial S(T,V)}{\partial T}\right)_V dT
+ \left[\left(\frac{\partial S(T,V)}{\partial V}\right)_T -p(S(T,V),V) \right]dV \
&=T\left(\frac{\partial S(T,V)}{\partial T}\right)_V dT
+ \left[\left(\frac{\partial S(T,V)}{\partial V}\right)_T -p(T,V) \right]dV \
&=C_V dT
+ \left[\left(\frac{\partial S(T,V)}{\partial V}\right)_T -p(T,V) \right]dV
\end{align}

と変形できます。内部エネルギーを全微分を利用して\(T, V\)で表現すれば、
\begin{align}
dU(T,V)=\left(\frac{\partial U(T,V)}{\partial T}\right)_V dT + \left(\frac{\partial U(T,V)}{\partial V}\right)_T dV
\end{align}
なので、\(C_V\)は
\begin{align}
C_V=\left(\frac{\partial U(T,V)}{\partial T}\right)_V
\end{align}
と導けます。すなわち定積熱容量は、体積が変わらない時の内部エネルギーの温度変化に等しい、ということです。ここまで分かりましたが内部エネルギーの明確な格好が全く分かっていないのでこれ以上の計算は出来ないままです。

定圧熱容量


続いて圧力が変わらない時の熱容量、すなわち定圧熱容量\(C_p\)を計算しましょう。数式で求めたい量を書けば、
\begin{align}
C_p=\left(\frac{TdS}{dT}\right)_p
\end{align}
という量です。同様に考えていくと、

  • \(S\)の\(T\)による変化量を考えたい
  • \(p\)は変化しない

という事に注目したいので、\(S=S(T,p)\)と書けていれば都合が良さそう、ということが分かります。
すなわち、\(S(T,p)\)について
\begin{align}
dS(T,p)=\left(\frac{\partial S(T,p)}{\partial T}\right)_p dT + \left(\frac{\partial S(T,p)}{\partial p}\right)_T dp
\end{align}
の関係式を用いることが出来ます。左辺の\(dS(T,p)\)を気体の持つエネルギーの表現に代入したいと考えましょう。定積熱容量の時と同じように、\(dS\)を消去するという方針で行きます。\(S,p\)を変数に取る最適な気体のエネルギーはエンタルピー\(H\)ですので、
\begin{align}
dH(S,p)=T(S,p)dS +V(S,p)dp
\end{align}
と書けます。\(dS\)に代入することで\(S,p\)から\(T,p\)に変数変換して
\begin{align}
dU(S(T,p),p)&=T(S(T,p),p)\left[\left(\frac{\partial S(T,p)}{\partial T}\right)_p dT + \left(\frac{\partial S(T,p)}{\partial p}\right)_T dp \right] -p(S(T,p),p)dp \
&=T(S(T,p),p)\left(\frac{\partial S(T,p)}{\partial T}\right)_p dT
+ \left[\left(\frac{\partial S(T,p)}{\partial p}\right)_T -p(S(T,p),p) \right]dp \
&=T\left(\frac{\partial S(T,p)}{\partial T}\right)_p dT
+ \left[\left(\frac{\partial S(T,p)}{\partial p}\right)_T -p(T,p) \right]dp \
&=C_p dT
+ \left[\left(\frac{\partial S(T,p)}{\partial V}\right)_T -p \right]dp
\end{align}

と変形できます。エンタルピーを全微分を利用して\(T, p\)で表現すれば、
\begin{align}
dU(T,p)=\left(\frac{\partial H(T,p)}{\partial T}\right)_p dT + \left(\frac{\partial H(T,p)}{\partial p}\right)_T dV
\end{align}
なので、\(C_p\)は
\begin{align}
C_p=\left(\frac{\partial H(T,p)}{\partial T}\right)_p
\end{align}
と導けます。すなわち定圧熱容量は、体積が変わらない時のエンタルピーの温度変化に等しい、ということです。ここまで分かりましたが定積熱容量の時と同様にエンタルピーの明確な格好が全く分かっていないのでこれ以上の計算は出来ないままです。

定積熱容量と定圧熱容量の関係


定積熱容量と定圧熱容量の関係を考えてみましょう。現在までにそれぞれ
\begin{align}
C_V=\left(\frac{\partial U(T,V)}{\partial T}\right)_V,~C_p=\left(\frac{\partial H(T,p)}{\partial T}\right)_p
\end{align}
であることが導かれています。\(C_V, C_p\)は別のエネルギーで議論されているので同じエネルギーで議論したいと考えます。\(U\)に統一するとすれば\(C_p\)は、
\begin{align}
C_p&=\left(\frac{\partial H(T,p)}{\partial T}\right)_p \
&=\left(\frac{\partial (U(T,p)+p V(T,p))}{\partial T}\right)_p \
&=\left(\frac{\partial U(T,p)}{\partial T}\right)_p
+p\left(\frac{\partial V(T,p)}{\partial T}\right)_p \
&=\left(\frac{\partial U(T,p)}{\partial T}\right)_p+p\left(\frac{\partial V(T,p)}{\partial T}\right)_p
\end{align}
と変形できます。すると、\(U\)を含む項があるので
\begin{align}
\left(\frac{\partial U(T,V)}{\partial T}\right)_V,~~\left(\frac{\partial U(T,p)}{\partial T}\right)_p
\end{align}
の関係が導ければよいと分かります。とりあえず内部エネルギー\(U\)の全微分をそれぞれの変数で考えてみましょう。
\begin{align}
dU(T,V)&=\left(\frac{\partial U(T,V)}{\partial T}\right)_V dT+\left(\frac{\partial U(T,V)}{\partial V}\right)_T dV \
dU(T,p)&=\left(\frac{\partial U(T,p)}{\partial T}\right)_p dT+\left(\frac{\partial U(T,p)}{\partial p}\right)_T dp
\end{align}
です。それぞれの微小量で比較することが出来れば良い訳です。

\(dp\)の消去

\(dp=dp(V,T)\)であるので、\(dp\)の中に温度依存が含まれていることになります。そこで\(dp\)を\(dV,dT\)として表す変数変換を行います。すなわち
\begin{align}
dp(T,V)=\left(\frac{\partial p(T,V)}{\partial T}\right)V dT+\left(\frac{\partial p(T,V)}{\partial V}\right)_T dV
\end{align}
を代入して
\begin{align}
dU(T,p(T,V))&=\left(\frac{\partial U(T,p(T,V))}{\partial T}\right)
{p(T,V)} dT+\left(\frac{\partial U(T,p(T,V))}{\partial p(T,V)}\right)T dp(T,V) \
dU(T,V)&=\left(\frac{\partial U(T,V)}{\partial T}\right)
{p(T,V)} dT+\left(\frac{\partial U(T,V)}{\partial p(T,V)}\right)T \left[\left(\frac{\partial p(T,V)}{\partial T}\right)_V dT+\left(\frac{\partial p(T,V)}{\partial V}\right)_T dV\right] \
&=\left[ \left(\frac{\partial U(T,V)}{\partial T}\right)
{p(T,V)}+\left(\frac{\partial U(T,V)}{\partial p(T,V)}\right)T\left(\frac{\partial p(T,V)}{\partial T}\right)_V \right] dT +
\left[\left(\frac{\partial U(T,V)}{\partial p(T,V)}\right)_T\left(\frac{\partial p(T,V)}{\partial V}\right)_T\right] dV \
&=\left[ \left(\frac{\partial U}{\partial T}\right)
{p}+\left(\frac{\partial U}{\partial p}\right)T\left(\frac{\partial p}{\partial T}\right)_V \right] dT +
\left[\left(\frac{\partial U}{\partial p}\right)_T\left(\frac{\partial p}{\partial V}\right)_T\right] dV
\end{align}
となります。最後の式ではすべてが\(T, V\)の関数として表現されています。そのため、\(dT, dV\)に掛かる係数同士を比較して
\begin{align}
\left(\frac{\partial U}{\partial T}\right)_V&=\left(\frac{\partial U}{\partial T}\right)
{p}+\left(\frac{\partial U}{\partial p}\right)_T\left(\frac{\partial p}{\partial T}\right)_V \
\left(\frac{\partial U}{\partial V}\right)_T&=\left(\frac{\partial U}{\partial p}\right)_T\left(\frac{\partial p}{\partial V}\right)_T
\end{align}
となります。2番目の式の右辺は、結局左辺と同じになるので意味はない式です。1番目の式は内部エネルギーの定圧変化、定積変化の関係を述べています。
よって、定圧熱容量は
\begin{align}
C_p&=\left(\frac{\partial U}{\partial T}\right)_p+p\left(\frac{\partial V(T,p)}{\partial T}\right)_p \
&=\left(\frac{\partial U(T,V)}{\partial T}\right)_V-\left(\frac{\partial U}{\partial p}\right)_T\left(\frac{\partial p}{\partial T}\right)_V +p\left(\frac{\partial V}{\partial T}\right)_p\
&=C_V-\left(\frac{\partial U}{\partial p}\right)_T\left(\frac{\partial p}{\partial T}\right)_V +p\left(\frac{\partial V}{\partial T}\right)_p
\end{align}
となることが導かれます。

\(dV\)の消去

\(dV=dV(T,p)\)であるので、\(dV\)の中に温度依存が含まれていることになります。そこで\(dV\)を\(dp,dT\)として表す変数変換を行います。すなわち
\begin{align}
dV(T,p)=\left(\frac{\partial V(T,p)}{\partial T}\right)p dT+\left(\frac{\partial V(T,p)}{\partial p}\right)_T dp
\end{align}
を代入して
\begin{align}
dU(T,V(T,p))&=\left(\frac{\partial U(T,V(T,p))}{\partial T}\right)
{V(T,p)} dT+\left(\frac{\partial U(T,V(T,p))}{\partial V(T,p)}\right)T dV(T,p) \
dU(T,p)&=\left(\frac{\partial U(T,p)}{\partial T}\right)
{V} dT
+\left(\frac{\partial U(T,p)}{\partial V}\right)T
\left[
\left(\frac{\partial V(T,p)}{\partial T}\right)_p dT+\left(\frac{\partial V(T,p)}{\partial p}\right)_T dp
\right] \
&=\left[\left(\frac{\partial U}{\partial T}\right)
{V}+\left(\frac{\partial U}{\partial V}\right)T\left(\frac{\partial V}{\partial T}\right)_p\right] dT
+\left(\frac{\partial U}{\partial V}\right)_T
\left(\frac{\partial V}{\partial p}\right)_T dp
\end{align}
となります。最後の式ではすべてが\(T, p\)の関数として表現されています。そのため、\(dT, dp\)に掛かる係数同士を比較して
\begin{align}
\left(\frac{\partial U}{\partial T}\right)_p&=\left(\frac{\partial U}{\partial T}\right)
{V}+\left(\frac{\partial U}{\partial V}\right)T\left(\frac{\partial V}{\partial T}\right)_p \
\left(\frac{\partial U}{\partial p}\right)_T&=\left(\frac{\partial U}{\partial V}\right)_T\left(\frac{\partial V}{\partial p}\right)_T
\end{align}
となります。2番目の式の右辺は、結局左辺と同じになるので意味はない式です。1番目の式は内部エネルギーの定圧変化、定積変化の関係を述べています。よって、定圧熱容量は
\begin{align}
C_p&=\left(\frac{\partial U}{\partial T}\right)_p+p\left(\frac{\partial V(T,p)}{\partial T}\right)_p \
&=\left(\frac{\partial U}{\partial T}\right)
{V}+\left(\frac{\partial U}{\partial V}\right)_T\left(\frac{\partial V}{\partial T}\right)_p+p\left(\frac{\partial V}{\partial T}\right)_p\
&=C_V+
\left[p+\left(\frac{\partial U}{\partial V}\right)_T\right]\left(\frac{\partial V}{\partial T}\right)_p
\end{align}
となることが導かれます。

\(dV\)と\(dp\)を消去した形式、どちらでも良いですが後者を採用します。少し変形して
\begin{align}
C_p-C_V=
\left[p+\left(\frac{\partial U}{\partial V}\right)_T\right]\left(\frac{\partial V}{\partial T}\right)_p
\end{align}
と書けます。

ここまでで熱容量の話は終わりです。理想気体の場合、マイヤーの関係式と呼ばれる関係を書くことが出来ますが、混乱するので本稿では載せません。

エントロピーについて


エントロピーについての詳細は統計力学を元に考えたほうが分かりやすいのでここでは詳しくは書きません。
気体は原子・分子の集団です。エントロピーは気体が取り得る微視的な状態数(量子力学の固有状態数)によって決まる量です。

例えば気体を構成する原子数と、気体の温度、気体の持つ総エネルギーが与えられたとき、それらの条件を満たすような原子の組み合わせ数、という感じです。
具体例を挙げましょう。原子数が5個、総エネルギーが5Jでどのような分配の仕方が考えられるか?という問題です。イメージだけつかむことを目標にするので、温度は無視します。量子力学では1つの原子が持つエネルギーは離散的にしか存在できないのでここでは1J刻みにしか取れないと仮定します。この時、

a | 5 0 0 0 0
b | 1 1 1 1 1
c | 0 0 0 1 4
d | 0 0 0 0 5

といったような様々な組み合わせが考えられます。数字は原子1~5が持つエネルギーを意味します。エントロピーとは、この組み合わせの数の総数に関係する数です。組み合わせ数が多ければ多いほどエントロピーが大きいと言いますし、組み合わせ数が少なければエントロピーが小さいと表現します。

断熱系では、気体のエントロピーの変化はありません。\(d’Q=TdS\)であり、\(T\gt 0\)ですから、断熱的\(d’Q=0\)ならば\(dS=0\)となります。これが意味するのは、熱が加わらない準静的な過程ではエントロピーが保存量になっているということです。

最後に


問題をみてすぐさま変数が思い浮かぶことは無いので、練習を重ねたりトライアンドエラーを重ねましょう。上の例では一本道のように変数を決定していますが、私も別の変数を選んでしまいうまく表現できなくで何度も失敗しています。たまたまうまくいったのをあたかも最初から分かったように説明しているだけにすぎません。

注釈


※1
 イメージをする上では、気体を構成する粒子ひとつひとつが、すべて電気的に正に帯電していると考えても良いでしょう。正に帯電しているのだから、新しい粒子を入れるためにはエネルギーが必要です。この過程が断熱的で、体積変化が無いと想定すると、熱エネルギーでも力学的なエネルギーでもない別のエネルギーによって粒子が押し込められたということですよね。このエネルギーが化学ポテンシャルなのです。
 詳しく言うと、熱力学第一法則で言われる気体は中性気体が想定されており、この場合に限り化学ポテンシャルと呼ぶようです。電気的に帯電していることによるエネルギー変化は、電気化学ポテンシャルと呼ばれています。
 また、化学ポテンシャルの由来は何から来るか?という問いの答えは場合による、です。場合によっては例に示したように分子間力であるし、原子であれば分極が関係しているのかもしれません。場合によっては電気化学ポテンシャルも化学ポテンシャルと読んでしまう場合もあるでしょう。
 なので、粒子の流入に関係する様々な要因をひっくるめて化学ポテンシャルと呼んでいるようです。

※2
 「共役な」という意味は熱力学の世界ではエネルギーの次元を作る二つの変数の関係のことを言いますが、分野が異なると別の物を作る二つの変数の関係を指す言葉となります。複素解析では複素共役という言葉が用いられ、ある複素数とその共役な複素数を掛けた時に大きさになるような複素数の関係を指します。古典力学や量子力学では、位置と運動量の関係や、時間と周波数の関係のような関係を言ったりします。

※3
 \(p, V\)の中から一つを選択するのではなく、例えば新たな変数として\(p+cV, p-cV\)を選んだ時の熱容量を定義することも可能でしょう(\(c\)は定数)。しかし\(p+cV, p-cV\)は人間が分かりやすい量ではないので通常は用いません。また、\(p, V\)を選んでいっても、適当な操作で\(p+cV, p-cV\)の張る空間全てに辿り着くので、\(p, V\)のみに対する熱容量を考えるだけで良いです。

※4
定積熱容量なので\(T,V\)の変数を用いているヘルムホルツの自由エネルギー\(F(T,V)\)を元に議論を進めていくのが良いと思いますが、うまい具合に\(S\)について解けないのです。なので、\(dS\)を消去するという考えで行うと都合よく比較ができるのでそのように行っています。

参考文献


伊藤敏雄著『な~るほど!の力学』学術図書出版社(第1版11刷2009年, 初版1995年)
小出昭一郎著『物理学(改訂版)』裳華房(第25版1984年, 初版1975年)

動かない壁に対する束縛運動と反射

動かない壁に対する束縛運動と反射を考えます。
例えば、初め跳ねてた運動が、壁に沿って動く運動に変化する、という状況です。
あんまり見たことが無いので、面白そうだと思いました。

束縛されている時と反発するときは、式(1),(2)によってあらわされます。それらは
壁に束縛されている場合

壁と反発する場合

です。壁と反発する場合、反発後の速度は式(3)に沿って動きます。

※式(1), (2)では壁が時間依存しても良いように定式化しています。この定式化は恐らく正しいです。また、本稿に載せているプログラムも壁が時間依存しても良いように作成していますが、動く壁の場合、プログラムではうまく計算が出来ません。

ここで、\(e\)は反発係数、\(\mathbf{n}\)は壁の法線ベクトルであり、

\(\nabla\)はナブラ演算子であり

と与えられます。また、\(\hat{H}f\)は関数\(f\)のヘッセ行列であり,

と与えられます。

定式化や数値計算手法の詳細は以下のページを参照してください。

壁との反発と束縛運動の定式化

質点と壁との反発を表す運動方程式
束縛条件下の運動 – 自由度がうまく落とせない運動

数値計算手法

ルンゲ=クッタ法の説明と刻み幅制御
Hyper-dual numbersによる二階偏微分の計算
ゼロ点を探す(二分法、挟み撃ち法、Anderson-Björk法、Brent法、Newton法、Steffensen法)


解放・束縛判定


ここでいう”解放”とは、束縛されていなく、反射を繰り返している状態を表します。また、”束縛”は壁に沿って運動している状態です。

前提として、壁を通り抜けることは無いと考えます。

すなわち、時刻\(t=t_0\)で位置\(\mathbf{r}=\mathbf{r}_0\)の時、任意の時刻\(t\)について

が満たされるとします。
もし、\(f(\mathbf{r}_0,t_0) =0\)ならば判断がつかないため、プラスかマイナスはこちらから与えます。

解放→解放判定、解法→束縛判定


解放状態から壁によって単に反発する場合

関数\(f(\mathbf{r},t)\)の符号が変化した時、壁との反発を考えます。
壁の法線方向の速度が十分に大きい場合、壁と反発し、そうでない場合、壁に束縛されると考えます。
その前提の元、

を満たす\(t=t_i, \mathbf{r}=\mathbf{r}_i\)を見つけます。
その後、式(6)に従い、速度を変化させます。

数値計算的には、衝突直前の時刻を採用します。すなわち、
ゼロ点を探す際にある範囲\(t_a \leq t_i \leq t_b\)で挟み込んでいくのですが、\(t_b\)は壁を越えてしまうので採用しません。

解放状態から壁沿いに束縛される場合

もし、壁の法線方向の速度が十分に小さい場合(ある閾値を下回った場合)、壁に束縛されると考えます。
この時、壁の法線方向の速度はゼロに変更します。
すなわち、速度は時刻\(t=t_i\)で

を持ちますが、\(v_{\parallel}=0\)にしてから、束縛運動に移行するということです。
これは、\(\mathbf{v}\)と壁の法線方向のベクトル\(\mathbf{e}_n\)の内積を取ることで得られます。
また、束縛された瞬間(\(t=t_i\))の束縛力\(C(t_i)\)を計算し、その符号を記録しておきます。
束縛力\(C(t)\)は、

です。
この符号が変化した瞬間が壁からの束縛が無くなる時なので、そのために記録します。

束縛→解放判定

束縛力\(C(t)\)の符号と\(C(t_i)\)の符号が変わるまで、式(1)に従い、時間発展させます。
すなわち、

を満たす\(t=t’, \mathbf{r}=\mathbf{r}’\)を見つけます。
その後、式(2)に従い運動します。
式(2)の運動では束縛力は働かないので、符号は自然と初期条件の符号と同一になる(※)。

※この条件はあまり良くありません。この判別方法のせいで、壁が時間依存している場合、束縛力が働いていない一瞬で質点が壁を超えてしまいます。プログラム自体は壁は時間依存しても良いことになっていますが、この条件分岐は上手く動きません。以下に示す計算プログラムは、質点が束縛されている場合に壁が時間依存しなければ正しいです。

プログラム


プログラムは以下のリンク先においておきます。

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

適当なディレクトリで展開し、lag_ver1.0というディレクトリに移動してから以下のコマンドで実行できます。

 $ sh cm.sh
 $ ./a.out
&INPUT
 MASS=  1.0000000000000000     ,
 G=  1.0000000000000000     ,
 TA=  0.0000000000000000     ,
 TB=  20.000000000000000     ,
 NT=        101,
 ELS= 0.59999999999999998     ,
 RX0= -1.0000000000000000     ,
 RY0=  1.0000000000000000     ,
 VX0=  3.0000000000000000     ,
 VY0=  0.0000000000000000     ,
 RKTOL=  1.0000000000000000E-008,
 ZRTOL=  1.0000000000000000E-008,
 TRTOL=  1.0000000000000000E-008,
 REGION=          1,
 /
           0
 $ gnuplot

    G N U P L O T
    Version 4.6 patchlevel 4    last modified 2013-10-02
    Build System: Linux x86_64

    Copyright (C) 1986-1993, 1998, 2004, 2007-2013
    Thomas Williams, Colin Kelley and many others

    gnuplot home:     http://www.gnuplot.info
    faq, bugs, etc:   type "help FAQ"
    immediate help:   type "help"  (plot window: hit 'h')

Terminal type set to 'wxt'
gnuplot> load "anime.plt"
gnuplot>

動かすと以下のような動画が描画されます。

デフォルトでは
ポテンシャル\(V=mgy\)(サブルーチンfp2d)
壁\(f(x,y)=x^2+y^2-4\)(サブルーチンfw2d)
に書かれています。
摩擦、空気抵抗は入っていません。唯一、反発係数(els)だけがinputファイルの中に書かれています。

初期条件の

rx0  = -2d0,  ! Initial condition
ry0  =  2d0,  !     position and velocity
vx0  =  1d0,  !  \mathbf{r} = (rx, ry)
vy0  =  0d0,  !  \mathbf{v} = (vx, vy)

だけを上の通り変更すると以下のような振る舞いをします。

確かめ


確かめを行います。
重力\(g\)の下で、質量\(m\)の物体が、半径\(r\)の円形の壁の内側に沿ってボールが進み、壁からの抗力が無くなり、壁から離れる状況を考えます(下の図を参照)。

円形の壁に沿っている時、垂直抗力\(N\)は

と書けます。エネルギー保存則より、

が成り立っています。ここで、\(v_0\equiv v(t=0)\)と置きました。
垂直抗力\(N\)がゼロになる点が壁から離れる条件ですので、\(v_0\)を用いて

と書けます。初速度が分かっている時、壁から離れる角度は

です。もしくは、壁から離れる角度が分かっている場合、初速度は

と与えられます。重力加速度, 半径を\(g=1, r=2\)とし、\(y=1\)と壁との交点、すなわち\(\theta=2\pi/3\)の場合、初速度は

です。実際に本稿のプログラムを動かしますと、\(y=1\)でちょうど離れていることが確認できます。

ここで、青線は壁に沿って動いて運動している状態であり、赤線は壁から離れている運動している状態です。

剛体に対する散乱

ポテンシャルを無くし(g=0)、弾性散乱(els=1)を考えると
古典的な散乱問題的なものもできます。

RLC並列回路の過渡現象

RLC並列回路を考えます。

良くあるページではインピーダンスを考えるだけでよし、としていますが、過渡現象が知りたいですよね。私はそうです。ですので、ラプラス変換を用いて解いていきましょう。

考える回路は以下の図の通りです。

回路方程式を立てれば、



なります。この式(1)と式(2)の計4本の連立方程式を解いて、未知の関数\(i(t),i_1(t),i_2(t),i_3(t)\)を求めることが目標です。

電圧、電流のラプラス変換を

と書くことにします。式(2)を式(1)に代入して、ラプラス変換を施せば、

を得ます。もっと厳密に書けば、式(1b)のラプラス変換には\(i_2(0^-)\)という項が含まれますが、\(i_2(0^-)=0\), すなわち電源がオンになるまでは電流は存在しないと仮定します。

行列表示にすれば

となります。関数\(F_k(s)\)にとって線形の問題です。表記を簡単にするために、

と書くことにします。ここで、\(x, y, z\)は

を意味します。ただの定数です。
クラメルの公式を用いて式(6)を解きます。関数

を定義すれば、式(6)の解は

と書くことが出来ます。あらわに\(G(x,y,z)\)を書けば

と表すことが出来ます。ここで、\(\alpha, \beta\)は

を満たすような解として書きました。
以降、\(\alpha, \beta\)は同じでは場合を考えていきます。
具体的に式(9)に現れる量を計算していけば、

となります。それぞれ、



という結果が得られます。ここにいたるまでに、



という関数を定義しました。
すると、それぞれの電流は

と書くことができ、具体的に

と求められます。
この結果を導くにあたって使用した仮定をまとめますと、
1. \(\alpha\ne \beta\)
2. \(i(t\lt 0)=0\)
という仮定の下、導き出された結果です。

導いた式(20)が合っているのか、直流電圧源を考えて考察してみましょう。
直流電圧源が時刻\(t=t_0\geq 0\)にスイッチオンする場合、電圧は

と書くことが出来ます。式(20)に出てくる積分は

と計算できますので、それぞれの素子を流れる電流は

と書きあらわすことが出来ます。
では、\(t\to\infty\)の漸近形を考えてみましょう。
\(t\to\infty\)の振る舞いは指数関数の型の\(t\)にかかる係数の実部の符号によって支配されます。式(11)より、具体的に\(\alpha, \beta\)を書くと

と求められます。今、\(R, C, L\)は全て正の実数ですので、ルートの項がその前の項の絶対値よりも大きくなることはありません。
よって、\(\text{Re}(\alpha)\lt 0, \text{Re}(\beta)\lt 0\)が導けます。すなわち、

という結果です。
これを式(23)に代入すれば、

という結果が得られます。
電源が入った後に系が十分に落ち着いた定常状態では、コンデンサーは開放、コイルは短絡されたものと見なしてよいので、その結果に見合った振る舞いであることが分かります。

メモとして書いておきますが、式(11)を計算するにあたって\(\alpha \beta\)を計算する必要があります。これは式(24)から求めるのではなく、式(11)の右辺を展開して、\(s^0\)の項を比較すれば

という結果が得られます。

ラプラス変換による回路方程式の解

  1. 回路方程式の組立
  2. ラプラス変換、畳み込み、超関数の積分について
  3. RC直列回路
  4. RL直列回路
  5. RLC直列回路

回路方程式の組立


閉じた回路があって、既知の電圧源が与えられたとき、回路に流れる電流を考えます。
色々考えた結果、間違えないようにするためにキルヒホッフの法則を信じて、
以下のように組み立てれば良いと思いました。

キルヒホッフの法則
閉じた回路内の起電力の和が、閉回路内の電圧降下の和に等しい
という法則を考えます。

ですが、起電力と電圧降下を別に考える、すなわち符号を取り換えるのは良く分かりません。
本稿では、上の起電力と電圧降下を区別しないことにします。
すなわち、キルヒホッフの法則を単純に
\(\displaystyle
\sum_k E_k =0
\)

と書くことにします。
素子の電圧降下分\(E_k\)はそれぞれ

で表されるとします。

電圧源の符号は、考えた閉経路の方向(電流の方向)に沿って電位が上がる場合は正、
反対方向に沿う場合は負として考えます。
素子の電圧降下は何も考えずに、上の表に対応する電圧降下分を書けば良いです。

(補足)回路方程式と運動方程式の対応

回路方程式とニュートンの運動方程式は物理的には全く異なりますが、数学的に似ています。すなわち

電荷\(q(t)\)は、質点の位置\(x(t)\)に対応し、
電流\(i(t)\)は、質点の速度\(v(t)\)に対応、
電圧\(v(t)\)(電位差)は、保存力\(F(t)\)に対応するとみることが出来ます。

ですが、キルヒホッフの法則は、古典力学で対応する法則がありません。
強いていうのであれば、元の位置に戻ってきたときに電位差がゼロなわけですから、
回路方程式は必ず保存力である、もしくはポテンシャルが必ず定義できなければならない、みたいなことでしょうか?
あんまりしっくりきません。

テレゲンの法則と呼ばれる法則が元にあると考えるのが良さそうです。これは、
\(\displaystyle
\sum_k E_k i_k =0
\)

と書かれるので、電圧と電流の積、すなわち単位時間当たりのエネルギー保存則に対応するわけです。
誤解を恐れず、単なるイメージ的に言えば、テレゲンの法則はハミルトニアン的な物として捉えることが出来るでしょう。

ラプラス変換、畳み込み、超関数の積分について


本稿では、ある時刻\(t\)から突然電源が入る、とかそういう現象を扱いたいので、
ラプラス変換を用いて回路方程式を解いていきます。

詳細は述べませんが、関数\(f(t)\)のラプラス変換\(F(s)\)を\(\mathcal{L}[f(t)](s)\)、逆ラプラス変換は\(\mathcal{L}^{-1}[F(s)](t)\)と書くと、

と定義されます。ラプラス変換は変換表を用いて計算することが殆どなので、詳細について知らなくても良いでしょう。
別の言い方をすれば、初等的な範囲でラプラス変換表に乗ってない関数場合、数値的に解くしかないと考えても良いです。

ラプラス変換表はLaplace transform -wikipedia や、wolfram のlaplace transform of exp(-a*t) -wolfram alphaで調べれば良いでしょう。逆変換もinverseとか付け加えれば可能です。

関数同士の掛け算の逆ラプラス変換は畳み込みと関係しています。なので、畳み込みも書いておきましょう。

また、ヘヴィサイド関数\(\theta(t)\)とディラックのデルタ関数\(\delta(t)\)を含んだ積分も頻出しますので、書いておきます。必要になったら参照してください。

RC直列回路


キルヒホッフの法則から出発して、RC直列回路を考えます。

回路方程式は

と書けます。この方程式を解いて、電流\(i(t)\)を求める事が目標です。

辺々にラプラス変換を施せば、

を得ます。ここで、

と置きました。式(2)を\(F(s)\)について解けば、

と書くことが出来ます。ここで、\(f\ast g(x)=\int_{0^-}^t f(t-\tau)g(\tau)d\tau\)は関数\(f(t)\)と\(g(t)\)の畳み込みを表し、さらに

を定義しました。式(4c)に現れる\(w(t)\)を求めるには、\(W(s)\)の逆ラプラス変換を行えばよいです。計算すれば、

となります。ここで、\(\delta(x)\)はディラックのデルタ関数、\(\theta(x)\)はヘヴィサイド関数です。

式(4d)の両辺を逆ラプラス変換を施すことで、電流\(i(t)\)を得ることが出来ます。計算すれば、

となり、求めたかった電流のあらわな式(7e)が得られます。
ここに至るまで、電圧源が直流であるとか、交流であるとかそういう条件は用いていません。

直流電圧源の場合


直流電圧源を考えます。状況としては、時刻\(t=t_0\)にスイッチがオンになり、電圧\(E_0\)の定電圧が掛かり、電流が流れ始めるという状況を考えます。
すると、電圧源\(v(t)\)は
\(\displaystyle
v(t)=E_0\theta(t-t_0)
\)
と書くことが出来ます。実際に式(7e)に代入して、

となります。もしも\(t_0>0\)ならば、

という結果を得ます。抵抗\(R\), キャパシタンス\(C\)は正ですので、\(t\to\infty\)の漸近で電流は流れなくなることが分かります。

RL直列回路


続いて、RL直列回路を考えます。

回路方程式は

と書くことが出来ます。式(10)の両辺にラプラス変換を施せば

を得ます。\(F(s)\)は\(i(t)\)のラプラス変換を意味しています。\(F(s)\)について解き、RC回路の時と同様に\(i(t)\)を求めれば、

と書くことが出来ます。すなわち電流\(i(t)\)は式(12c)に逆ラプラス変換を施して

と得ることが出来ます。具体的に\(w(t)\)について解けば、

と解けてしまうので、

と解を求める事が出来ます。

直流電圧源の場合


RC回路の時と同様に直流電圧源を考えます。電圧源\(v(t)\)は
\(\displaystyle
v(t)=E_0\theta(t-t_0)
\)
と書くことが出来ます。式(15)に代入して、

と書けます。抵抗\(R\), インダクタンス\(L\)は正ですので、\(t\to\infty\)の漸近でコイルの影響はなくなる、ということを表しています。

RLC直列回路


最後に、RLC直列回路を考えます。

回路方程式は

です。今までと同様に、ラプラス変換を施して

を得ます。\(i(t)\)のラプラス変換である\(F(s)\)について解けば、

を得ますので、\(i(t)\)は

と書けます。\(w(t)\)は

と書くことが出来ます。
これから、\(w(t)\)を具体的に求めるのですが、(式(21b)の括弧内の分母)=0が異なる2つの解を持つか、重解を持つかで場合分けをしなければなりません。
まず、異なる2つの解を持つ場合について計算し、その後、重解を持つ場合について考えましょう。

\(s^2+\frac{R}{L}s+\frac{1}{CL}=(s-\alpha)(s-\beta),~~\alpha\ne \beta\)の場合


(式(21b)の分母)=0 が異なる2つの解\(\alpha, \beta\)を持つ場合を考えます。すると、\(W(s)\)は

と部分分数分解することが出来ます。ここで、式が長くなるので\(x=\frac{\alpha}{\alpha-\beta},~~y=\frac{-\beta}{\alpha-\beta}\)と置きました。

すると、式(21b)は計算できて、

となります。
よって、電流\(i(t)\)は

と求められます。

\(s^2+\frac{R}{L}s+\frac{1}{CL}=(s-\alpha)^2\)の場合


続いて、(式(21b)の分母)=0 が重解\(\alpha\)を持つ場合を考えます。すると、\(W(s)\)は

と部分分数分解することが出来ます。
よって、\(w(t)\)は

と求められるので、電流\(i(t)\)は

と求める事が出来ました。