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

熱力学

  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)\)は

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

質点と壁との反発を表す運動方程式

質点が壁に衝突し、反発することを数式で表現します。

  1. 壁の定義
  2. 運動方程式の導出
  3. プログラム
  4. 実行結果
  5. 参考文献

壁の定義


壁とは、壁の法線方向に平行な質点の速度成分を反転させるデバイスである。
と定義します。

運動方程式の導出


まず議論を簡単にする為に一次元の運動について考え、その後多次元の運動について定式化をしていきます。

1次元の壁との衝突を表す運動方程式


ポテンシャル\(V(x)\)保存力下の質点の運動を考えます。
壁が
\(f(x,t)=0\)
で表現されているとします。
壁との衝突では位置は連続、速度は不連続な振る舞いをすると考えると、壁に衝突した時、力はデルタ関数の振る舞いをしていることが予想できます。
よって、未知の定数\(c\)を用いて、

と書くことが出来ます。
ここで、時刻\(t=t_i\)は、
\(
f(x(t),t)=0~~\to~~t=t_i
\)

を満たす時刻(壁との衝突時刻)です。\(f(x,t)\)の時間微分は、

と書くことが出来ます。
\(c\)を定めるために式(2)の両辺を、時刻\(t=t_j\)周りを微小時間\(\Delta\)で積分します。すると

を得ます。ここで、保存力の時間変化は連続であることを仮定します。すなわち、

がいかなる時刻\(t\)で成立するとします。計算を進めると、

という結果が導かれます。ここで、\(\delta_{i,j}\)はクロネッカーのデルタを表します。式(6)を式(2)の右辺に代入し、\(c\)を消去すると、運動方程式

を得ます。式(7)はこのままでは解くことが出来ません。なぜなら、未来の時刻の速度\(v(t_i+0)\)が含まれているからです。これをどうにかするには、新たな条件式、すなわち反発前後の条件式が必要になります。

壁と衝突する場合、質点の衝突前後の速度\(v(t_i-0), v(t_i+0)\)と、壁の衝突前後の速度\(v_\text{w}(t_i-0), v_\text{w}(t_i+0)\)の間には、反発係数\(e\)を用いて

の関係があると予想します。今、壁の質量が無限大であり、衝突前後で速度変化がない場合を考えると、壁の速度は、質点との衝突によって変化しないと考えます。すなわち

が成立していると考えます。すると、衝突後の質点の速度は、\(e\)を用いて

と書けます。よって、式(7)に代入して、

を得ます。

最後に、位置\(x(t_i)\)における壁の速度\(v_\text{w}(t_i)\)を考えましょう。衝突の前後のごく短い時間では、質点の動きは壁に追従すると考えます。すなわち、

が成り立っているとします。衝突時には質点の位置\(x\)と衝突位置は同じであることを注記します。式(12)を書き換えると、衝突が起こる時刻\(t=t_i\)の周りで

が成立しています。このことから、衝突前後のごく短い時間では

が成立します。よって、

から、壁の速度

を得ます。よって、壁との衝突を記述する運動方程式は

となります。

多次元の壁との衝突を表す運動方程式




と質点が衝突することを考えます。運動方程式は

であり、1次元の場合と同様に時刻\(t=t_j\)周りの微小時間で積分して

を得ます。式(20)に代入すれば運動方程式

を得ます。次の節で述べる結果(多次元の反発)を先取りすると、式(22)は

と変形することが出来ます。ここで、\(e\)は反発係数、\(\mathbf{v}_\text{w}(t)\)は壁の速度ベクトル、\(\mathbf{n}\)は壁の法線ベクトルであり、\(\mathbf{n}\)は

と表されます。
実際に数値計算を行う際には、衝突時刻\(t_i\)と位置\(\mathbf{r}=\mathbf{r}_i\)を求めた後、

に従って速度ベクトルを変更すると良いと思います(根拠は特にありません)。

多次元の反発



衝突時刻\(t=t_i\), 点\(\mathbf{r}(t_i)\)で質点が\(f(\mathbf{r},t)=0\)で表される壁に衝突することを考えます。
衝突前、後の質点の速度ベクトルをそれぞれ\(\mathbf{v}(t_i-0), \mathbf{v}(t_i+0),\)と置きます。
位置\(\mathbf{r}_i\)における壁の単位法線ベクトルを\(\mathbf{n}, (\mathbf{n}^2=1)\)と書くと、

と書くことが出来ます。ここで、\(c\)は未知の定数でこれから定めていきます。
衝突では壁の法線方向の速度成分のみが変化すると仮定しているので、

が成立します。ここで、\(e\)は反発係数で\(0\leq e\leq 1\)である(\(e=0\):完全非弾性衝突、\(e=1\):完全弾性衝突。
また、\(v_\perp, v_{\text{w}\perp}\)は\(\mathbf{v}, \mathbf{v}_\text{w}\)の、壁の法線方向の成分であり、

と書きあらわすことが出来ます。
式(27)の\(c\)を求めましょう。式(27)の両辺に\(\mathbf{n}^\mathsf{T}\)を掛けて

のように得ます。よって衝突後の質点の速度は

と表すことが出来ます。壁が\(f(\mathbf{r},t)=0\)を満たす線と表現されていれば、時刻\(t=t_i\)で\(\mathbf{n}\)は

と求める事が出来ます。
続いて壁面の速度\(\mathbf{v}_\text{w}(t)\)を考えます。1次元の場合と同様に衝突前後のごく短い時間では質点の動きは壁に追従すると予想します。すなわち、

が成立すると考えます。ここで、式(35)の\(\mathbf{r}(t)\)は質点の位置ベクトルです。すなわち、衝突が起こる時刻\(t=t_i\)の周りで

が成立しており、衝突のごく短い時間では

が得られます。よって

を得ますので、壁の速度

を得ることが出来ます。

プログラム


2次元平面の運動に対するFortran90のプログラムはこちらです。
時間発展は、刻み幅制御陽的ルンゲクッタ法
衝突時刻を求める際の根の探索には、Anderson-Björk’s法を用いています。

壁の関数(サブルーチン fw)とその偏微分(サブルーチン pwf)は手で入力しています。

壁との衝突は、関数の”符号が変わった時”で判定しているので、法線の方向はどちらでも構いません。
すなわち、
\(
f(x,y,t)=\pm g(x,y,t)
\)

は同じです。

プログラム

gnuplot用のスクリプト

実行結果


プログラムを実行した結果です。
適当な壁を定義して実行しています。

重力なしの運動、動かない壁

楕円の焦点から放たれた質点の軌跡

\(
\begin{align}
x(0)&=4,~~x'(0)=(\text{const}) \\
y(0)&=0,~~y'(0)=(\text{const}’)
\end{align}
\)


\(
\begin{gather}
f(x,y,t)=\left(\frac{x}{l_x}\right)^2+\left(\frac{y}{l_y}\right)^2-1=0,\\
l_x=5,~l_y=3
\end{gather}
\)

重力下の運動、動かない壁


\(
\begin{align}
x(0)&=0,~~x'(0)=2 \\
y(0)&=0,~~y'(0)=7
\end{align}
\)


\(
\begin{align}
f(x,y,t)=y – \sin(x)=0
\end{align}
\)

重力下の運動、動く壁


\(
\begin{align}
x(0)&=4,~~x'(0)=0 \\
y(0)&=4,~~y'(0)=5
\end{align}
\)


\(
\begin{gather}
f(x,y,t)=y – \sin(kx+\phi_x)\sin(\omega t+\phi_t)=0,\\
k=1,~w=1,~\phi_x=-\pi/2,~\phi_t=0
\end{gather}
\)

参考文献


3. 壁との衝突 -物理学の見つけ方

4. 動く壁との衝突 -物理学の見つけ方

量子テレポーテーションのざっくりとした説明

量子テレポーテーションとは、量子状態を異なる場所に移動させる方法です。

良くある間違いとして、
・光速を超えて情報通信が出来る(EPRパラドックス)
・物体そのものが移動する(SFのテレポーテーション、ど〇でもドア等)
・物体のコピーが生成される
が挙げられます。上の3点は正しくありません。

本稿の内容が間違えていたら申し訳ありません。

量子テレポーテーションとは


量子テレポーテーションとは、量子もつれを利用して「地点Aの量子状態Xを地点Bに移すこと」です。
量子力学の基本として、未知である任意の量子状態に対し、それと全く同じ複製を作る事は不可能です(量子複製不可能定理 -wikipedia)。
量子状態の複製は出来ませんが、量子状態の移動は行うことが出来ます。これが量子テレポーテーションなのです。
重要なのは、移動するのは物質ではなく、量子状態であることが重要です。

量子テレポーテーションのアイデア


・送りたい側 
 地点Cから送られてきたもつれた量子1’と送りたい量子状態Sを相互作用させる→弱く観測
・受け取り側 
 地点Cから送られてきたもつれた量子2’に対し、弱く観測した結果に基づいた演算をする。

  1. 地点Cで量子1, 2をもつれ(相互作用)させる
    → 量子1,2は区別できないので、量子1′,2’と書くことにする。
  2. 量子1’を地点A、量子2’を地点Bに届ける。
  3. 地点Aで送りたい量子状態Sと量子1’を相互作用させて”弱く”観測する。
    →送りたい量子状態Sはこの時点で量子1’と区別が出来なくなる。
  4. 地点Aの観測結果を古典的な光の速さで地点Bに伝える。
  5. 地点Bの量子2’に、地点Aの観測結果に依存した操作をする。
  6. 地点Bの量子2’は量子状態Sになっている。

という手順です。

地点Aの観測結果を地点Bに送る、という操作が入るので、意味のある情報(量子状態S)を光速を超えて通信する、ということは出来ないわけです。

また、”弱く”観測するの部分では、量子状態を決めることはしません。しっかり観測してしまうと確定的な情報を得てしまい、ダメだそうです。

物体の移動をともなわないにもかかわらず、”テレポーテーション”と呼ばれる由来は、量子が区別できないことに由来します。
区別することができないので、地点Aの量子状態Sと、全てが終わった後の地点Bに存在する量子状態Sは同じものとして捉えられ、あたかもテレポートしたように見える訳です。

同種の粒子が2つあり、それらが同じ量子状態であるのならば、どちらがオリジナルなのか区別する方法はありません。

2017年には中国のグループが、1000km程度の地上-衛星間の量子テレポーテーションに成功した、との報告があります。
Ji-Gang Ren et al, “Ground-to-satellite quantum teleportation”, arXiv:1707.00934 [quant-ph], (2017) https://arxiv.org/abs/1707.00934

参考文献


量子テレポーテーション https://www.m-nomura.com/st/qteleport.html

ーーー観測するって何なんでしょうね。

束縛条件下の運動 – 自由度がうまく落とせない運動

このページは
束縛条件下の運動 – ホロノミックな束縛と非保存力
の続きです。

拘束がある場合の2次元の運動


ラグランジュの方程式を解く際に適切な座標変換が見付からず自由度が落とせない場合を考えます。その時の運動方程式は

の形で止まります。式(75)において、未知の関数は\(x(t), y(t), \lambda(t)\)の3つであり、条件式は3つなので解くことは出来るはずです。
方針は、式(75)から\(\lambda(t)\)を消去すれば良いのです。

束縛条件\(f(x,y,t)=0\)が成立する点\(x=x_0, y=y_0\)周りで、微小時間\(\Delta t\)の変化を考えます。
点\((x_0, y_0)\)まわりでテーラー展開すれば、

を得ます。よって

が満たされる変化でなければなりません。すなわち、\(x,y\)の時間変化は式(77)をいつも満たしていなければならないのです。両辺を\(\Delta t\)で割って極限をとれば、

が成立することが分かります。式(75)を書きかえれば、

を得ます。式(79c)もまた時間変化しても成立し続けなければなりません。よって、左辺の時間微分を取れば、

を得ます。整理して

を得ます。ここで関係式

を用いました。式(79a), (79b)から未定乗数\(\lambda(t)\)を消去すれば

を得ます。よって、式(81)と式(84)から連立方程式

を立てることが出来ます。これをあらわに解けば、運動方程式

を得ます。ここで、式(86), (87)の右辺の\(-\frac{\partial V}{\partial x}, -\frac{\partial V}{\partial y}\)は\(x, y\)方向の保存力であり、右辺の残りの項は束縛条件によって決まる束縛力であることが分かります。

拘束がある場合の3次元の運動


スカラー関数\(f\)で表される拘束がある条件のもとで, 運動方程式

またはベクトル表記で

について考えます。ここで表記を略すために

の略記を用いました(\(y, z\)方向についても同様)。
また、\(\nabla\)はナブラ演算子で、

です。2次元の場合と同様に、式(88d)の時間微分から

が導けます。ここで、\(\hat{H}f\)は関数\(f\)のヘッセ行列であり、

を表します。式(92)は、式(81)を3次元に拡張した形になっています。

式(88a),(88b),(88c)から未定乗数を消去すると従属な3つの関係式

を得ることが出来ます。この3式で独立なのは2つです。
独立な式として式(94a)と式(94b)を選ぶことにし、式(92)と共に書けば、連立方程式

を得ます。この式をあらわに解けば、運動方程式

を得ることが出来ます。ここで、\(n_x, n_y, n_z\)は壁の法線ベクトル\(\mathbf{n}=\frac{\nabla f}{|\nabla f|}\)の\(x,y,z\)成分です。偏微分ではないことに注意してください。

ベクトル表記であれば

と書くことが出来ます。

単振り子の例


ここでは単振り子を例に挙げ、式(86)に代入した時本当に振り子を表す方程式になっているのか、また束縛条件の選び方に依らず、同じ運動を記述しているのかを例に挙げます。
保存力は重力で

を考えます。

考える束縛条件は
\(
f(x,y,t)=x^2+y^2-l^2=0
\)


\(
f(x,y,t)=\sqrt{x^2+y^2}-l=0
\)

です。

束縛条件1


束縛条件として

を考えます。この束縛条件に対する偏微分を計算して、

を得ます。これを式(86), (87)に代入すれば、

であり、整理すれば

を得ます。

束縛条件2


束縛条件として

を考えます。この束縛条件に対する偏微分を計算して、

を得ます。これを式(86), (87)に代入すれば、

であり、整理すれば

を得ます。

考察


同じ束縛条件なので、式(105)と(109)の示す運動は同じになるはずです。
それを示すには、2式の差異である

を示せればよい、ということです。そのために拘束力に対して座標変換

を考えます。右辺と左辺をそれぞれ計算すれば、

となり、同じ運動を記述していることが分かります。

Fortran90によるプログラム


Fortran90で2次元平面を動く質点の束縛運動のプログラムを作ります。
そのプログラムは
https://slpr.sakura.ne.jp/qp/supplement_data/constraint.tar.gz
に置いてあります。

一応補足しておきますが、初期状態の位置は\(f(x, y, t)=0\)を満たすx,yでなければなりませんし、初期速度も\(f(x, y, t)=0\)の傾きと一致していなければなりません。
上の条件を満たさずとも何も警告は出さず、計算自体は行われますが、その計算は束縛運動ではなくなります。

上のプログラムの中で、変更すると思われるのは以下の2つです。
1. fp2d
2. fw2d

1. fp2dは、ポテンシャル\(V(x,y)\)を記述しており、デフォルトのプログラムでは、重力による位置エネルギー
\(
V(x,y)=mgy
\)
に設定しています。

2. fw2dは、束縛条件\(f(x,y, t)=0\)を記述しており、デフォルトのプログラムでは、
\(
f(x,y,t)=y – [2(cos(1.5*x)-1) + 0.1x^2]=0
\)
に設定しています。

もしも非保存力を入れたい場合は、grkの最後の方を変更してください。

数値計算では、
時間発展は刻み幅制御の陽的ルンゲクッタ、
gradient, Hesse行列、束縛条件に関する時間の二階微分は、Hyper-dual numbersを利用してほぼ厳密な値を得ています。

program main
  ! Author : sikino
  ! Date   : 2019/10/12 (yyyy/mm/dd)
  ! URL    : http://slpr.sakura.ne.jp/qp/constraint-condition3/
  use Parameters
  use Hyperdualmod
  implicit none
  integer::i,N,info,Nt
  double precision::t,h,tol,ta,tb
  double precision,allocatable::x(:),tgrid(:)
  external::grk
  double precision::E,Ek,Ep

  N=4
  allocate(x(1:N))

  ! Time range [ta, tb]
  ta = 0d0
  tb = 30d0
  ! Divide time range equally among Nt
  Nt = 201
  ! Tolerance of the RK method
  tol=1d-8

  ! Initial condition
  x(1)=0d0 ! x (ta)
  x(2)=4d0 ! x'(ta)
  x(3)=0d0 ! y (ta)
  x(4)=0d0 ! y'(ta)

  ! Time grid
  allocate(tgrid(1:Nt))
  do i=1,Nt
     tgrid(i) = (i-1)*(tb-ta)/dble(Nt-1)+ta
  enddo
 
  do i=2,Nt
     info = 0
     t = tgrid(i-1)
     h = (tgrid(i) - t)*0.1d0
     do while(info.le.0)
        call drkf45(grk,t,h,N,x,tgrid(i),info,tol)
     enddo
     
     ! Kinetic energy
     Ek = 0.5d0*mass*(x(2)**2+x(4)**2)
     ! Potential energy
     Ep = mass*g*x(3)
     ! Total energy
     E  = Ek + Ep
     write(10,'(8e25.10e3)')t,x(1),x(2),x(3),x(4),E,Ek,Ep
  enddo
 
  stop
end program main

subroutine fp2d(N,xH,fH)
  use Parameters, only:mass, g
  use Hyperdualmod
  implicit none
  integer,intent(in)::N
  type(Hyperdual),intent(in)::xH(1:N)
  type(Hyperdual),intent(out)::fH

  ! Potential

  type(Hyperdual)::x,y,t
  x = xH(1)
  y = xH(2)
  t = xH(3)
  ! If conservative force,
  !  fH doesn't depend on the time.
  !-----------------------------
 
  fH = mass * g * y
 
  return
end subroutine fp2d

subroutine fw2d(N,xH,fH)
  use Hyperdualmod
  implicit none
  integer,intent(in)::N
  type(Hyperdual),intent(in)::xH(1:N)
  type(Hyperdual),intent(out)::fH

  ! Wall of the function

  type(Hyperdual)::x,y,t
  x = xH(1)
  y = xH(2)
  t = xH(3)
  !-----------------

  fH = y - (2d0*(cos(1.5d0*x)-1d0) + 0.1d0*x*x)
!  fH = (y-5d0)**2 + x**2 - 25d0
 
  return
end subroutine fw2d

デフォルトのまま計算すると、fort.10に以下の出力が現れます。アニメーションは、同じファイルに入っているanime.pltで実行できます。

エネルギーは計算の範囲ではほとんどありません。全力学的エネルギーは、計算時間内でおおよそ5桁目が変化していました。

参考文献


6. 振り子 -物理学の見つけ方
9. 自由な座標 -物理学の見つけ方

束縛条件下の運動 – ホロノミックな束縛と非保存力

このページは
束縛条件下の運動 – ホロノミックな束縛 1
の続きです。

単振り子



ポテンシャルが存在し、ホロノミックな束縛が課されている代表的な問題である、単振り子を考えましょう。
単振り子は、束縛条件下の運動 – ホロノミックな束縛 1で考えた円に束縛されている問題に対してポテンシャルが加わった問題です。

座標系は以下のような極座標系を選びます。

通常の極座標系とは異なります。\(\theta=0\)で鉛直下方向にとる座標系にしたかったので、通常の極座標系とは変えています。

ホロノミックな束縛条件\(f(x,y,t)=0\)を含めたラグランジアンは、デカルト座標系から極座標系に変えて

と書くことが出来ます。注記しますが、ラグランジュの未定乗数は時間に依存していて構いません。
ラグランジアンの偏微分をそれぞれ計算すれば、

ですので、運動方程式と束縛条件

が得られます。連立微分方程式(31)の3本の方程式に対して未知の関数は\(r(t),\theta(t),\lambda(t)\)の3つなので、方程式を立てて\(r(t),\theta(t)\)を求めることが出来ます。

では、伸び縮みしない振り子を考えましょう。束縛条件は

をして表現することが出来ます。束縛条件の一般化座標による偏微分は

ですので、運動方程式

を得ます。式(34)の第2式を変形すると

ですので、これはまさに振り子を表す運動方程式です。

実際に解いてみると、こんなグラフが得られます。

・パラメータ
 \(m=3, l=1, g=9.8\)
\(\theta(0)=\frac{4\pi}{5}, \theta'(0)=0\)
・計算手法
 陰的ルンゲクッタ法

支点が動く振り子



次に、束縛条件が2つある場合を考えましょう。その例題として、
振り子の支点が摩擦の無いレールの上に載っていることを考えます。
ラグランジアンは

です。2つのホロノミックな束縛条件

を考えます。すると、ホロノミックな束縛条件を含めたラグランジアンは、

と書くことが出来ます。
これから伸び縮みしない振り子を考えることを見越して、座標系

を考えます。この新たな座標系でのラグランジアンは、

と書けます。それぞれの偏微分である



を用いると、運動方程式

を得ます。
連立微分方程式(46)の4本の方程式と2本の束縛条件に対して、未知の関数は\(x(t), y(t), r(t), \theta(t), \lambda_0(t), \lambda_1(t)\)の6つなので、方程式を立ててそれぞれの関数を求めることが出来ます。

では、支点が直線状のレールの上に置かれている状況を考えましょう。
支点に関する束縛条件\(f_0\)と, 振り子の紐の長さが変わらないという束縛条件\(f_1\)を考えると

と書くことが出来ます。なので、

とそれぞれ計算できるので、運動方程式

を得ます。求めたい未知の関数は\(x(t), \theta(t)\)です。
式(51a)と(51b)から、\(x(t), \theta(t)\)に関する運動方程式

を得ます。変形すると

ですので、あらわに求めれば、

を得ます。
特に支点の質量が振り子の質量に比べて非常に重い場合、すなわち\(m_0 \gg m_1\)の場合、

になります。第2式はまさに単振り子の方程式となっています。

実際に解いてみると、こんなグラフが得られます。

・パラメータ
 \(m_0=3, m_1=1, l=1, g=9.8\)
\(x(0)=0, x'(0)=0, \theta(0)=\frac{4\pi}{5}, \theta'(0)=0\)
・計算手法
 陰的ルンゲクッタ法

支点が放物線上に束縛されている振り子



座標系は1つ前の問題と同じく、

にとることにします。すると運動方程式は式(46)で表されるので、運動方程式の導出過程を省略することが出来ます。
問題設定から、束縛条件は

と書くことが出来ます。なので、それぞれの偏微分は

と書けますので、運動方程式は

です。式(60a)と(60b)から\(\lambda_0(t)\)を消去すると、

を得ます。更に、束縛条件から座標\(y(t)\)の時間変化は\(x(t)\)を用いて

と書くことが出来ます。式(61)と式(60d)から、\(x(t), \theta(t)\)の運動方程式

を得ます。あらわに解くと式が長くなってしまうので、ここで止めておきます。

実際に解きますと、

というような振る舞いになります。

・パラメータ
 \(m_0=3, m_1=1, l=1, g=9.8\)
 \(x(0)=0, x'(0)=5, \theta(0)=0, \theta'(0)=0\)
・計算手法
 陰的ルンゲクッタ法

単振り子と空気抵抗


束縛条件と、非保存力である空気抵抗が存在する場合を考えます
座標系は

です。ホロノミックな束縛条件を含めたラグランジアンは、

と書けます。

さて、非保存力である空気抵抗が存在することを考えましょう。
非保存力なのでラグランジアンで書くことは出来ませんので、個別に求める必要があります。
空気抵抗の力\(F_x, F_y\)はデカルト座標系において

で掛ける事を既知とします。今知りたいことは、式(64)で表される一般化座標\(r, \theta\)で式(67)がどのように表されるのか?ということです。
\(r, \theta\)において、一般化座標における空気抵抗\(Q_x, Q_y\)は、

で変換することが出来ます。具体的にそれぞれの力を示せば、

ですので\(Q_x, Q_y\)は、

と表せられます。非保存力がある場合、ラグランジュの方程式の右辺にそれを入れれば良いので、

と表されます。よって運動方程式

を得ます。
長さ\(l\)の伸び縮みしない紐に繋がれた振り子をを考えれば、束縛条件は

ですので、運動方程式

を得ます。

実際に解きますと

・パラメータ
 \(m=3, l=1, g=9.8, c_1=0.3, c_2=0.3\)
 \(\theta(0)=0.8\pi, \theta'(0)=0\)
・計算手法
 陰的ルンゲクッタ法

このページはここで終わりです。
続いて、自由度がうまく落とせない、うまい座標系が見付からない場合を考えます。

次のリンク
束縛条件下の運動 – 自由度がうまく落とせない運動

参考文献


田辺行人、品田正樹著『理・工基礎 解析力学』 裳華房

束縛条件下の運動 – ホロノミックな束縛

ホロノミックな束縛を受けている質点の運動を考えます。

ホロノミックな束縛とは、
”束縛条件が一般化座標と時間のみで決まる方程式によって解析的に表現できるもの”
です。空気抵抗のない単振り子などは、半径(一般化座標)が一定という条件のみで決まるので、ホロノミックな束縛に相当します。

ラグランジュの方程式を基本として考えていきます。

スタートはラグランジアン、ゴールは運動方程式を立てる事です。

2次元平面上に束縛されている質点


1つの質量を持った質点が、ある軌道上しか動けない運動に対する運動方程式を導きましょう。
ある軌道上しか動けないという条件は、時刻\(t\)における質点の位置\((x,y)\)が、関数\(f(x,y,t)=0\)を満たす動きしか許されないということで表現することが出来ます。

束縛条件を含めた実効的なラグランジアンは, 未定乗数\(\lambda(t)\)を用いて

と書けます。ラグランジュの方程式

に代入すれば運動方程式

を得ます。

円に束縛されている運動


運動が半径\(l\)の円に束縛されている場合、束縛条件は

として表現されます。もちろん、\(f(x,y,t)=\sqrt{x^2+y^2}-l=0\)としても構いません。運動は束縛条件の形には依らないからです。

式(3)に代入すると、運動方程式

を得ます。方程式(3),(4),(5)の3本の方程式に対して未知の関数は\(x(t),y(t),\lambda(t)\)の3つなので、方程式を立てて\(x(t),y(t)\)を求めることが出来ます。

式(5)の第1式に\(x\)、第2式に\(y\)を掛けて両式を引くと

です。変形すれば

を得ます。どういう解き方でも構わないのですが、半径\(l\)が変わらないことから極座標

を用いることを考えます。すると、式(7)は

と書けます。よって

ここで、\(\omega,\theta_0\)は初期条件から来まる定数です。よって、解

を得ます。ちなみに、\(\lambda(t)\)は

です。

角周波数\(\omega\)で等速回転する筒に入れられた質点が描く軌道



続いて、ホロノミックな束縛条件が時間\(t\)に依存する問題を考えてみましょう。
そこで、角周波数\(\omega\)で等速回転する筒に入れられた質点の軌道を考えます。

極座標表示の方が扱いやすいので、極座標で考えていきます。

ホロノミックな束縛条件を含めた実効的なラグランジアンは、デカルト座標系のラグランジアンから出発して

を得ます。ここで、\(\lambda(t)\)は未定乗数を表します。
計算していけば、運動方程式

を得ます。方程式(18a),(18b),(18c)の3本の方程式に対して未知の関数は\(r(t),\theta(t),\lambda(t)\)の3つなので、方程式を立てて\(x(t),y(t)\)を求めることが出来ます。

角周波数\(\omega\)で等速回転する筒を表現する関数を考えてみましょう。
筒の上の任意の点\((x,y)\)は、

と書くことで表現することが出来ます。よって、束縛条件から

という関係式を得ます。式(21)を運動方程式(18)に代入すれば、

を得ます。式(22a)を解けば、

であるので、解

を得ます。ここで、\(A,B\)は初期条件から決まる定数です。

実際に解いてみると、こんなグラフが得られます。

・パラメータ
\(A=0.01, B=0, w=2\)

このページはここで終わりです。
続いて、ポテンシャルが加わる問題である単振り子と束縛条件を考えていきます。
次のリンク
束縛条件下の運動 – ホロノミックな束縛と非保存力

参考文献


田辺行人、品田正樹著『理・工基礎 解析力学』 裳華房

グリーン関数による1次元の運動方程式の解法

グリーン関数を用いて1次元のニュートンの運動方程式を解きます。

  1. グリーン関数
  2. 重力下の質点
  3. 重力下のばね
    1. 1つ目の選び方
    2. 2つ目の選び方
  4. 注釈
  5. Thanks(フォントについて)

グリーン関数は非斉次微分方程式の特解を求める為に良く利用されます。
やっていることは、
微分方程式の満たす関数を、あたかも解いた風に書くこと
です。
なので実際、グリーン関数を使う、とは単なる式変形(微分方程式→積分方程式)をしているにすぎません。グリーン関数を本当に求め終わった時、はじめて解けた、と言えるのです。

物理現象では、グリーン関数が表現するのは「無限小時間に力が働く衝撃を受けた結果」を表現します。
どんな力も「複数の衝撃を積算したもの」と捉えることが出来ます。
衝撃に対する応答を考えてから、任意の力を衝撃の足し合わせと考え、問題を解きます。
すると、各々の力に対してその都度微分方程式を解いていくよりも
グリーン関数を求めた後の結果を足し合わせて任意の力に適応する、という考えの方が便利そうだと想像できます。

グリーン関数


微分方程式

をグリーン関数を用いて解くことを考えます。
ここで、左辺の\(\hat{D}\)は演算子を表し、例えば

のような形です。
グリーン関数で解くとは、解を(一般解)+(特解)で書くときの特解を与える手段です。
特解なので、何かしらの境界条件の下でグリーン関数を求める訳です。
物理の問題では、”因果律”と呼ばれる、結果は原因の後の時間に起こる、ということを仮定した境界条件を採用します。

全ての演算子を左辺に移してしまいましたが、右辺に時間に関する微分演算子が入っていても問題ありません。
式(1.1)の形式的な解は、グリーン関数\(G(t,t’)\)を用いて

と、式(1.3a),(1.3b)を満たす解を用いて書くことが出来ます。
ここで、\(\delta(x)\)はディラックのデルタ関数です。
実際に解(1.3)は式(1.1)を満たすことを確かめます。
式(1.3)に左から演算子\(\hat{D}\)を作用させます。
ここで、\(\hat{D}\)は\(t\)に関する演算子なので、積分と交換できます。よって

と変形でき、まさに解(1.3)は式(1.1)を満たす解として表されます。

重力下の質点


さて、質量\(m\)を持つ質点が、重力加速度\(g\)の重力下で運動する問題を考えます。運動方程式は

です。質量\(m\)で割って

を得ます。
もちろん、この問題の解は
\(
\displaystyle x(t)=-\frac{1}{2}gt^2+at+b
\)

です。この問題はグリーン関数を用いて解く必要は全くありませんが、解析解が簡単に求められる練習問題として良いでしょう。
グリーン関数は単なる形式ですが、一般的な解法であり、非常に強力です。
しかし、こんな簡単な問題も解けないのであればお話になりません。

さて、式(1.6)の形式的な解は

と書くことが出来ます。ここで、\(x_0(t), G(t,t’)\)は以下の微分方程式を満たします。

\(x_0(t)\)については簡単に解けてしまい、

です。ここで、\(\alpha, \beta\)は\(x(t)\)の初期条件によって決まる定数です。

続いてグリーン関数の微分方程式を解きましょう。
方程式の右辺を見ますと、\(t= t’\)だけで変なことが起こっています。
ですが、\(t\ne t’\)であれば右辺は完全にゼロです。
すなわち、\(G(t,t’)\)は\(t\ne t’\)で、微分方程式

を満たしています。これは簡単に解けてしまい、

です。ここで、\(A_0, A_1, B_0, B_1\)を適当な定数としておきました。
これらの係数間には、とある関係式が成り立っています。それを調べましょう。

デルタ関数を発散する点を含んで積分すると定数1になることを利用します。
式(1.18)の下の式の両辺を時刻\(t=t”\)周りで積分します。
すなわち、微小な時間\(\varepsilon\)を用いて、

を計算します。
もしも、たまたま\(t”=t’\)であれば、積分を実行して

を得ます。これは、時刻\(t=t’\)周りでは、グリーン関数の導関数は、右から極限と左から極限は一致しないことを言っており、不連続になると言っています。
\(t”\ne t’\)であれば、この極限は一致し、グリーン関数の導関数は連続であるということを言っています。

式(1.11)を代入すれば、係数間の関係式

を得ます。
式(1.7)のグリーン関数を含む項は特解ですので、不明な定数は存在しません。
今、4つの未知の定数\(A_0, A_1, B_0, B_1\)に対して式の数は2つなので、全ての係数を決めるには2つ条件式が足りません。
そこで、境界条件を加えます。この条件は式(1.7)の解の形を見て決められます。
運動方程式を解こうとする場合、通常は物理的な要請から、ある時刻の運動はその時刻より過去の運動のみから決定される、という因果律が課されます(※1)。
さて、因果律を仮定すると、事が起こる時刻\(t=t’\)より前には何も起こらない、すなわち\(G(t,t’)=0\)を意味しており、\(A_0=A_1=0\)となります。よって、

を得ます。よって、グリーン関数

を得ます。ここで、\(\theta(x)\)はヘヴィサイド関数を表します。
グリーン関数(t’=5に固定)を書いてみるとこんな感じ。

黒線がグリーン関数、赤線がその導関数です。導関数の右からと左からの極限の差がちょうどデルタ関数の積分、すなわち1になっていることが分かります。

これでグリーン関数が得られました。式(1.16)を使って式(1.7)の積分を計算すると、

を得ます。ここで、\(c\)は任意の値ではないことに注意してください。定まっている値です。

式(1.9)と式(1.17)の結果を代入すれば、運動方程式の解

を得ます。ここで、\(\beta+c\)を改めて\(\beta\)と書きました。
この結果は本当の解に一致します。

重力化のばね


続いて、鉛直方向に置かれたばねと重力の問題を考えましょう。
ばね定数を\(k\)と置けば、運動方程式は

です。質量\(m\)で割り、\(\omega=\sqrt{k/m}\)と置くと、

を得ます。
解は平衡点を中心とした振動運動です。平衡点\(x_0\)は力のつり合い
\(
-mg=kx_0~~~\to~~~x_0=-\frac{g}{\omega^2}
\)
から導くことが出来るので、
\(
\displaystyle x(t)=a\sin(\omega t +\varphi) -\frac{g}{\omega^2}
\)

です。これをグリーン関数を用いて導くことが出来れば目的達成です。

さて、グリーン関数の考え方として少なくとも2通りの方法が考えられます。
どちらの方法でも同じ解を与えるので、どちらをえらんでも構いません。
数学的に解きやすい方法を選ぶのが賢い方法です。

実際に試してみましょう。

1つ目の選び方

微分方程式を

として考えます。
解を

と考えるのです。\(x_0(t), G(t,t’)\)はそれぞれ

の微分方程式を満たします。式(1.9), 式(1.16)の結果を用いれば、容易に求められて、

です。式(1.22)の積分部分は

なので、解

を得ます。
さて、あらわな形ではないものの、解が得られました。
式(1.27)を満たす\(x(t)\)を探さなければならないのですが、
正直なところ解を見つけられる気がしません。これは自己無撞着の方程式であり、
例えば散乱問題等のように、グリーン関数を得るのが難しい時などでこの形式にとどめておく方法がとられます。
あらわに求める事が出来ずとも、他の式へ代入が出来たり、摂動展開なんかに用いられるので、微分方程式の形式よりその後の計算が扱いやすいのです。
ただ、今の目的からそれてしまうので、ここではこれ以上は話しません。

2つ目の選び方

微分方程式を

のように選びます。すると、解は

であり、\(x_0(t), G(t,t’)\)はそれぞれ

の微分方程式を満たします。
\(x_0(t)\)はただのばねの問題なので、

という解を得ます。一方、グリーン関数も\(t=t’\)を境に場合分けすると

を得ます。時刻\(t=t”\)周りの微小時間で積分すれば、

を得て、たまたま\(t”=t’\)であれば、

の関係式を得ます。ここで、グリーン関数自体はいかなる時刻でも右からと左から極限が一致するという事実をもちいました。
因果律を境界条件として使用すれば、係数は

であるので、グリーン関数は

と求められます。ここで、\(n\)は三角関数の周期性からくる任意の整数です。具体的なグリーン関数(t’=5に固定)の格好はこんな感じ。

黒線がグリーン関数、赤線がその導関数です。導関数の右からと左からの極限の差がちょうどデルタ関数の積分、すなわち1になっていることが分かります。

式(1.29)の積分を計算すると

と求められます。三角関数を時刻無限に飛ばした時の値ですが、拡張リーマン・ルベーグの補題から得られる事実を用います(※2)。これによると、

という事実が得られるので、式(1.37)の最後の項はゼロです。よって、

を得ます。
したがって、解

を得ます。これは求めたかった平衡点を中心に振動する運動です。

さて、1つ目の選び方で得た解(1.27)に解(1.42)を代入してみましょう。解(1.27)が本当にばねの問題なら、式は満たされるはずです。
式(1.27)の積分は

と計算できます。ここで、\(d\)を任意の定数ではない定数と置きました。
式(1.27)に代入すると、(RHS:right-hand side, 右辺)

を得ますが、今、\(\alpha’, \beta’\)は自由に決めることが出来ます。
もしも\(\alpha’=0,~~\beta’+d=0\)となるように定数を決めれば、解(1.27)は式(1.42)で表される解と同じで、微分方程式(1.20)の解となっていることが分かります。

注釈

※1) 因果律によるグリーン関数の境界条件は、妥当っぽく思われますが、人間が勝手にそうだろう、として経験的に与えている条件にすぎません。導くことは出来ない境界条件です。数学的には特解が何でも1つ与えられればいいので、時刻が未来であろうが特定の時間であろうが何でも構わないのです。
 今回の場合は解きたい微分方程式の右辺が定数(重力のみ)という非常に特別な場合ですが、一般的にはグリーン関数を用いた特解の被積分関数には、左辺の解である\(x(t)\)が含まれます。
この時、因果律を仮定しないと、『ある時刻の運動を計算するために、その時刻より未来の運動を知らなければならない』という理解しがたい状況が生まれます。言い換えれば、『現時刻の運動は、未来の振る舞いが決定されているので、どんな力を新たに加えようとも決して運動を変えることは出来ない』となってしまい、非物理的に感じてしまいます。
 もう一度言いますが、これは人間が「時間は必ず過去から未来に進む」という経験的に得られた事実を使用しているのであって、この根拠以外に因果律を採用する理由は存在しません。数学的に解こうとする際にはどんな条件でも構わないので、因果律を境界条件としても特解として正しい答えが得られるので良いのです。
 位置に対する微分方程式である場合は、例えば位置が無限大でゼロである条件だったり、ある点で位置と微分が決まっていたり、漸近で特定の形を持っていなければならない、という条件が当てはめられます。

(※2)これいいんですかね、使って。あんまり考えてはいません。もしも時刻が負の∞ではなかったらと思うと不安で仕方ありません。都合がよかったので使いましたが…ダメな気がしないでもないです。

Thanks

アイキャッチ画像のフォントはmagicringを使用しています。
配布ページ
http://inatsuka.com/extra/magicring/
配布元のサークル様
森の中の猫の小屋 http://inatsuka.com/