sikino のすべての投稿

複素関数のメモ

複素関数論のメモです。

  1. 微分
  2. 積分
  3. 参考文献

微分


コーシー・リーマンの関係式(まとめ)


複素関数がコーシー・リーマンの関係式を満たすと、その複素関数は通常の変数通りに微分することが出来ます。複素数\(z=x+iy\)を引数に持つ複素関数\(f(z)=f(x+iy)\)が、実数値関数\(u(x,y), v(x,y)\)を用いて

と書けているとき、コーシー・リーマンの関係式

を満たすならば、\(f(z)\)は正則関数と呼ばれます。ここで、
\(f_x(x,y)=\frac{\partial f}{\partial x}, f_y(x,y)=\frac{\partial f}{\partial y}\)を意味します。式(2)を満たす場合、複素関数\(f(z)\)の複素数\(z\)による微分は

と書くことができ、式(3), (4)は同じ値となります。

コーシー・リーマンの関係式(導出)


複素関数の引数の数は1つであり、それは複素数です。すなわち1変数関数です。故に、複素関数の微分を定義しようとすれば、
\begin{align}
\frac{df(z)}{dz}
\end{align}
という量は1つに定まるはずです。しかし実数の空間で生活している我々には、1つの複素数を実部と虚部に分けて考えてしまう癖のようなものがあります。そのせいで、複素関数は平面に記述される関数であり、実部方向と虚部方向をそれぞれ独立に考えることが出来る2変数関数だと勘違いしてしまいますが、そうではありません。

すなわち、2変数の実数値関数と1変数の複素関数は全く異なる関数であり、同じ様に扱ってはいけません。しかし、我々は1次元上に複素数を表現する方法を持っておらず、本来、1変数の複素関数を2変数の実数値関数のように表現しなければイメージすることが出来ません。この場合、複素関数の表現に何らかの制約がかかることがイメージできるでしょう。その制約が正則である、という条件であり、コーシー・リーマンの関係式なのです。

では実際に本来1変数関数の複素関数を実軸と虚軸という、あたかも2変数関数のように拡張して考えて、微分を考えてみましょう。

前提として、複素関数の微分
\begin{align}
\frac{df(z)}{dz}
\end{align}
は一意に定まるとします。この下で2変数関数として\(f(z)=f(x,y)\)を拡張した場合、微分したい点への近づき方が二通りあることに気が付きます。それは、\(x\)方向と\(y\)方向ということです。2変数関数であれば、偏微分の事です。どんな方向から近づいても微分値は一意に定まると考えているので、両結果は等しくなければなりません。
微分は実数値関数と同様に差分の極限として定義して、
\begin{align}
\frac{df(z)}{dz}\equiv \lim_{\Delta z\to 0}\frac{f(z+\Delta z)-f(z)}{\Delta z}
\end{align}
と書けるとします。ここで、\(\Delta z = \Delta x+i\Delta y\)を意味しており、\(\Delta x, \Delta y\)は実数です。

点\(z=(x, y)\)に実数軸に沿って近づく場合、\(\Delta y=0\)を考えればよいので、微分は差分の極限として
\begin{align}
\frac{df(z)}{dz}=\lim_{\Delta x\to 0}\frac{f(x+\Delta x,y)-f(x,y)}{\Delta x}
\end{align}
と書かれます。ここで\(h\)は実数です。また、虚数軸に沿って近づく場合\(\Delta x=0\)を考えればよいので、
\begin{align}
\frac{df(z)}{dz}=\lim_{\Delta y\to 0}\frac{f(x,y+\Delta y)-f(x,y)}{i\Delta y}
\end{align}
と書かれるはずです。複素関数の微分は一意に決まることから、二つの近づき方を考えたとしても一致していなければなりません。複素関数が\(f(z)=u(x,y)+iv(x,y)\)と実部と虚部があらわに分離されて書かれていたとするならば、
\begin{align}
&\lim_{\Delta x\to 0}\frac{f(x+\Delta x,y)-f(x,y)}{\Delta x}\\
&=\lim_{\Delta x\to 0}\frac{[u(x+\Delta,y)+iv(x+\Delta,y)]-[u(x,y)+iv(x,y)]}{\Delta x} \\
&=\lim_{\Delta x\to 0}\frac{[u(x+\Delta,y)-u(x,y)]+i[v(x+\Delta,y)-v(x,y)]}{\Delta x} \\
&=\frac{\partial u}{\partial x}+i\frac{\partial v}{\partial x} \\
\end{align}
と書けます。また、虚軸方向は
\begin{align}
&\lim_{\Delta y\to 0}\frac{f(x,y+\Delta y)-f(x,y)}{i\Delta y}\\
&=\lim_{\Delta y\to 0}\frac{[u(x,y+\Delta y)+iv(x,y+\Delta y)]-[u(x,y)+iv(x,y)]}{i\Delta y} \\
&=\lim_{\Delta y\to 0}\frac{-i[u(x,y+\Delta y)-u(x,y)]+[v(x,y+\Delta y)-v(x,y)]}{\Delta y} \\
&=-i\frac{\partial u}{\partial y}+\frac{\partial v}{\partial y}
\end{align}
と書くことが出来ます。
両者はもともと同じものであるとしてきましたから、実部と虚部で比較します。
すると、

という関係式が導けるという訳です。これがコーシー・リーマンの関係式と呼ばれるもので、1変数関数の複素関数を実部、虚部という2変数関数に無理矢理拡張した場合に現れる条件です。
複素関数の微分が定義できるとき、複素関数の実部虚部は、この関係式を満たしていなければなりません。ある範囲においてすべての点で複素関数の微分が定義できることを、正則であると言ったり、滑らかである、解析的であると言います。
”複素関数の滑らか”は通常の”滑らか”と意味が異なることに注意しましょう。通常は不連続点が無いことを指しますが、複素関数の場合はそれ以上の制約(コーシー・リーマンの関係式)を含んでいます。

また、ここでは述べませんが、複素関数は1度微分できると何回でも微分することが出来ます。
さらに微分が定義できる正則関数であれば定義域を超えた領域へ関数が拡張できる場合があります。これが解析接続と呼ばれるテクニックです。コーシー・リーマンの関係式を制約として今まで捉えてきたのですが、この制約があるからこそまだ見ぬ定義域外をこの関係式を手掛かりに関数を拡張できるのです。

例として、究極的に滑らかそうに見えるけど微分できない関数を上げましょう。
\begin{align}
f(z)=z=x+iy
\end{align}
はコーシー・リーマンの関係式を満たすため、複素関数として微分が定義でき、正則です。しかし、その複素共役を取った関数
\begin{align}
f(z)=z^*=x-iy
\end{align}
はコーシー・リーマンの関係式を満たさないため、複素関数として微分が定義できず、正則ではありません。近づき方によって微分の値が変わってしまうのです。

通常の微分が可能という意味を明確にするため、例を挙げます。\(f(z)\)が\(z\)の多項式で書かれている場合、微分は

として求める事が出来ます。通常と同じですね。これは、複素数の成分を持たない場合には実軸上で考えたものと同等になるので、実数軸上の微分を自然に拡張したようになっています。

積分


特異点周りの積分


端的にまとめます。
点\(a\)を囲うように複素平面上の閉じた積分経路の方向を\(C\)で指定するとします。\(z\)の多項式で表現される関数\(f(z)=(z-a)^n\)の積分は、

となります。点\(a\)の周りを反時計回りに1周回る場合、\(C\)はパラメータ\(t\)を用いて\(z(t)\)と書くと、
\begin{align}
z(t)=a+re^{it},~(0\le t \le 2\pi)
\end{align}
と書けます。\(r\)は積分経路の半径を意味します。式(6)の計算をすれば、

となります。\(n=-1\)の時、右辺の積分は\(2\pi\)に等しい事は明らかです。それ以外の\(n\ne -1\)ならば、

となります。よって、

が導かれます。

留数定理


一つ前の、ある点周りの積分をもう少し一般的に考えましょう。\(z=a\)を中心とした閉曲線に沿った複素関数\(f(z)\)の複素積分\(\oint_C f(z)dz\)は

と書けます。ここで\(\text{Res}[f,a]\)は留数と呼ばれる量です。これから留数について述べていきます。

\(\text{Res}[f,a]\)は、\(f(z)\)が点\(a\)周りで

とローラン展開出来るとします。この時、\(a_{-1}\)を留数と呼びます。すなわち、

を意味します。

では、留数を実際に求めていきましょう。\(f(z)\)が点\(a\)周りでローラン展開可能な場合、\(n\)の下限が重要になります。下限を極の位数\(k\)と表現します。\(k\)は、

を満たすときの値として定義できます。

もし、関数\(f(z)\)が\(z=a\)で\(k\)位の極であれば、\(z=a\)周りで

と展開できます。
留数は\(a_{-1}\)なので、\(a_{-1}\)について式変形を行っていけば、

と計算できます。

参考文献


原惟行、松永秀章著、『複素解析入門』(第3刷2010, 初版2007)

導出した答えが妥当なものか確かめる方法

自分で導出した式があったとして、それは正しいのでしょうか?
どうやって合っているのだ、と自信を持てばよいのでしょう。
物理の理論研究を行う上で重要な能力となります。

教育の範疇では、教科書なり参考書なりを見て答えを見つければ良いのですが、研究では、貴方が導出した式は世界中探しても何処にもない式なのです。答えの無い解が合っていることを確認する。これほど研究で必要になる能力は無いと思います。この問いに対する答えは、

その式は当たり前を記述しているのか

を確かめることです。
 物理学を行う上では、実験や数値実験などから、ある程度どういう振る舞いをする答えなのか分かっていることがほとんどです。理論を組み立てて、導いた結果から初めて分かることはあまりないのです。
 しかし論文や発表、教科書等においては発表する際には、分かっていない現象を定式化して実験orシミュレーションと一致しました!凄いでしょ!という順序で結果を見せるので誤解するかもしれません。この構成の方が見栄えがしますし、発表が物語みたいに立てられるので人を引き込めるのです。
 実際の理論研究では結果が分かっていて、どうやって基本方程式からその結果を導出するか?その間を埋める作業が研究となっていることが多いのです。それっぽいのが導出できたら、式が妥当なものなのかを精査していくわけです。この過程で、本来なるべき式の振る舞いからずれていて、計算過程に間違いが無いのだとしたら、それが新しい発見に繋がったりするわけです。また、答えがどういう状況で現れるのか分かっているので、適切な近似が想いつくのです(本当の天才は違うでしょうが…)。

なので、物理の理論研究では、分からない答えをどうやって評価するのかが重要になるわけです。この指標を本稿で取り上げたいと思います。
数学の厳密性はさほど考えないので、数学者からは怒られてしまうかもしれませんけれど…

  1. 極端な場合を考える
  2. 特別な点を考える
  3. 特異点を評価する
  4. 次元を確認する
  5. 漸近形を考える
  6. 当たり前の振る舞いをする関数系を持っているか確かめる
  7. 保存則を満たしているか
  8. 複数の関数の和で書けている場合、それぞれの項に意味を見出せるか
  9. 数値計算を行う
  10. 近くの人に、自分の導出に沿って解いてもらう

極端な場合を考える


もしも答えAを導いた場合、極端な場合を考えましょう。
極端な場合、物理的な描像が思い浮かびやすく答えが想像できることが多いです。
そしてその振る舞いが、Aを導いた際に使用した仮定に反しない範囲で説明できるか、ということを確認しましょう。
極端な場合とは、あるパラメータを0または無限大にした時などです。

古典的な単振り子の減衰振動なんかよい例かもしれません。
減衰振動を頑張って解いた結果、振動の抑制を表す正のパラメータ\(\gamma\)(例えば粘性率のようなパラメータ)を使って
\begin{equation}
f(t)=Ce^{-\gamma t}\cos(\omega t),~~\omega =(\text{const})
\end{equation}
と完全に書ける、と導いたとしましょう(上式は間違っていますが、減衰振動の解析解を知らないとして間違っていることを確認する作業を行います)。

まずは初期条件です。時刻\(t=0\)である位置\(f(t=0)\)から手を放す時を考えると、それが表現できている出来るでしょうか?
これは、解に対して\(f(t=0)\)が自由度を含んで\((-\infty, \infty)\)の範囲にいることができるのかという問いです(問題によっては\(\infty\)の位置にいることなど出来ないかもしれません。その場合は適切な範囲で考えなければなりません)。
\(t=0\)を代入すれば、
\begin{equation}
f(t=0)=C
\end{equation}
となって表現できるわけです。まずは問題ないですね。初速度は
\begin{equation}
f'(t=0)=0
\end{equation}
となります。ということは、この解は初速度ゼロを仮定した場合の解ということになります。初速度がある場合、この解は正しくないということが分かります。

反対に、解が
\begin{equation}
f(t)=Ce^{-\gamma t}\sin(\omega t)
\end{equation}
と書かれていた場合は\(t=0\)で\(f(t=0)=0\)になるわけですから、原点にいなければなりません。その代わり、速度に任意のパラメータを含むようになるわけです。

ここで、解を初速度ゼロで放すという過程を入れて導いたのであればあっている解ですが、そうでなければこの時点で間違いであることが分かります。

続いて\(\gamma\)を考えます。\(\gamma = 0\)であれば減衰などない訳ですから、単なる単振動の振る舞いをすることが容易に判明します。
これを導いた式に代入すると、
\begin{equation}
\lim_{\gamma\to 0}f(t)= C \cos(\omega t)
\end{equation}
となるので、想定に見合った振る舞いをしていることが分かります。

反対に\(\gamma \to \infty\)の時、すなわち非常に粘性の高い中にいる場合を考えますとこの場合も容易に想像が出来て、振動はしないで無限の時間を掛けてゆっくりと平衡位置に向かうだろう、ということが想像できるでしょう。
この振る舞いを式が記述しているかを確かめます。
\(\gamma \to \infty\)を代入しますと、
\begin{eqnarray}
\lim_{\gamma \to \infty} f(t)&=& Ce^{-\gamma t}\cos(\omega t)
&\to & 0
\end{eqnarray}
となることが分かります。粘性の高い物質の中にいるのですから振動はせず、\(f(t)=0\)に行く、という事を言っていますが、そこに至るまでに位置\(y=C\)から無限小の時間で\(y=0\)に向かうと言っているのです。速度がめちゃくちゃ早くなることはおかしい振る舞いです。これは想定には合わず、どこかで間違いがあることに気が付けます。

別の例として光がスリットを通過する問題を挙げましょう。
\((x,y)\)の2次元平面を考えて、\(x\to -\infty\)から平面波がやってくる場合を考えます。ここで、\(x=0\)に幅\(a\)のスリットを平面波が通る問題があるとします。この時、\(a\)の関数として\(x\gt 0\)の領域が書きあらわされて導出できたとしましょう。
すると、もし\(a\to \infty\)であれば、スリットは無いと同じなので、関数の格好は平面波として書けていなければなりません。また、トンネル効果やエバネッセント波による染み出しを定式化していない基本方程式から出発したのであれば、\(a\to 0\)の振る舞いは波が全く存在しない、すなわち0になっていなければなりません。
トンネルを記述していたならば、全領域を積分した結果が確率密度に等しいのか?を調べます。エバネッセント波であるならば\(x\to \infty\)の無限遠方でゼロになる振る舞いをするか?を調べるべきです。

更に別の例であれば、気体や物性物理学で用いられる考えでしょう。
例えば有限温度の問題を考えているのであれば、絶対零度の場合には良く知られた解に落ち着くはずです。もし明らかにゼロにできないような近似をしているのならば、\(T\to 0\)もしくは\(k_B T\to 0\)で導出した答えを級数展開しましょう。\(exp(aT)\to 1+aT+O(T^2)\)のように、です。もしも出来ないのであればなにかしら間違えた式変形をしている可能性があります。もちろん、特別な有限温度を中心に展開するような近似をしているのであれば、絶対零度を考えるのは意味がありません。別のパラメータがゼロにすることが出来るのか考えましょう。

特別な点を考える


特別なパラメータが何かしら定義されている時、例えばそのパラメータがエネルギー保存などから導かれているのであれば、それを代入することで妥当性を評価できます。また、初期条件を与えた時、無理のない値となるかで判断することが出来ます。
例えば光学の世界では焦点が分かっている時、その近傍では何か特徴的な振る舞いが起こるはずです。光の強度が発散するとか、光の強度が極大値を持つのでgradがゼロになるとか、そんな地点です。そういった特別な点で確かめるのです。

また、電磁気学や量子力学では点電荷や、デルタ関数型のポテンシャル等の超関数を考えたりします。
これは、もし積分の格好で答えの関数が表現されていたとき、デルタ関数であれば積分が(ほとんどの場合簡単に)計算できてしまうので便利なのです。現実にはあり得なくてもデルタ関数はダメ!という仮定を入れずに定式化してきたのであれば、代入してみるべきです。デルタ関数が扱えない場合は、極微小な幅を持った離散的な点電荷、デルタ関数でもいいでしょう。積分はその(中心の値)×(幅)として大雑把に計算を進めて行けばいいのです。この時は積分の値が適当な値に規格化しておいて考える必要があります。

振り子が良い例かもしれませんね。角度が大きい振り子を真面目に解いて、結果を得たとします。
最下点から初速度を持たせて計算するとき、ポテンシャルエネルギーと同じ運動エネルギーで射出するならば、頂上で止まるはずです。
この時の初速度はエネルギー保存則から
(初期運動エネルギー)=(頂上にある時のポテンシャルエネルギー)
で導出できるはずで、そうなっているのかを確かめるのです。
結果は頂上で永久に止まる解になるはずで、それが定式化されているかを確かめるのです。こういった特別な積分は完全楕円積分など、名前がついていて有名な積分として解ける場合が非常に多いです。

特異点を評価する


例えば、運動方程式を解いた結果、答えが
\begin{eqnarray}
f(t)&=& \frac{g(t)}{t-a}
\end{eqnarray}
みたいな形だったとしましょう。
この場合、\(g(t)\)が分母とキャンセルしない限り、時刻\(t=a\)で\(f(t)\to \infty\)になります。考えている近似の範囲内で起こる現象で\(f(t)=\infty\), もしくは限りなく大きくなる振る舞いがあり得るのか?を調べます。そんなことが起こらないのであれば導出が間違えている、もしくは\(t=a\)近傍の時刻では用いた近似の適応範囲を超えており、使えないことを示しています。

また、複数の関数の積や和で書けている時、とこかの関数がたまたま1やゼロになるような点を見つけてその値を代入してみる、ということも有用です。そういった点は系を特徴づける値になっていることが多く、特徴的な時刻や位置、状況になっていることが多いです。

次元を確認する


これもかなり有効な手段であり、計算ミスが有ったかどうかを判別する手助けになります。
例えば指数関数が出てきた場合、指数の肩は必ず無次元でなければなりません。そのような形式になっているかを確認します。

無次元化して解いている場合は少し判別が難しくなりますが、同じ量の比なはずなのに、分子はxの2次式、分母は1次式みたいな状況になっていると何故2次係数が無くなっている/現れているのかを確認すると良いでしょう。

漸近形を考える


ある微分方程式を解いて求めた答えがあったとします。
その時、微分方程式の時点で漸近形を考えるという近似を入れると簡単に解けてしまう場合があります。
例えば原点近傍で展開してあるパラメータのオーダーで評価するとか、漸近ではある項が落ちてしまうので簡単に解ける、とかです。

当たり前の振る舞いをする関数系を持っているか


例えば振動する現象を表したいのに、導出した結果が\(\exp(i\omega t), \sin(\omega t), \cos(\omega t)\)等を持っていないのであれば、間違えている可能性が高いです。
しかし、指数関数や三角関数は\(\omega\)が虚数をとり得るのであれば、振る舞いが変わるため、十分な精査が必要です。

保存則を満たしているか


確認するために、系の保存則が満たされているのか調べることも良い方法です。
中心力しか存在しない系を計算しているならば角運動量を、外界とのエネルギーのやり取りが無いのであればエネルギー保存則を満たしているのかを調べる等です。もし、2つの粒子を考えているならば、2つの系が持っているエネルギーの和が変わらないや、与えたエネルギー分上昇しているのかなどです。

数値的に確かめることも可能なので、ある程度の指標として用いることが出来ます。
この手助けとしては系の対称性と保存量の関係であるネーターの定理が大きな手助けとなります。

複数の関数の和で書けている時、それぞれの項の持つ意味に妥当性を見出せるか


とくに物理学の場合に役立ちます。もっとも幅広く使われる例は線形結合でしょう。
解の線形結合として、導出した解が書けているのであれば、それぞれの項だけが存在する場合ももちろん解です。
量子力学で変分原理を用いているのであれば、各々の項は固有状態となっているはずです。もしも固有状態を利用して展開しているにもかかわらず、固有状態ではない項が出てくるのであればそれは間違いです。

また、電気回路の過渡現象でも現れるでしょう。パルスの立ち上がりと共にゆっくり消えてしまう項と振動する項が現れるなど、それぞれの項の持つ意味があっても不思議ではないか考えることが出来ます。

数値計算を行う


式とにらめっこしてもどうしても近似が見付かりそうもなくて、数値計算が得意ならば、ぜひ数値計算を行うと良いでしょう。
典型的なパラメータを代入して計算するのです。少なくとも何かの傾向は見えるはずです。
それをみて妥当なのか判断するのも良いでしょう。

近くの人に、自分の導出に沿って解いてもらう


近くにある程度同じことをやっている人がいるならば非常に幸運です。もしそのが信頼できる人で、関係が悪くなければ、焼肉や寿司を奢るから導出仮定を辿ってみて!と言ってみましょう。
アイデアをあげたくなければ、ある程度式をぼかすか不安な一部分だけでも良いでしょう。
貴方の労力がゼロで、ある程度指摘がもらえることが多いです。確認する方法がもらえるかもしれません。
勿論、その人の言うことを全て信用してはいけません。

「導出した」と「式を追う」では、天と地ほどの理解度の差があります。数か月かけて導いた結果を当たり前と言われるかもしれませんが、こらえましょう。
心の中で、(でもあなたはそんな当たり前なことに今まで気が付けなかったんですね)とでもほくそ笑むだけにしましょう。逆もしかりですから。

熱力学

  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]。この結果、
 ・ホップを掛けるほど初速が上がる
 ・破裂音が鳴る
という効果が得られるらしいです。これらの効果を物理的な言葉に置き換えます。

「ホップを掛ける」ということは、BB弾がホップを掛ける為のゴムの所で止まりそこを抜けるために力が必要、ということです。
すなわち、内圧がある程度高まらないと動き出さないということを意味します。
「破裂音が鳴る」ということは、射出時に外気圧と内圧の差が大きいことを意味します。
また、短いバレルにするのは初速が上がりすぎるためだそうです。

以上の振る舞いがシミュレーションで再現できるでしょうか。
グラフの横軸を射出に必要な内圧にとり、バレル長を10cmに設定した時の射出速度と射出時の内圧をプロットしたのが以下の図です。

確かによく言われている通り、BB弾が動き出すまでに力が必要なほど(ホップを掛けるほど)、初速が上がる振る舞いが再現できています。
ある程度のホップの強さからは頭打ちになるようです。このある程度の強さが実際のエアガンで掛けられる強さなのかはわかりません。
射出時の内圧も3気圧と非常に高くなります。外気圧は1気圧なので、これだけ差があれば破裂音が鳴るのも納得でしょう。

流速チューン?のパラメータ設定は以下の通りです。通常の設定と比べ、ピストンヘッド+ばねの重さを22→30g, ばね定数を431→800kg s-2に変更しました。
————————-
BB弾直径 5.95mm
BB弾重さ 0.25g
バレルの直径 6.05mm
バレル長 10cm
シリンダー直径 24mm
ピストンヘッド+ばねの重さ 30g
ばね自然長 150mm
ばね定数 800kg s-2
————————-
流速チューンだ!と言っている方々のブログ等を拝見しましたが、なかなかばねの強さやピストンヘッドの重さ等を詳細に詳しく述べているページはありませんでしたので、上記パラメータはなんとなく設定しています。
また、流速チューン?とはてなマークを付けているのは、具体的なパラメータ例が無いことや、定義が個人によって違うからです。

追記)
通常のばねや、ピストンヘッドの重さでもホップを掛けるほど初速が上がるという傾向は変わらないことが分かりました。
すなわち、上記の傾向はバレル長が短いことによる普遍的な傾向だ、ということです。バレルが短ければホップを掛けるほど初速が上がる、という意味です。

上の計算は以下のパラメータで計算を行いました。
————————-
BB弾直径 5.95mm
BB弾重さ 0.25g
バレルの直径 6.05mm
バレル長 10cm
シリンダー直径 24mm
ピストンヘッド+ばねの重さ 22g
ばね自然長 150mm
ばね定数 431kg s-2
————————-

以上の計算は、下記ページで配布しているエクセルシート(ver1.3以降)にて計算することが出来ます。
弾道計算のコード(Excel)

参考資料


[1]流速チューンの原理と製作例、注意点

弾道計算のコード(Excel)

Excelのマクロ(VBA)で弾道計算を行うコードを作成しました。

https://slpr.sakura.ne.jp/qp/supplement_data/BBullet_course_excel/BBbullet_course_ver1.41.xlsm
上記プログラムの利用の際は連絡先と著作権についてをご覧ください。

以下の微分方程式を計算します。

エクセルシート”bullet”
\(
\begin{align}
m\frac{d^2 \vec{r}}{dt^2}&=
-mg\vec{k}-\frac{1}{2}C_d \rho \pi R^2 |\vec{V}|\vec{V}
-C_l \frac{4}{3}\pi R^3 2\rho |\vec{\omega}| |\vec{V}|\frac{\vec{V}\times\vec{\omega}}{|\vec{V}\times\vec{\omega}|} \\
\frac{d\mathbf{\omega}}{dt}&=N/I \\
N&=\frac{\rho C_f}{2}R\frac{4\pi R^2}{2}\left(|v_{up}|v_{up}+|v_{down}|v_{down}\right)
\end{align}
\)
回転の減衰は本来、弾道計算(BB弾)の理論で導いた通り積分を行わなければなりませんが、簡単にするためにBB弾を真四角の箱として近似して導いています。
この近似は本来の積分値と1パーセント程度の違いしか出てこないので、BB弾の計算には非常に良い近似です。
回転の減衰を取り入れたくない場合、入力パラメータの”回転の減衰”項目をNOにすることでそのように設定が出来ます。
ゼロインを固定した時の弾道の探索機能は実装していないので出来ません。
Excelの設定項目では、Fortranで実装したコードから求めたパラメータをあらかじめ用意しておきました。
設定されたものと書いてあっても、高精度に積分計算をして求めた結果と上の近似をもちいた本エクセルの結果は違くなります。数10cm程度のずれは出てきますので、その点は注意してください。

エクセルシート”barrel”
\(
\begin{align}
m\frac{d^2}{dt^2}x(t)&=\left[P(t)-P_0\right]S_{BB}-\frac{1}{2}C_d \rho \pi R^2 |v(t)|v(t) \\
m_s\frac{d^2}{dt^2}x_0(t)&=-k\left[x_0(t)-x_B-l\right]-[P(t)-P_0]S_s -F_f\label{bbin2} \\
\frac{d}{dt}V(t)&=v(t)S_{\text{barrel}}-v_0(t)S_s+c'(S_{barrel}-\pi R^2)\mathrm{sgn}(P(t)-P_0)\sqrt{\frac{2}{\rho}|P(t)-P_0|}
\end{align}
\)

内部は以下の構造をしていると想定しています。詳細はバレル内部でのBB弾の方程式をご覧ください。

計算手法はどちらも陽的4次ルンゲクッタ法を用いています。

使い方


  1. ファイルのダウンロード
  2. 最新
    https://slpr.sakura.ne.jp/qp/supplement_data/BBullet_course_excel/BBbullet_course_ver1.41.xlsm

    をクリックでダウンロード、もしくは
    右クリックして、”名前を付けてリンク先を保存”からダウンロードしてください。

  3. 計算パラメータの入力
  4. 入力する場所は、色付きのセルです。

    • 弾道計算について

    • 出力結果は大きく分けて2つの項目に分かれています。
      1つ目は定数、射出高さに着弾、地面に着弾の項目は重要と思われるパラメータの抜粋、2つ目は弾道の時間ごとのデータです。
      グラフは、2つ目のデータを元に表示しているだけですので、必要に応じて追加、変更を行ってください。

    • バレル内部の計算について

    • 出力結果は大きく分けて2つの項目に分かれています。
      1つ目は定数、内圧が大気圧に等しい時に射出される時と指定したバレル長におけるパラメータの抜粋、2つ目はBB弾、ピストンの位置、速度の時間ごとのデータです。
      グラフは、2つ目のデータを元に表示しているだけですので、必要に応じて追加、変更を行ってください。

    手入力で入力することも可能ですが、参考設定と書かれたボタンを押すことで、予め設定されたパラメータを読み込ませることもできます。

  5. 計算開始
  6.   入力パラメータを入力した後、計算開始ボタンを押すことで計算が開始されます。

注意点

常識的なパラメータの範囲内でしか動作確認をしておりません。
この結果はあくまで一例です。本シミュレーションを信用して設計・組立を行い、本シミュレーションと結果が違う事があります。その場合私は責任を負いませんので、その点をご了承して使用してください。
BB弾で出来るからと言って、銃弾等の計算はそれなりの信用しかできないことを考慮してください。

追記


バレル内部の計算において、抜弾抵抗の考慮はしていません。すなわちホップが全く掛かっていない状況の計算です。
もし、ホップを掛けた状況を再現したいのであれば、以下のような考えで出来るでしょう。

まず、ホップを掛けた状態でBB弾をバレル内部から引き出す際に必要な力をどうにかして実験的に求めます。
このホップを掛けた際に必要な力\(F\)は、内部圧力\(P\)とBB弾の面積\(S_{BB}\)を用いて
\(
F=PS_{BB}
\)
と書けるはずです。すなわち、内圧がこの力に達した時にBB弾が射出されるとして考え、その内部圧力に達した際にBB弾の運動方程式を解き始める、という方法で可能なはずです。

参考文献

[1] 弾道計算の理論 -シキノート https://slpr.sakura.ne.jp/qp/bullet-course/
[2] BB弾の弾道学 -ハイパー道楽 http://www.hyperdouraku.com/colum/ballistics/index.html

極値を求めるニュートン法

ニュートン法に関するお話です。
”ニュートン法”と呼ばれる方法は
2種類(1.初期値近傍の極値を求めるニュートン法、2. ゼロ点を求めるためのニュートン法
があるようですが、ここでは1. 初期値近傍の極値を求める方のニュートン法についてのお話です。

まとめ

ベクトル\(\mathbf{x}=(x_1, x_2, \cdots, x_N)\)で指定される実数値関数\(f(\mathbf{x})\)の極値を求める事を考えます。
ここで、初期値近傍の点\(\mathbf{x}=\mathbf{x}^{(0)}\)が分かっている時、その近傍にある\(f(\mathbf{x})\)の極小値を与える\(\mathbf{x}=\mathbf{a}\)は適当な反復計算を行い、
\begin{equation}
\mathbf{a}=\lim_{k\to\infty}\mathbf{x}^{(k)}
\end{equation}
で与えられます。極値を求めるニュートン法を用いる場合、\(\mathbf{x}^{(k)}\)は漸化式として与えられ、
\begin{align}
\mathbf{x}^{(k+1)}&=\mathbf{x}^{(k)}+\mathbf{d}^{(k)} \\
\mathbf{d}^{(k)}&=-\left[H f\bigl(\mathbf{x}^{(k)}\bigr)\right]^{-1}\cdot \left[\nabla f\bigl(\mathbf{x}^{(k)}\bigr)\right]
\end{align}
として逐次的に与えられます。ここで、\(H f(\mathbf{x}), \nabla f\bigl(\mathbf{x}\bigr)\)は関数\(f(\mathbf{x})\)のヘッセ行列、傾き(Gradient)を表し、
\begin{equation}
\hat{H}f=
\begin{pmatrix}
\frac{\partial^2 f}{\partial x_1^2} & \frac{\partial^2 f}{\partial x_1 \partial x_2} &\cdots &\frac{\partial^2 f}{\partial x_1 \partial x_n} \\
\frac{\partial^2 f}{\partial x_2 \partial x_1} & \frac{\partial^2 f}{\partial x_2^2} &\cdots & \frac{\partial^2 f}{\partial x_2 \partial x_n} \\
\vdots & \vdots &\ddots & \vdots \\
\frac{\partial^2 f}{\partial x_n \partial x_1} & \frac{\partial^2 f}{\partial x_n \partial x_2} &\cdots &\frac{\partial^2 f}{\partial x_n^2} \\
\end{pmatrix}
,~~\nabla f=
\begin{pmatrix}
\frac{\partial f}{\partial x_1} \\
\frac{\partial f}{\partial x_2} \\
\vdots \\
\frac{\partial f}{\partial x_n}
\end{pmatrix}
\end{equation}
で与えられます。

特に\(f\)が\(x, y\)の2変数関数\(f(x,y)\)であるならば、
\begin{align}
\begin{pmatrix}
x^{(k+1)} \\
y^{(k+1)}
\end{pmatrix}
&=
\begin{pmatrix}
x^{(k)} \\
y^{(k)}
\end{pmatrix}
-\frac{1}{f_{xx}f_{yy}-f_{xy}f_{yx}}
\begin{pmatrix}
f_{yy} & -f_{xy} \\
-f_{yx} & f_{xx}
\end{pmatrix}
\begin{pmatrix}
f_x \\
f_y
\end{pmatrix}
\end{align}
ここで、
\begin{equation}
\hat{H}f
=
\begin{pmatrix}
f_{xx} & f_{xy} \\
f_{yx} & f_{yy}
\end{pmatrix}
,~~
\nabla f=
\begin{pmatrix}
f_x \\
f_y
\end{pmatrix}
\end{equation}
です。ここで、\(f_p\equiv\frac{\partial f}{\partial p}, ~f_{pq}\equiv\frac{\partial^2 f}{\partial p\partial q}, (p, q=x,y)\)と書きました。
微分を差分で近似するとします。\(x, y\)方向の刻み幅を\(h=h_x=h_y\)とすればそれぞれの偏微分は点\((x,y)\)近傍の7点を用いて
\begin{align}
f_x &= \frac{f(x+h, y) – f(x-h, y)}{2h} +O(h^2)\\
f_y &= \frac{f(x, y+h) – f(x, y-h)}{2h} +O(h^2)\\
f_{xx} &= \frac{f(x+h, y) – 2f(x, y)+ f(x-h, y)}{h^2} +O(h^2)\\
f_{xy} &=f_{yx}= \frac{1}{2h^2}\bigl[f(x+h, y) + f(x-h, y) + f(x, y+h) \\
&\hspace{1em}+ f(x, y-h)-2f(x,y)-f(x+h,y+h)-f(x-h,y-h)\bigr]+O(h^2) \\
f_{yy} &= \frac{f(x, y+h) – 2f(x, y)+ f(x, y-h)}{h^2} +O(h^2)
\end{align}
と書けます[2]。

導出


詳しくは述べませんが、1次元の問題における極小値を求めるニュートン法くらいは導いておきましょう。
多次元の場合は拡張すれば良く、考え方は変わりません。多次元では他の文献、例えば[1]を参考にして下さい。

関数\(f(x)\)を\(x=x^{(k)}\)周りでテイラー展開を行えば、
\begin{equation}
f(x)=f\bigl(x^{(k)}\bigr)+f’\bigl(x^{(k)}\bigr)\bigl(x-x^{(k)}\bigr)+\frac{1}{2}f^{\prime\prime}\bigl(x^{(k)}\bigr)\bigl(x-x^{(k)}\bigr)^2+O(\Delta^3)
\end{equation}
となります。ここで、\(\Delta=|x-x^{(k)}|\)とおきました。
今、欲しい解\(a\)は関数\(f(x)\)の極値なので、
\(\frac{df}{dx}=0\)を満たす\(x\)を探せばよいのです。よって、両辺を\(x\)で微分して、
\begin{align}
\frac{df(x)}{dx}=f’\bigl(x^{(k)}\bigr)+f^{\prime\prime}\bigl(x^{(k)}\bigr)\bigl(x-x^{(k)}\bigr)+O(\Delta^2)
\end{align}
を得ます。これが\(0\)になる\(x\)が解\(a\)なので、
\begin{equation}
\left. \frac{df(x)}{dx}\right|_{x=a}=f’\bigl(x^{(k)}\bigr)+f^{\prime\prime}\bigl(x^{(k)}\bigr)\bigl(a-x^{(k)}\bigr)+O(\Delta^2) = 0
\end{equation}
となります。今、\(O(\Delta^2)\)を無視すれば得られる解\(a\)は近似解\(a\approx x^{(k+1)}\)となるので、
\begin{equation}
f’\bigl(x^{(k)}\bigr)+f^{\prime\prime}\bigl(x^{(k)}\bigr)\bigl(x^{(k+1)}-x^{(k)}\bigr) = 0 \\
\end{equation}
を解いて
\begin{equation}
x^{(k+1)} = x^{(k)}-\frac{f’\bigl(x^{(k)}\bigr)}{f^{\prime\prime}\bigl(x^{(k)}\bigr)}
\end{equation}
となります。…ここまで来てしまうとニュートン・ラフソン法と同じですね。関数値がゼロになる点を探すのか、導関数値がゼロになる点を探すのかの違いなだけです。
多次元だと行列になってしまいシンプルではない形になり難しく見えますが、1次元では非常に単純です。

問題は1階(傾き)、2階導関数(ヘッセ行列)をどのように求めるかです。厳密な値が無い場合は近似的にしか得るしかありません。基本的にこれらの計算量は膨大になります。
最も簡単な方法は微分を差分に置き換えてしまうことですが、差分は余り良い方法とは言えません。

厳密なヘッセ行列は、例えば自動微分の考えで得ることが出来ます。例えば、Hyper-dual numbersによる二階偏微分の計算を参考にして得ることが出来ます。

また、引数が複素数で計算結果が実数で返ってくる関数\(f(z)\)の最小値を求めたいのであれば、\(z=x+iy\)として、\(f(x,y)\)の二変数関数としてその極致を探せば良いです。

2次元のプログラム


ここでは2次元の問題に限り、傾き、ヘッセ行列を差分で近似した場合のプログラムを記述します。

下のプログラムは関数
\begin{equation}
f(x,y)=x e^{-(x^2+y^2)/2}
\end{equation}
の、点\((x_0, y_0)=(-1.2, -0.3)\)近傍にある極値を求めるプログラムです。
グラフで示せば、以下の通りになります。
この関数の極値は黒点、初期値は赤点で示しています。

以下に実行例を示します。

$ gfortran main.f90
$ ./a.out
  -1.2000000000000000      -0.29999999999999999      -0.55840071716917605    
  -1.0000000016663813        3.4150669573473641E-013 -0.60653065971263342      4
$

1行目1,2,3列目はそれぞれ\(x,y\)の初期値, 関数\(f(x,y)\)の値を示しています。
また、2行目1,2,3,4列目は求まった\(x,y\)の値, 関数\(f(x,y)\)の値、そして収束に至るまでに行った繰り返しの回数を示しています。

参考文献


[1]塩浦昭義, 数理計画法 第13回
[2]Abramowitz and Stegun, “Handbook of Mathematical Functions”, http://people.math.sfu.ca/~cbm/aands/page_883.htm, http://people.math.sfu.ca/~cbm/aands/page_884.htm

反復計算における収束判定について

収束判定


漸化式
\begin{equation}
x^{(k+1)}=x^{(k)}+\Delta^{(k)}
\end{equation}
に従って、\(\lim_{k\to\infty}x^{(k)}=a\)が適当な定数に収束する問題を考えます。
絶対誤差\(\varepsilon_A\)、相対誤差\(\varepsilon_R\)とすると、数値計算を終了する時は、以下の不等式が満たされる時にすると良いです。
\begin{equation}
|x^{(k+1)}-x^{(k)}|\lt \varepsilon_A + \varepsilon_R(|x^{(k)}|+|x^{(k+1)}|)
\end{equation}

判定式の意味


上記判定式の意味を考えましょう。上の判定式は絶対誤差による評価と相対誤差による評価をまとめて評価する式になっています。

相対誤差が重要な場合

本当の解が\(a=100\)であるとします。その時、\(x^{(k)}=100.3, x^{(k+1)}=100.2\)であったと仮定すると、大雑把に計算して
\begin{align}
|x^{(k+1)}-x^{(k)}| &\lt \varepsilon_A + \varepsilon_R(|x^{(k)}|+|x^{(k+1)}|) \\
0.1 &\lt \varepsilon_A + \varepsilon_R\times 100
\end{align}
となります。右辺は本来2倍ですが、2倍はそんなに重要ではないので無視しています。もし仮に、\(\varepsilon\equiv\varepsilon_A=\varepsilon_R\)としていれば、
\begin{align}
0.1 &\lt 100\varepsilon_R
\end{align}
と近似できます。以上から、本当の解が1より十分大きい場合、相対誤差で収束を評価する式になっていると言えます。
つまり上記の問題設定の場合、\(\varepsilon_R\gt 10^{-3}\)にしてあれば”収束した”として判定します。

絶対誤差が重要な場合

本当の解が\(a=0.01000\)であるとします。その時、\(x^{(k)}=0.01002, x^{(k+1)}=0.01001\)であったと仮定すると、大雑把に計算して
\begin{align}
|x^{(k+1)}-x^{(k)}| &\lt \varepsilon_A + \varepsilon_R(|x^{(k)}|+|x^{(k+1)}|) \\
0.00001 &\lt \varepsilon_A + \varepsilon_R\times 0.01
\end{align}
となります。右辺は本来2倍ですが、2倍はそんなに重要ではないので無視しています。もし仮に、\(\varepsilon\equiv\varepsilon_A=\varepsilon_R\)としていれば、
\begin{align}
0.00001 &\lt 0.01\varepsilon_A
\end{align}
と近似できます。以上から、本当の解が1より十分小さい場合、絶対誤差で収束を評価する式になっていると言えます。
つまり上記の問題設定の場合、\(\varepsilon_A\gt 10^{-3}\)にしてあれば”収束した”として判定します。

参考


[1]杉山耕一朗, OBOROの収束判定条件の設定

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

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

束縛されている時と反発するときは、式(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)\)は

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