\(\def\bm{\boldsymbol}\)
\(\def\di{\displaystyle}\)
\(\def\ve{\varepsilon_0}\)
\(\def\dd#1#2{\dfrac{\partial #1}{\partial #2}}\)
\(\def\dda#1#2{\dfrac{\partial^2 #1}{\partial #2}}\)
光の速度で、ベクトル解析の公式を使ったので、まとめておきます。

１．はじめに

(1) ベクトルの内積
\(\quad\bm{a}\cdot\bm{b}=a_xb_x+a_yb_y+a_zb_z=\di\sum^3_ia_ib_i\)
(2) ベクトルの外積
\(\quad\bm{a}\times\bm{b}=(a_yb_z-a_zb_y,a_zb_x-a_xb_z,a_xb_y-a_yb_x)\)
\(\quad(\bm{a}\times\bm{b})_i=\di\sum_{j,k=1}^3\epsilon_{i,j,k}\,a_jb_k\)
\begin{align}\hspace{-20mm}\epsilon_{i,j,k}=1&\,:(i,j,k)=(1,2,3),(2,3,1),(3,1,2)\\
-1&\,:(i,j,k)=(1,3,2),(2,1,3),(3,2,1)\\
0&\,: other\end{align}
3. スカラー三重積
\(\quad\bm{A}\cdot(\bm{B}\times\bm{C})=\bm{B}\cdot(\bm{C}\times\bm{A})=\bm{C}\cdot(\bm{A}\times\bm{B})\)
\begin{align}
\bm{A}\cdot(\bm{B}\times\bm{C})&=A_x(\bm{B}\times\bm{C})_x+A_y(\bm{B}\times\bm{C})_y+A_z(\bm{B}\times\bm{C})_z\\[2pt]
&=A_x(B_yC_z-B_zC_y)+A_y(B_zC_x-B_xC_z)+A_z(B_xC_y-B_yC_x)\\[2pt]
&=A_xB_yC_z-A_xB_zC_y+A_yB_zC_x-A_yB_xC_z+A_zB_xC_y-A_zB_yC_x\\[2pt]
&=B_x(C_yA_z-C_zA_y)+B_y(C_zA_x-C_xA_z)+B_z(C_xA_y-C_yA_x)\\[2pt]
&=B_x(\bm{C}\times\bm{A})_x+B_y(\bm{C}\times\bm{A})_y+B_z(\bm{C}\times\bm{A})_z\\
&=\bm{B}\cdot(\bm{C}\times\bm{A})\\[-40pt]
\end{align}
\(\qquad\)また、\(\quad\bm{A}\cdot(\bm{B}\times\bm{C})=\left|\begin{array}{ccc}A_x & A_y & A_z\\ B_x & B_y & B_z\\C_x & C_y & C_z\end{array}\right|\qquad\)である。

4. ベクトル三重積
\(\quad\bm{A}\times(\bm{B}\times\bm{C})=(\bm{A}\cdot\bm{C})\bm{B}-(\bm{A}\cdot\bm{B})\bm{C}\)
\begin{align}
(\bm{A}\times(\bm{B}\times\bm{C}))_x&=A_y(\bm{B}\times\bm{C})_z-A_z(\bm{B}\times\bm{C})_y\\
&=A_y(B_xC_y-B_yC_x)-A_z(B_zC_x-B_xC_z)\\
&=(A_yC_y+A_zB_z)B_x-(A_yB_y+A_zB_z)C_x\\
&=(A_xC_x+A_yC_y+A_zB_z)B_x-(A_xB_x+A_yB_y+A_zB_z)C_x\\
&=(\bm{A}\cdot\bm{C})B_x-(\bm{A}\cdot\bm{B})C_x\\[6pt]
&\hspace{-12mm}\text{同様にして}\\[2pt]
(\bm{A}\times(\bm{B}\times\bm{C}))_y&=(\bm{A}\cdot\bm{C})B_y-(\bm{A}\cdot\bm{B})C_y\\[2pt]
(\bm{A}\times(\bm{B}\times\bm{C}))_z&=(\bm{A}\cdot\bm{C})B_z-(\bm{A}\cdot\bm{B})C_z\\[8pt]
&\hspace{-18mm}\therefore\quad \bm{A}\times(\bm{B}\times\bm{C})=(\bm{A}\cdot\bm{C})\bm{B}-(\bm{A}\cdot\bm{B})\bm{C}
\end{align}
5. 微分演算子
\(\quad\nabla=\left(\,\dd{}{x}\,,\,\dd{}{y}\,,\,\dd{}{z}\,\right)\)

２．勾配、発散、回転

勾配（ grad ）：grad\(\,\phi(\bm{r})=\nabla\,\phi(\bm{r})\equiv\left(\,\dd{}{x}f\,,\,\dd{}{y}f\,,\,\dd{}{z}f\,\right)\)

発散 （ div ）：div\(\,\bm{A}(\bm{r})=\nabla\cdot\bm{A}(\bm{r})\equiv\,\dd{}{x}A_x+\dd{}{y}A_y+\dd{}{z}A_z\)

回転 （ rot ）：rot\(\,\bm{A}(\bm{r})=\left(\dd{}{y}A_z-\dd{}{z}A_y\,,\,\dd{}{z}A_x-\dd{}{x}A_z\,,\,\dd{}{x}A_y-\dd{}{y}A_x\,\right)\)
\(\hspace{8mm}\)または \(\quad(\nabla\times\bm{A})_i=\di\sum_{j,k=1}^3\epsilon_{ijk}(\nabla)_j(\bm{A})_k=\di\sum_{j,k=1}^3\epsilon_{ijk}\,\partial_jA_k\)

３．計算例－１

\(\quad\bm{r}=(x,y,z)\quad r=|\bm{r}|\ne 0\quad\)とする。

(1)\(\,\nabla\,r\)
\begin{align}
\nabla\,r&=\nabla(\sqrt{x^2+y^2+z^2})\\[2pt]
&=\left(\dfrac{x}{\sqrt{x^2+y^2+z^2}}\,,\,\dfrac{y}{\sqrt{x^2+y^2+z^2}}\,,\,\dfrac{z}{\sqrt{x^2+y^2+z^2}}\,\right)\\[2pt]
&=\dfrac{\bm{r}}{r}
\end{align}
(2)\(\,\nabla\left(\dfrac{1}{r}\right)\)
\begin{align}
\nabla\left(\dfrac{1}{r}\right)&=\left(\dd{}{x}\dfrac{1}{\sqrt{x^2+y^2+z^2}}\,,\,\dd{}{y}\dfrac{1}{\sqrt{x^2+y^2+z^2}}\,,\,\dd{}{z}\dfrac{1}{\sqrt{x^2+y^2+z^2}}\,\right)\\[3pt]
&=\left(\,\dfrac{-x}{(x^2+y^2+z^2)^{3/2}}\,,\,\dfrac{-y}{(x^2+y^2+z^2)^{3/2}}\,,\,\dfrac{-z}{(x^2+y^2+z^2)^{3/2}}\,\right)\\[3pt]
&=-\dfrac{\bm{r}}{r^3}
\end{align}
(3)\(\,\nabla\cdot\bm{r}\)
\begin{equation}
\nabla\cdot\bm{r}=\dd{}{x}x+\dd{}{y}y+\dd{}{z}z=3
\end{equation}
(4)\(\,(\nabla\cdot\nabla)\,r\)
\begin{align}
(\nabla\cdot\nabla)r&=\left(\dda{}{x^2}+\dda{}{y^2}+\dda{}{z^2}\right)\sqrt{x^2+y^2+z^2}\\[2pt]
&=\dd{}{x}\dfrac{x}{\sqrt{x^2+y^2+z^2}}+\dd{}{y}\dfrac{y}{\sqrt{x^2+y^2+z^2}}+\dd{}{z}\dfrac{z}{\sqrt{x^2+y^2+z^2}}\\[2pt]
&=\dfrac{1}{\sqrt{x^2+y^2+z^2}}-\dfrac{x^2}{(x^2+y^2+z^2)^{3/2}}\\
& +\dfrac{1}{\sqrt{x^2+y^2+z^2}}-\dfrac{y^2}{(x^2+y^2+z^2)^{3/2}}\\
& +\dfrac{1}{\sqrt{x^2+y^2+z^2}}-\dfrac{z^2}{(x^2+y^2+z^2)^{3/2}}\\[2pt]
&=\dfrac{3}{r}-\dfrac{(x^2+y^2+z^2)}{r^3}=\dfrac{2}{r}
\end{align}
(5)\(\,\nabla\times\bm{r}\)
\begin{equation}
\nabla\times\bm{r}=\left(\dd{}{z}y-\dd{}{y}z\,,\,\dd{}{x}z-\dd{}{z}x\,,\,\dd{}{y}x-\dd{}{x}y\right)=\bm{0}
\end{equation}
(6)\(\,\nabla\times(\nabla\phi)\hspace{10mm}\mathrm{rot}\,\,\mathrm{grad}\,\,\phi\)
\begin{align}
\nabla\times(\nabla\phi)&=\left(\,\dd{}{y}(\dd{\phi}{z})-\dd{}{z}(\dd{\phi}{y})\,\,,\right.\\
&\hspace{20mm}\dd{}{z}(\dd{\phi}{x})-\dd{}{x}(\dd{\phi}{z})\,\,,\\
&\hspace{30mm}\left.\dd{}{x}(\dd{\phi}{y})-\dd{}{y}(\dd{\phi}{x})\,\right)\\
&=0\hspace{40mm}\because\,\,\dd{}{y}(\dd{\phi}{z})=\dd{}{z}(\dd{\phi}{y})
\end{align}
(7)\(\,\nabla\cdot(\nabla\times\bm{A})\hspace{10mm}\mathrm{div}\,\,\mathrm{rot}\,\,\bm{A}\)
\begin{align}
\nabla\cdot(\nabla\times\bm{A})&=\dd{}{x}(\dd{A_z}{y}-\dd{A_y}{z})+\dd{}{y}(\dd{A_x}{z}-\dd{A_z}{x})+\dd{}{z}(\dd{A_y}{x}-\dd{A_x}{y})\\[3pt]
&=\dfrac{\partial^2A_z}{\partial x\partial y}-\dfrac{\partial^2A_y}{\partial x\partial z}+\dfrac{\partial^2A_x}{\partial y\partial z}-\dfrac{\partial^2A_z}{\partial y\partial x}+\dfrac{\partial^2A_y}{\partial z\partial x}-\dfrac{\partial^2A_x}{\partial z\partial y}\\[3pt]
&=0\hspace{40mm}\because\,\,\,\dfrac{\partial^2A_z}{\partial x\partial y}=\dfrac{\partial^2A_z}{\partial y\partial x}
\end{align}
(8)\(\,\nabla\times(\nabla\times\bm{A})\hspace{10mm}\mathrm{rot}\,\,\mathrm{rot}\,\,\bm{A}\)
\begin{align}
(\nabla\times(\nabla\times\bm{A}))_x&=\dd{}{y}\left(\dd{A_y}{x}-\dd{A_x}{y}\right)-\dd{}{z}\left(\dd{A_x}{z}-\dd{A_z}{x}\right)\\[3pt]
&=\dfrac{\partial^2}{\partial y\partial x}A_y-\dda{}{y^2}A_x-\dda{}{z^2}A_x+\dfrac{\partial^2}{\partial z\partial x}A_z\\[3pt]
&=\dda{}{x^2}A_x+\dfrac{\partial^2}{\partial y\partial x}A_y+\dfrac{\partial^2}{\partial z\partial x}A_z-\left(\dda{}{x^2}+\dda{}{y^2}+\dda{}{z^2}\right)A_x\\[3pt]
&=\dd{}{x}\left(\dd{A_x}{x}+\dd{A_y}{y}+\dd{A_z}{z}\right)-\Delta A_x\\[3pt]
&=\dd{}{x}(\nabla\cdot\bm{A})-\Delta A_x\\[3pt]
&\hspace{-20mm}\text{同様に}\\
(\nabla\times(\nabla\times\bm{A}))_y&=\dd{}{y}(\nabla\cdot\bm{A})-\Delta A_y\\[3pt]
(\nabla\times(\nabla\times\bm{A}))_z&=\dd{}{z}(\nabla\cdot\bm{A})-\Delta A_z\\[3pt]
&\hspace{-20mm}\text{よって}\\[2pt]
\nabla\times(\nabla\times\bm{A})&=\nabla(\nabla\cdot\bm{A})-\Delta\bm{A}
\end{align}

\(\def\bm{\boldsymbol}\)
\(\def\di{\displaystyle}\)
\(\def\ve{\varepsilon_0}\)
\(\def\dd#1#2{\dfrac{\partial #1}{\partial #2}}\)
\(\def\dda#1#2{\dfrac{\partial^2 #1}{\partial #2}}\)
オーソドックスな形で、光（電磁波）の真空中の速度を求めてみる。

［１］マックスウェル方程式

光も電磁波なので、マックスウェルの方程式に従う。マックスウェル方程式は４つに分けられる。
(1)電場の発散
\[\mbox{div }\bm{E}=\nabla\cdot\bm{E}=\dfrac{\rho}{\ve}\tag{1.1}\]
電場の源が電荷であり、電場が電荷から放射状であることを表す。
ここで \(\,\rho\,\)は電荷密度で、\(\ve\,\)は真空の誘電率である。
(2)電場の回転
\[\mbox{rot }\bm{E}=\nabla\times\bm{E}=-\dd{\bm{B}}{t}\tag{1.2}\]
磁束密度の変化により電場が生じる　Faradayの法則=電磁誘導の法則
(3)磁場の発散
\[\mbox{div }\bm{B}=\nabla\cdot\bm{B}=0\quad\tag{1.3}\]
磁場には源がない
(4)磁場の回転
\[\mbox{rot }\bm{B}=\nabla\times\bm{B}=\mu_0(\bm{j}+\ve\dd{\bm{E}}{t})\tag{1.4}\]
電流及び電場の変化が磁場を生むことを表す。Ampére-Maxwellの法則
ここで\(\,\mu_0\,\)は真空の透磁率である。

［２］マックスウェル方程式の展開

真空中を伝わる電磁波について考えたいので, 電荷密度はいたるところで\(\,0\,\)であるとする. よって電流密度も\(\,0\,\)であるので、次の４つの式が得られる。
\begin{align}
\mbox{div }\bm{E}&=0 \tag{2.1}\\
\mbox{rot }\bm{E}&=-\dd{\bm{B}}{t} \tag{2.2}\\
\mbox{div }\bm{B}&=0 \tag{2.3}\\
\mbox{rot }\bm{B}&=\mu_0\ve\dd{\bm{E}}{t} \tag{2.4}
\end{align}
場の回転を改めて場とみなして、回転を調べる。式\(\,(2.2)\,\)を回転すると
\begin{align}
\nabla\times(\nabla\times\bm{E})&=-\dd{}{t}(\nabla\times\bm{B})\\
&=-\mu_0\dd{}{t}\left(\bm{j}+\ve\dd{\bm{E}}{t}\right)\\
&=-\mu_0\dd{\bm{j}}{t}-\mu_0\ve\dda{\bm{E}}{t^2}\tag{2.5}
\end{align}
となり、電流密度の他には電場のみの閉じた形となる。ここでは、ベクトル解析のベクトルの回転の回転からラプラシアンを導く、以下の公式を用いる。
\[\nabla\times(\nabla\times\bm{A})=\nabla(\nabla\cdot\bm{A})-\Delta\bm{A}\qquad\because\,\,\Delta\equiv\nabla\cdot\nabla\tag{2.6}\]
この公式を\(\,(2.5)\,\)式の右辺に適用する。
\begin{align}
\nabla\times(\nabla\times\bm{E})&=\nabla(\nabla\cdot\bm{E})-\Delta\bm{E}\\[3pt]
&=0-\Delta\bm{E}\qquad\,\because\,\,\nabla\cdot\bm{E}=0\quad(2.1)\,\,\text{より}\\[3pt]
&=-\Delta\bm{E}\tag{2.7}
\end{align}
\(\,(2.5)\,\)式と\(\,(2.6)\,\)をあわせると
\[\Delta\bm{E}-\mu_0\ve\dda{\bm{E}}{t^2}-\mu_0\dd{\bm{j}}{t}=0\tag{2.8}\]
ここの電流は\(\,0\,\)なので、第3項は消せる。ここで、\(x\,\)成分で考えると。
\[\dda{\bm{E}}{x^2}-\mu_0\ve\dda{\bm{E}}{t^2}=0\tag{2.10}\]
の形の波動方程式が得られる。

［3］波動方程式の一般解

この方程式の解は次のような形で表される。例えば\(\,\bm{E}\,\)の振幅を\(\,z\,\)軸にとると、解は
\[E_z=f(x-vt)+g(x+vt)\tag{3.1}\]
の形で表せる。ここで、\(f\,\)と\(\,g\,\)は任意の関数であり、\(f(x − vt)\,\,\)は\(\,x\,\)軸の正の向きに速さ\(\,v\,\)で
進む波動を表し、\(g(x + vt)\,\,\)は\(\,x\,\)軸の負の向きに速さ\(\,v\,\)で進む波動を表す。
ここで、式\(\,(2.2)\,\)を成分に分解する。
\begin{align}
\dd{E_z}{y}-\dd{E_y}{z}&=-\dd{B_x}{t}\\
\dd{E_x}{z}-\dd{E_z}{x}&=-\dd{B_y}{t}\\
\dd{E_y}{x}-\dd{E_x}{y}&=-\dd{B_z}{t}
\end{align}
ここで、\(E_x=0\,\,,\,\,E_y=0\,\)なので、次のような結果となる
\begin{align}
\dd{B_x}{t}&=0\\
\dd{B_y}{t}&=f^{\prime}(x-vt)+g^{\prime}(x+vt) \tag{3.2}\\
\dd{B_z}{t}&=0
\end{align}
式\(\,(3.2)\,\)を\(\,t\,\)で積分すると
\[B_y=-\dfrac{1}{v}f(x-vt)+\dfrac{1}{v}g(x+vt)\tag{3.3}\]
となる。つまり、\(E_z\,\)と\(\,B_y\,\)は、同じ関数\(\,f\,\)と\(\,g\,\)で表され、両者が互いに組み合って、離れることなく、同じ速さ\(\,v\,\)を持つ波動となって\(\,x\,\)軸を伝搬する。

電場と磁場の進行波

［4］電磁波の伝搬する速さ

波動方程式の一般解と波動方程式を再掲する。
\begin{align}
&E_z=f(x-vt)+g(x+vt)\tag{3.1}\\
&\dda{\bm{E}}{x^2}-\mu_0\ve\dda{\bm{E}}{t^2}=0\tag{2.10}
\end{align}
この\(\,(3.1)\,\)式を\(\,(2.10)\,\)に代入する。このとき、\(x-vt=p\,\,,\,\,x+vt=q\,\,\)とする。
\[\dd{E_z}{x}=\dd{f}{x}+\dd{g}{x}=\dd{f}{p}\dd{p}{x}+\dd{g}{q}\dd{q}{x}=\dd{f}{p}+\dd{g}{q}\quad\text{なので}\]
\[\dda{E_z}{x^2}=\dda{f}{p^2}+\dda{g}{q^2}\qquad\text{となる。次に時間で微分する。}\]
\[\dd{E_z}{t}=\dd{f}{t}+\dd{g}{t}=\dd{f}{p}\dd{p}{t}+\dd{g}{q}\dd{q}{t}=-v\dd{f}{p}+v\dd{g}{q}\quad\text{なので}\]
\[\dda{E_z}{t^2}=v^2\dda{f}{p^2}+v^2\dda{g}{q^2}=v^2\dda{E_z}{x^2}\qquad\text{となるので}\quad (2.10)\,\text{より}\]
\[v^2=\dfrac{1}{\mu_0\ve}\qquad\text{よって}\qquad v=\dfrac{1}{\sqrt{\mu_0\ve}}\tag{4.1}\]
真空の誘電率\(\,\ve\,\)と透磁率\(\,\mu_0\,\)の値

\(\qquad\ve=\dfrac{10^7}{4\pi c^2}\quad[\,\mathrm{C}\cdot\mathrm{N}^{-1}\cdot\mathrm{m}^{-2}\,\,]\hspace{20mm}\mu_0=4\pi\times 10^{-7}\quad[\,\,\mathrm{N}\cdot\mathrm{A}^{-2}\,\,]\quad\)を代入すると

\(\qquad\quad v=\dfrac{1}{\sqrt{\mu_0\ve}}=c\qquad\)となり、電磁波の速度が光速と一致する。

まあ、誘電率の定義からして、光速を使っているので当然の結果である。

1.経緯

放送大学の「物理と化学のための数学演習」という授業のレポートの課題のひとつに「二重振り子」があったので、忘れないようにここにメモっておきます。

2.レポート課題

天井から長さ \(\ell\) の紐で吊るされた質量 \(M\) の錘に、もう一つの同じ長さ \(\ell\) の紐で同じ質量の錘が垂れ下がった、二重振り子の微小振動を考える。この問題は自由度２の連成振動となっている。以下の問いに答えよ。
\(\textbf{[1]}\) それぞれの紐が鉛直線となす角度を \(\theta_1\,,\,\theta_2\,\)とし、近似式、\(\sin\theta\simeq\theta\,\)を用いて、２つの錘の運動方程式を導き、この連成振動の基準座標と基準振動数を求めよ。紐の質量は小さいとして無視する。
\(\textbf{[2]}\) ２つの紐を変形しない同じ質量 \(m\) の棒で置き換えるとどうなるか。棒の重心はその真ん中にあるとする。

※　2. の問題は、私のやり方が拙かったせいか、解析解が得られませんでした。そこで、1. だけをここにメモります。

3.方程式の導出

最初のオモリの座標を \((x_1,y_1)\,\) , そこから垂れ下がったオモリの座標を \((x_2,y_2)\,\) とする。

\(x_1=\ell\sin\theta_1\qquad y_1=-\ell\cos\theta_1\)
\(x_2=x_1+\ell\sin\theta_2=\ell(\sin\theta_1+\sin\theta_2)\)
\(y_2=y_1-\ell\cos\theta_2=-\ell(\cos\theta_1+\cos\theta_2)\)
\(\dot{x}_1=\ell\cos\theta_1\dot{\theta}_1\simeq\ell\dot{\theta}_1\)
\(|\dot{x}_1|^2=\ell^2\dot{\theta}_1^2\)
\(\dfrac{1}{2}|\dot{x}_1|^2=\dfrac{\ell^2}{2}\dot{\theta}_1^2\)
\(\dot{x}_2=\ell(\cos\theta_1\dot{\theta}_1+\cos\theta_2\dot{\theta}_2)\simeq\ell(\dot{\theta}_1+\dot{\theta}_2)\)
\(\dfrac{1}{2}|\dot{x}_2|^2=\dfrac{\ell^2}{2}(\dot{\theta}_1+\dot{\theta}_2)^2\)

\(\theta_1\,,\,\theta_2\,\)が小さいとするので、運動エネルギーは \(\quad T=\dfrac{M\ell^2}{2}\dot{\theta}_1^2+\dfrac{M\ell^2}{2}(\dot{\theta}_1+\dot{\theta}_2)^2\)

ポテンシャル・エネルギーは\(\quad U=-Mg\ell(\cos\theta_1+\cos\theta_1+\cos\theta_2)\)

よって、ラグラジアン\(\,L\,\)は

\(L=T-U=\dfrac{M\ell^2}{2}(2\dot{\theta}_1^2+2\dot{\theta}_1\dot{\theta}_2+\dot{\theta}^2)+Mg\ell(2\cos\theta_1+\cos\theta_2)\)

\(\dfrac{\partial L}{\partial\theta_1}\,=-2Mg\ell\sin\theta_1\simeq -2Mg\ell\theta_1\hspace{12mm}\dfrac{\partial L}{\partial\dot{\theta}_1}=M\ell^2(2\dot{\theta}_1+\dot{\theta}_2)\)

\(\dfrac{\partial L}{\partial\theta_2}\,=-Mg\ell\sin\theta_2\simeq -Mg\ell\theta_2\hspace{16mm}\dfrac{\partial L}{\partial\dot{\theta}_2}=M\ell^2(\dot{\theta}_1+\dot{\theta}_2)\)

オイラー・ラグランジュ方程式は\(\quad\dfrac{d}{dt}\left(\dfrac{\partial L}{\partial\dot{q}_1}\right)=\dfrac{\partial L}{\partial q}\,\,\)なので

\(2\ddot{\theta}_1+\ddot{\theta}=-2\dfrac{g}{\ell}\theta_1\qquad \ddot{\theta}_1+\ddot{\theta}_2=-\dfrac{g}{\ell}\theta_2\quad\)の関係が得られる。
これより
\(\left(\begin{array}{@{\,}c@{\,}}\ddot{\theta}_1\\[2pt] \ddot{\theta}_2\end{array}\right)=\dfrac{g}{\ell}\left(\begin{array}{@{\,}rr@{\,}}-2&1\\[2pt]2&-2\end{array}\right)\left(\begin{array}{@{\,}c@{\,}}\theta_1\\[2pt] \theta_2\end{array}\right)\qquad\)の連立微分方程式が導かれる。
\(\hspace{28mm}\)行列\(\,\,K\,\)

4.方程式の解法

この行列 \(K\) を対角化することで、方程式を解く。まず \(K\) の固有値と固有ベクトルを求める。
固有値 \(\lambda\) は、固有ベクトル \(\textbf{x}\) と \(K\textbf{x}=\lambda\textbf{x}\) の関係にあるので、\(\det{(K-\lambda I)}=0\) より
\(\lambda^2+4\lambda+2=0\qquad \lambda=-2\pm\sqrt{2}\qquad\)が得られる。

\(\textbf{(4-1)} \,\,\lambda_1=-2+\sqrt{2}\) のとき

\(\quad(K-\lambda_1 I)=\left(\begin{array}{@{\,}rr@{\,}}-\sqrt{2}&1\\[2pt]2&-\sqrt{2}\end{array}\right)\left(\begin{array}{@{\,}c@{\,}}x_1\\[2pt] x_2\end{array}\right)=0\qquad\) より

\(\quad x_1=\dfrac{1}{\sqrt{3}}\qquad x_2=\dfrac{\sqrt{2}}{\sqrt{3}}\qquad\quad \textbf{x}_1=\left(\begin{array}{@{\,}r@{\,}}\frac{1}{\sqrt{3}}\\[2pt]\frac{\sqrt{2}}{\sqrt{3}}\end{array}\right)\)

\(\textbf{(4-2)} \,\,\lambda_2=-2-\sqrt{2}\) のとき

\(\quad(K-\lambda_2 I)=\left(\begin{array}{@{\,}rr@{\,}}\sqrt{2}&1\\[2pt]2&\sqrt{2}\end{array}\right)\left(\begin{array}{@{\,}c@{\,}}x_1\\[2pt] x_2\end{array}\right)=0\qquad\) より

\(\quad x_1=\dfrac{1}{\sqrt{3}}\qquad x_2=-\dfrac{\sqrt{2}}{\sqrt{3}}\qquad\quad \textbf{x}_2=\left(\begin{array}{@{\,}r@{\,}}\frac{1}{\sqrt{3}}\\[2pt]-\frac{\sqrt{2}}{\sqrt{3}}\end{array}\right)\)

5.正則行列の逆行列

固有ベクトルを縦に並べた行列 \(P\) を用いて \(P^{-1}KP\) として \(K\) を対角化する。

\(\quad P=\left(\begin{array}{@{\,}rr@{\,}}\frac{1}{\sqrt{3}}&\frac{1}{\sqrt{3}}\\[3pt]\frac{\sqrt{2}}{\sqrt{3}}&-\frac{\sqrt{2}}{\sqrt{3}}\end{array}\right)\qquad\quad P^{-1}=\left(\begin{array}{@{\,}rr@{\,}}\frac{\sqrt{3}}{2}&\frac{\sqrt{3}}{2\sqrt{2}}\\[3pt]\frac{\sqrt{3}}{2}&-\frac{\sqrt{3}}{2\sqrt{2}}\end{array}\right)\)

この \(P^{-1}\) を用いて、新しい変数 \(\Theta_1\,\,,\,\,\Theta_2\,\,\)を定義する。

\(\quad\left(\begin{array}{@{\,}c@{\,}}\Theta_1\\[2pt] \Theta_2\end{array}\right)=P^{-1}\left(\begin{array}{@{\,}c@{\,}}\theta_1\\[2pt]\theta_2\end{array}\right)=\left(\begin{array}{@{\,}rr@{\,}}\frac{\sqrt{3}}{2}&\frac{\sqrt{3}}{2\sqrt{2}}\\[3pt]\frac{\sqrt{3}}{2}&-\frac{\sqrt{3}}{2\sqrt{2}}\end{array}\right)\left(\begin{array}{@{\,}c@{\,}}\theta_1\\[2pt]\theta_2\end{array}\right)\)

元の連立微分方程式は、変数　\(\Theta_1\,\,,\,\,\Theta_2$\,\,\)に対する以下の微分方程式と等価である。

\(\quad\dfrac{d}{dt}\left(\begin{array}{@{\,}c@{\,}}\Theta_1\\[2pt] \Theta_2\end{array}\right)=\dfrac{g}{\ell}\left(\begin{array}{@{\,}rr@{\,}}-2+\sqrt{2}&0\\[3pt]0&-2-\sqrt{2}\end{array}\right)\left(\begin{array}{@{\,}c@{\,}}\Theta_1\\[2pt] \Theta_2\end{array}\right)\)

最終的に、運動は２つの固有振動の足し合わせで表され、それぞれの固有振動数 \(\omega_1\,\,,\,\,\omega_2\,\,\)は、

\(\quad\omega_1=\sqrt{2+\sqrt{2}}\sqrt{\dfrac{g}{\ell}}\quad,\quad \omega_1=\sqrt{2-\sqrt{2}}\sqrt{\dfrac{g}{\ell}}\quad\) となる。

放送大学のオンライン授業で「力学演習」を履修している。そこでの練習問題で一次元の弾性衝突問題があり、すっかり忘れているので再確認のため、計算してみました。

図のような２球で考える。向きは\(\,x\,\)軸方向に正として、初速度をそれぞれ\(\,v_1\,\,,\,v_2\,\)とする。衝突後の速度をそれぞれ\(\,v_3\,\,,\,v_4\,\)とする。弾性衝突なので、運動量とエネルギーが保存する。

・運動量保存
\(\qquad m_A\,v_1+m_B\,v_2=m_A\,v_3+m_B\,v_4\quad\cdots\,\)(1)
・エネルギー保存
\(\qquad \dfrac{1}{2}m_A\,v_1^2+\dfrac{1}{2}m_B\,v_2^2=\dfrac{1}{2}m_A\,v_3^2+\dfrac{1}{2}m_B\,v_4^2\quad\cdots\,\)(2)

式(1)から、\(\quad v_4=\dfrac{m_A(v_1-v_3)+m_Bv_2}{m_B}\quad\cdots\,\)(3)
式(2)から、\(\quad v_4^2=\dfrac{m_A(v_1^2-v_3^2)+m_Bv_2^2}{m_B}\quad\cdots\,\)(4)
式(3)を２乗すると、
\(\quad v_4^2=\dfrac{m_A^2v_1^2-2m_A^2v_1v_3+m_A^2v_3^2+m_B^2v_2^2+2m_Am_Bv_1v_2-2m_Am_Bv_2v_3}{m_B^2}\quad\cdots\,\)(5)
式 (4)\(\,=\,\)(5) なので、
\(\quad m_Am_Bv_1^2-m_Am_Bv_3^2+m_B^2v_2^2\)
\(\qquad\qquad =m_A^2v_1^2-2m_A^2v_1v_3+m_A^2v_3^2+m_B^2v_2^2+2m_Am_Bv_1v_2-2m_Am_Bv_2v_3\quad\cdots\,\)(6)

\(\quad m_A^2v_1^2-2m_Av_1v_3+m_A^2v_3^2+2m_Am_Bv_1v_2-2m_Am_Bv_2v_3-m_Am_Bv_1^2+m_Am_Bv_3^2=0\,\cdots\,\)(7)

\(\quad(m_A+m_B)v_3^2-2(m_Av_1+m_Bv_2)v_3+m_Av_1^2+2m_Bv_1v_2-m_Bv_1^2=0\quad\cdots\,\)(8)

この式(8) を\(\,v_3\,\)の２次方程式として解くと
\(\quad v_3=\dfrac{m_Av_1+m_Bv_2\pm\sqrt{m_B^2(v_1-v_2)^2}}{m_A+m_B}=\dfrac{m_Av1+m_Bv_2\pm m_B|v_1-v_2|}{m_A+m_B}\quad\cdots\,\)(9)

ここで、球 A の速度が B の速度より大きいケースで考えてみる。 \(\,v_1\ge v_2\)
\(\,\pm\,\)の\(\,+\,\)では、\(\,v_3=v_1\,\)となり、すり抜ける形なので\(\,-\,\)を採用する。

\(\quad v_3=\dfrac{m_Av_1-m_Bv_1+2m_Bv_2}{m_A+m_B}\quad\cdots\,\)(10)

ここで、両方の球の質量が同じ場合、\(\,v_3=v_2\,\)となり、A と B の運動量が入れ替わることになる。

◇ 大小球の正面衝突

大小の球を正面衝突させた場合を考える。
Ａは\(\,m_A=3m\,\)とし、Ｂは\(\,m_B=m\,\)とする。
速度はともに\(\,v_1=v_2=v_0\,\)(符号は正負となる)として、正面からの衝突とする。

式(10)と式(1)により

\(\quad v_3=\dfrac{3mv_0-mv_0-2mv_0}{3m+m}=0\quad v_4=\dfrac{3mv_0-mv_0}{m}=2v_0\,\)となる。

結局、Ａは停止しＢは速度\(\,2v_0\,\)で反対側に跳ね返ることになる。

◇ ３球での衝突問題

『中川のビジュアル物理学教室』のページの中の「解説：一次元衝突」に次のような３球の衝突が紹介されていました。(画像は中川氏のページのものを使わせてもらいました)

静止している球Ｂ，Ｃに，その２倍の質量をもつ球Ａを\(\,v_0\,\)で弾性衝突したらどうなるだろうか？

この問題は，過去に青学大の入試問題として出題されたことがあるという事でした。

まず、Ａが静止しているＢに衝突する時点を考える。式(10) に\(\,m_A=2m\,,\,m_B=m\,,\,v_1=v_0\,,\,v_2=0\,\)を代入すると、
\(\quad v_A=v_3=\dfrac{2mv_0-mv_0}{3m}=\dfrac{1}{3}v_0\quad\cdots\,\)(11)
静止していたＢの速さ\(\,v_B\,\)は、式(11)の値を式(1)に代入すると\(\,v_B=v_4=\dfrac{4}{3}v_0\,\)となる。
Ｂがこの速度\(\,v_B\,\)で、静止しているＣに衝突すると
ＢとＣは同質量なので、式(11)より、Ｂは停止してＣは速度\(\,v_C=v_B\,\)で動き出すことになる。
停止した(１回動いて停止した形になる)ＢにＡが\(\,v_A=\dfrac{1}{3}v_0\,\)で衝突することになる。
同じく、式(11)に代入すると、再度のＡの速度\(\,v_A’\,\)は
\(\quad v_A’=\dfrac{2m\frac{1}{3}v_0-m\frac{1}{3}v_0}{3m}=\dfrac{1}{9}v_0\quad\cdots\,\)(12)
これを式(1)に代入すると\(\quad2m\dfrac{1}{3}=2m\dfrac{1}{9}v_0+mv_B’\,\)より
Ａは\(\,\dfrac{1}{9}v_0\,\)で動き出し、Ｂは\(\,\dfrac{4}{9}v_0\,\)で動き出す。
Ｃは先程、Ｂの衝突で\(\dfrac{4}{3}v_0\,\)で動き出しているので、３球はバラバラになることになる。

1.経緯

私が通学（通信制なので、通ってはいないが）している放送大学のゼミで、自宅での実験を動画撮影して皆で共有するという企画があったので、振り子の実験をやってみました。
振り子は小学\(5\)年生用の実験セットをAmazonから仕入れました。
（「ふりこのはたらき」 \(730\)円）

2.実験

ひもの長さは\(30\)cmとして、振れ角測定用に角度を書いた紙を付けました。大きく振った位置から重りを降ろした。
(\(2020.09.05\) 実施)
撮影は iPhone をそのままの設定で(\(30\)f)三脚無しで行いました。時間は約\(2\)分間でした。
動画は DropBox 経由で、Windows マシンに落とした。再生で簡単にコマ落としが出来、時間の表示が \(1/100\) 秒単位で出てくるので MovieMaker を使用した。（最近は MovieMaker はサポートされていないようですね）

撮影した動画の抜粋

3.実験結果

動画をMovieMaker でコマ落とし再生して、最大振幅の角度と時間を記録したのが下図です。

動画が\(30\)fなので、時間刻みが\(1/30\)sec単位となり、グラフは階段状になるが、振れ角が大きいと周期が長いという結果となった。

4.検討1：振幅の小さい場合

糸の張力\(\,T\,\)、重力\(\,W\,\)の合力\(\,f^{\prime\prime}\)が重りに掛かる。
\(\quad W=mg\)
\(\quad f^{\prime\prime}=mg\sin\theta\)
円周方向の重りの動く長さを\(\,s\,\)とすると
\(\quad m\dfrac{d^2s}{dt^2}=-mg\sin\theta\)
\(\quad s=l\theta\,\)なので\(\quad\dfrac{d^2\theta}{dt^2}=-\dfrac{g}{l}\sin\theta\)
振幅が小さいときは
\(\quad\theta\approx\sin\theta\quad\)とすることが出来るので
\(\quad\dfrac{d^2\theta}{dt^2}=-\dfrac{g}{l}\theta\quad\)とすることが出来、次のような一般解が得られる。

\(\quad\theta=A\sin\sqrt{\dfrac{g}{l}}t+B\cos\sqrt{\dfrac{g}{l}}t\)
ここで初期条件として\(\quad t=0\quad d\theta/dt=0\quad\theta=\theta_0\quad\)とすると
\(\quad\theta=\theta_0\cos\sqrt{\dfrac{g}{l}}t\quad\)が得られる。
周波数\(\,f\,\)は\(\quad f=\dfrac{1}{2\pi}\sqrt{\dfrac{g}{l}}\quad\)となる。
周期\(\,T\,\)は\(\quad T=\dfrac{1}{f}=2\pi\sqrt{\dfrac{l}{g}}\quad\)である。

今回は紐の長さを\(\,30\,\)cm としたので、

\(\quad T=2\times 3.141593\times\sqrt{\dfrac{0.3}{9.80655}}=1.099..\,\)(s) となった。結果のグラフを見ると、振れ角\(20\,\)度以下ではほぼこの周期となっているようである。

5.検討2：振幅が大きい場合

エネルギー保存則を用いて、振幅の大きい振り子の運動を計算する。振り子の最下点を原点とし、振りだし高さを\(\,h_0\,\)、\(\Delta h\,\)だけ降りた点の円周方向の速度を\(\,v\,\)とすると
\(\quad\dfrac{1}{2}mv^2=mgl(\cos\theta-\cos\theta_0)\quad\)なので
\(\quad\therefore v=\sqrt{2gl(\cos\theta-\cos\theta_0)}\)
ここで、\(\,v=\dfrac{ds}{dt}=l\dfrac{d\theta}{dt}\quad\)なので
\(\quad l\dfrac{d\theta}{dt}=\sqrt{2gl(\cos\theta-\cos\theta_0)}\quad\)より
\(\quad dt=\sqrt{\dfrac{l}{2g}}\dfrac{d\theta}{\sqrt{\cos\theta-\cos\theta_0}}\)

これを積分すると周期を求めることが出来る。
\(\,t\,\)の積分範囲を\(\,4\,\)分の\(\,1\,\)周期とすると、\(\theta=0\,\to\,\theta_0\)
\(\displaystyle\int^{\frac{T}{4}}_0dt=\dfrac{T}{4}=\sqrt{\dfrac{l}{2g}}\displaystyle\int_0^{\theta_0}\dfrac{d\theta}{\sqrt{\cos\theta-\cos\theta_0}}\)

ここからの積分は振り子の計算その２を参考にしてください。

プロセスとしては\(\quad\sin\phi=\sin(\theta/2)/\sin(\theta_0/2)\quad\)とおいて、いくつかの計算を経て以下の積分となる。
\(\quad T=4\sqrt{\dfrac{l}{g}}\displaystyle\int_0^{\frac{\pi}{2}}\dfrac{d\phi}{\sqrt{1-\sin^2(\theta_0/2)\sin^2\phi}}\)
この積分は楕円積分と呼ばれる。

5.1 楕円積分の計算

楕円積分を参照のこと。

第１種の楕円積分は以下のような式となる。

\(\quad I=\displaystyle\int_0^{\frac{\pi}{2}}\dfrac{d\phi}{\sqrt{1-a^2\sin^2\phi}}=\dfrac{\pi}{2}\displaystyle\sum_{n=0}^{\infty}\left\{\dfrac{2n-1)!!}{2n!!}\right\}^2a^{2n}\)
ここで\(\quad n!!\quad\)は二重階乗である。

これより、\(\,T\quad\)は以下のようになる。
\(\quad T=4\sqrt{\dfrac{l}{g}}\displaystyle\int_0^{\frac{\pi}{2}}\frac{d\phi}{\sqrt{1-\sin^2(\theta_0/2)\sin^2\phi}}\)

\(\qquad=2\pi\sqrt{\dfrac{l}{g}}\displaystyle\sum_{n=0}^{\infty}\left\{\dfrac{(2n-1)!!}{2n!!}\right\}^2\left(\sin\dfrac{\theta_0}{2}\right)^{2n}\)

具体的な計算は\(\quad b=\sin\dfrac{\theta_0}{2}\quad\)とおいて、振り子の長さ\(\,l=0.3\,\)m を使うと

\(\quad　T=2\pi\sqrt{\dfrac{l}{g}}\left\{1+\left(\dfrac{1}{2}\right)^2b^2+\left(\dfrac{1}{2}\dfrac{3}{4}\right)^2b^4+\left(\dfrac{1}{2}\dfrac{3}{4}\dfrac{5}{6}\right)^2b^6\right.\)

\(\qquad\qquad\qquad \left.+\left(\dfrac{1}{2}\dfrac{3}{4}\dfrac{5}{6}\dfrac{7}{8}\right)^2b^8+\left(\dfrac{1}{2}\dfrac{3}{4}\dfrac{5}{6}\dfrac{7}{8}\dfrac{9}{10}\right)^2b^{10}+\cdots\,\right\}\)

\(\qquad=1.099\times\left\{1+\dfrac{1}{4}b^2+\dfrac{9}{64}b^4+\dfrac{25}{256}b^6+\dfrac{1225}{16384}b^8\right.\)

\(\qquad\qquad\qquad\left.+\dfrac{3969}{65536}b^{10}+\dfrac{53361}{1048576}b^{12}+\cdots\,\right\}\)

この値と実験結果を併せて表示する。

測定データが \(\quad 1/30\,\)sec 単位であるが、まあ妥当な結果となったようである。

積分その１

【１】\(\,\,\displaystyle\int_0^{\infty}xe^{-ax}dx\quad\)の計算

\(\quad\displaystyle\int_0^{\infty}xe^{-ax}dx\quad\cdots\,\)(1)

まず、\(\,\displaystyle\int xe^{-ax}dx\,\,\)の不定積分から計算する。

部分積分で考えると

\(\quad\displaystyle\int xe^{-ax}dx=x(-ae^{-ax})+\dfrac{1}{a}\displaystyle\int e^{-ax}dx\)
\(\qquad=-axe^{-ax}+\dfrac{1}{a}\displaystyle\int e^{-ax}dx\quad\cdots\,\)(2)

(1)式の定積分を考える時、まず次の極限を計算する。

\(\quad\displaystyle\lim_{x\to\infty}\dfrac{x}{e^x}\quad\cdots\,\)(3)

\(\quad x\ge 0\,\,\)のとき\(\,\,e^x\ge\dfrac{x^2}{2}\,\,\)を考える。

\(\quad f(x)=e^x-\dfrac{x^2}{2}\quad\cdots\,\)(4)

とすると、\(x\ge 0\quad\)のとき、\(f'(x)\gt 0\,\,,\quad f^{\prime\prime}\gt 0\quad\)なので、\(f(x)\gt 0\quad\)である。よって

\(\,\,e^x\ge\dfrac{x^2}{2}\,\,\)が成り立つ。そこで

\(\quad 0\ge \dfrac{x}{e^x}\ge \dfrac{x}{\frac{x^2}{2}}\quad\)つまり\(\quad 0\ge \dfrac{x}{e^x}\ge \dfrac{2}{x}\quad\)となり、

\(\quad\displaystyle\lim_{x\to\infty}\dfrac{2}{x}=0\quad\)なので\(\quad\displaystyle\lim_{x\to\infty}\dfrac{x}{e^x}=0\quad\cdots\,\)(5)

となる。

この結果 (2)式を\(\,0\,\)から\(\,\infty\,\)まで積分すると

\(\quad\displaystyle\int_0^{\infty}xe^{-ax}dx=\left[-axe^{-ax}\right]_0^{\infty}+\dfrac{1}{a}\displaystyle\int_0^{\infty}e^{-ax}dx\)

\(\qquad =0-0+\dfrac{1}{a}\left[-\dfrac{1}{a}e^{-ax}\right]_0^{\infty}=\dfrac{1}{a^2}\quad\cdots\,\)(6)

グラフは\(\,y=xe^{-x}\,\)である。\(x\,\to\,\infty\,\)まで積分すると斜線部分の面積は\(\,1\,\)となります。

1-1. クラウジウスのビリアル定理

ビリアル定理は １９世紀に Clausius によって考案された。
「系の平均活力は、その（平均）ビリアル（の大きさ）に等しい。」
ここで登場する「活力」（”vis viva”）は、今日の運動エネルギーに相当します。 そもそも「活力」は古典力学の草創期にライプニッツが導入した量で、今日の運動エネルギーの２倍に相当する量 \(\,mv^2\,\)ですが、クラウジウスは \(\,mv^2/2\,\)とちょうど運動エネルギーと同じ量として用いています。ビリアル virial はラテン語のvis（「力」）からクラウジウスが作った造語です。ここでは、（力）×（位置ベクトル）をビリアルと呼ぶことにします。

距離の逆２乗則に従う重力やクーロン力の中心力で相互作用しあっている多体系では、長時間平均した運動エネルギー\(\,\langle\,T\,\rangle\,\)と平均の全ポテンシャルエネルギー\(\,\langle\,V\,\rangle\,\)は次の関係式を満たす。

\(\quad 2\langle\,T\,\rangle+\langle\,V\,\rangle=0\quad\cdots\,\)(1)

1-2. ビリアル定理の古典力学的証明

ここで位置ベクトル\(\,\vec{r}\,\)と運動量\(\,\vec{p}\,\)の内積の総和を以下のように考える。

\(\quad G=\displaystyle\sum_i\vec{p}_i\,\cdot\,\vec{r}_i\quad\cdot\,\)(2)

(2)式を時間\(\,t\,\)で微分する。

\(\quad\dfrac{dG}{dt}=\displaystyle\sum_i\vec{p}_i\cdot\dfrac{d\vec{r}_i}{dt}+\displaystyle\sum_i\dfrac{d\vec{p}_i}{dt}\cdot\vec{r}_i\)
\(\qquad\quad=\displaystyle\sum_i\,m_i(\vec{v}_i)^2+\displaystyle\sum_i\vec{F}_i\cdot\vec{r}_i\)
\(\qquad\quad=\displaystyle\sum_i2\times\dfrac{1}{2}m_i(\vec{v}_i)^2+\displaystyle\sum_i\vec{F}_i\cdot\vec{r}_i\)
\(\qquad\quad=2T+\displaystyle\sum_i\vec{F}_i\cdot\vec{r}_i\quad\cdots\,\)(3)

ここでは、次の関係式を使用している。

\(\quad \vec{p}_i=m_i\vec{v}_i\,,\qquad\vec{v}_i=\dfrac{d\vec{r}_i}{dt}\,,\qquad\vec{F}_i=\dfrac{d\vec{p}_i}{dt}\quad\cdots\,\)(4)

(3)式の両辺を\(\,0\,\)から時間\(\,t\,\)の範囲で積分して\(\,t\,\)で割り、\(t\to\infty\,\)の極限をとって長時間平均する。すると粒子が動き得る範囲は有限なので\(\,G\,\)も有限だから、左辺は\(\,0\,\)に収束する。

\(\quad\displaystyle\lim_{t\to\infty}\dfrac{1}{t}\displaystyle\int_0^t\dfrac{dG}{dt}=\displaystyle\lim_{t\to\infty}\dfrac{G(t)-G(0)}{t}=0\quad\cdots\,\)(5)

したがって、

\(\quad 0=2\langle\,T\,\rangle+\left\langle\displaystyle\sum_i\vec{F}_i\cdot\vec{r}_i\right\rangle\quad\cdots\,\)(6)

つまり、ビリアル定理を得る。

\(\quad\langle\,T\,\rangle=-\dfrac{1}{2}\left\langle\displaystyle\sum_i\vec{F}_i\cdot\vec{r}_i\right\rangle\quad\cdots\,\)(7)

次にポテンシャルエネルギー\(\,V\,\)が中心力ポテンシャルで、粒子間の距離に反比例する形で、系のポテンシャル\(\,V\,\)が各粒子対の相互作用の和によって書き表される場合、以下のように表される。

\(\quad V(r_1,\cdots,r_N)=\displaystyle\sum_{i\lt j}\dfrac{a_{ij}}{|\vec{r}_i-\vec{r}_j|}\quad\cdots\,\,\)(8)

粒子\(\,i\,\)に作用する力の合計\(\,\vec{F}_i\,\)は、ポテンシャルを位置座標で微分することで以下のように表すことが出来る。

\(\quad\vec{F}_i=-\nabla_{r_i}V=\displaystyle\sum_{j\ne i}\dfrac{a_{ij}(\vec{r}_i-\vec{r}_j)}{|\vec{r}_i-\vec{r}_j|^3}=\displaystyle\sum_{j\ne i}\vec{F}_{ij}\quad\cdots\,\)(9)

ここで、\(\,\vec{F}_{ij}\,\)は、粒子\(\,j\,\)から粒子\(\,i\,\)に働く力である。

\(\quad\vec{F}_{ij}=a_{ij}\dfrac{\vec{r}_i-\vec{r}_j}{|\vec{r}_i-\vec{r}_j|^3}\quad\cdots\,\)(10)

この\(\,\vec{F}_{ij}\,\)を用いると、ビリアル（力と位置ベクトルの内積）は次のように表せる。

\(\quad\displaystyle\sum_i\vec{F}_i\cdot\vec{r}_i=\displaystyle\sum_{i,j(i\lt j)}\vec{F}_{ij}\cdot\vec{r}_i+\displaystyle\sum_{i,j(j\lt i)}\vec{F}_{ij}\cdot\vec{r}_i\)
\(\qquad\qquad\quad=\displaystyle\sum_{i,j(i\lt j)}\vec{F}_{ij}\cdot\vec{r}_i+\displaystyle\sum_{i,j(i\lt j)}\vec{F}_{ji}\cdot\vec{r}_j\qquad\because\,i\leftrightarrow j\)
\(\qquad\qquad\quad=\displaystyle\sum_{i,j(i\lt j)}\vec{F}_{ij}\cdot(\vec{r}_i-\vec{r}_j)\quad\because\,\vec{F}_{ji}\!=\!-\vec{F}_{ij}\,\,\cdots\,\)(11)

(10)式を(11)式に代入すると

\(\quad\displaystyle\sum_i\vec{F}_i\cdot\vec{r}_i=\!\displaystyle\sum_{i,j(i\lt j)}\dfrac{a_{ij}}{|\vec{r}_i-\vec{r}_j|}\)
\(\qquad\qquad\qquad=\displaystyle\sum_{i,j(i\lt j)}\!V_{ij}=\!\displaystyle\!\sum_{i,j(i\lt j)}\!V_{ji}\quad\cdots\,\)(12)

よって(6)式は以下のように表すことが出来る。

\(\quad 0=2\langle\,T\,\rangle+\left\langle\displaystyle\sum_i\vec{F}_i\cdot\vec{r}_i\right\rangle=2\langle\,T\,\rangle+\left\langle\displaystyle\sum_{i,j(i\lt j)}V_{ij}\right\rangle\)
\(\qquad\qquad=2\langle\,T\,\rangle+\langle\,V\,\rangle\quad\cdots\,\)(13)

tamami.tech というドメインを取得したけど、更新料が 3,400 円と高額なので、tamami8.net という 更新料 1,480 円のドメインに変更しました。
まあ変更と行っても、”*.tech” のドメインは 2021年の6月まで有効ですが、HTML のソースを順次書き換え中です。

最近、ロリポップではサーバーレンタルとドメイン取得をまとめてやると、更新料無料らしい。まあ、しょうがないと諦めています。

ひもの両端を持って垂らしたときにできる曲線を懸垂曲線( カテナリー catenary )といいますが、この曲線の方程式を導いてみる。

ひもの底を原点\(O\,\)とし、水平方向に\(x\,\)軸、鉛直方向に\(y\,\)軸をとる。ひもの曲線を\(\,y=f(x)\,\)とおいて関数\(\,f(x)\,\)を求める。

ひもの\(\,x>0\,\)の部分に点\(\,(x,f(x))\,\)をとり、これを\(P\,\)とします。また、ひもの\(O\,\)から\(P\,\)の部分(両端を含む)を\(C\,\)とします。

【\(\,f(x)\,\)に対する微分方程式を立てる】

ひもの一部\(\,C\,\)に加わる力の釣り合いの式から\(\,f(x)\,\)についての微分方程式を立てる。\(C\,\)には次の３つの力が加わっている。

• 点\(\,O\,\)に水平方向のひもの張力\(\,T_0\)
• 点\(\,P\,\)に接線方向上向きにひもの張力\(\,T\)
• \(\,C\,\)の重心に鉛直下向きに重力\(\,G\)

ここで、原点にはたらく張力\(\,T_0\,\)は\(\,x\,\)に関係なく一定なので定数とする。

\(T\,\)について、\(P\,\)における\(\,f(x)\,\)の接線と\(\,x\,\)軸の正の部分とのなす角を\(\,\theta\,\)とすると、\(T\,\)の\(\,x\,\)成分、\(y\,\)成分はそれぞれ\(\,T\cos\theta\,,\,T\sin\theta\,\)となる。

この角度\(\,\theta\,\)と\(\,f(x)\,\)は、\(x\,\)における微分が接線なので
\(\quad f'(x)=\tan\theta\quad\)となる。

\(C\,\)の長さを\(\,l\,\)とすると曲線の長さの式から
\(\quad l=\displaystyle\int_0^x\sqrt{1+\{f'(x)\}^2}dx\quad\)で求められる。

線密度を\(\,\rho\,\)、重力加速度を\(\,g\,\)とすると
\(\quad G=\rho lg=\rho g\quad\)となる。

これらより釣り合いの式は以下のようになる。

• 鉛直方向：\(\,T\sin\theta=\rho lg\)
• 水平方向：\(\,T\cos\theta=T_0\)

\(T\,\)を消去して、定数をまとめて \(\,k=\dfrac{\rho g}{T_0}\,\,\)とおくと

\(\quad f'(x)=k\displaystyle\int_0^x\sqrt{1+\{f'(x)\}^2}dx\quad\)となる。

両辺を\(\,x\,\)で微分すると

\(\quad f^{\prime\prime}(x)=k\sqrt{1+\{f'(x)\}^2}\quad\)となる。

境界条件は\(\quad f(0)=0\,\,,\quad f'(0)=0\quad\)となる。

【微分方程式を解く】

この微分方程式を解く。\(\quad z=f'(x)\quad\)とおくと

\(\quad\dfrac{dz}{dx}=k\sqrt{1+z^2}\)
\(\quad\dfrac{dz}{\sqrt{1+z^2}}=k\,dx\quad\)となる。

先程の条件で\(\quad f'(x)=\tan\theta=z\quad\)なので
\(\quad dz=\dfrac{1}{\cos^2\theta}d\theta\quad\)となる。よって、

\(\quad \dfrac{dz}{\sqrt{1+z^2}}=\dfrac{1}{\sqrt{1+\tan^2\theta}}\dfrac{d\theta}{\cos^2\theta}=\dfrac{d\theta}{\sqrt{\dfrac{\cos^2\theta+\sin^2\theta}{\cos^2\theta}}\cos^2\theta}\)
\(\qquad =\dfrac{d\theta}{\cos\theta}\quad\)となる。

このままでは積分出来ないので、\(\quad\dfrac{d}{d\theta}\sin\theta=\cos\theta\quad\)を利用して

\(\quad\dfrac{d\theta}{\cos\theta}=\dfrac{d(\sin\theta)}{\cos^2\theta}=\dfrac{d(\sin\theta)}{(1-\sin\theta)(1+\sin\theta)}\)

\(\quad =\dfrac{1}{2}\left(\dfrac{1}{1-\sin\theta}+\dfrac{1}{1+\sin\theta}\right)\,d(\sin\theta)\quad\)となる。

この形なら積分可能なので\(\,\,\sin\theta\,\,\)で積分する。

\(\quad\displaystyle\int\dfrac{dz}{\sqrt{1+z^2}}=\dfrac{1}{2}\displaystyle\int\left(\dfrac{1}{1-\sin\theta}+\dfrac{1}{1+\sin\theta}\right)\,d(\sin\theta)\)

\(\quad=\dfrac{1}{2}\left\{-\log(1-\sin\theta)+\log(1+\sin\theta)\right\}\)

\(\quad=\dfrac{1}{2}\log\left(\dfrac{1+\sin\theta}{1-\sin\theta}\right)\)

よって\(\quad z=\tan\theta\quad\)として

\(\quad f'(x)=\displaystyle\int\dfrac{dz}{\sqrt{1+z^2}}=\dfrac{1}{2}\log\left(\dfrac{1+\sin\theta}{1-\sin\theta}\right)=kx+C_1\quad\) となる。

境界条件\(\quad f'(x)=0\quad x=0\,\,,\,\theta=0\quad\)より\(\quad C_1=0\,\,\)。

\(\quad\dfrac{1+\sin\theta}{1-\sin\theta}=e^{2kx}\quad\)より、まず\(\quad\sin\theta\quad\)について解く。

\(\quad 1+\sin\theta=e^{2kx}(1-\sin\theta)\qquad (e^{2kx}+1)\sin\theta=e^{2kx}-1\)

\(\quad \sin\theta=\dfrac{e^{2kx}-1}{2^{2kx}+1}\quad\)となる。ここから以下の公式を使う。

\(\quad 1+\dfrac{1}{\tan^2\theta}=\dfrac{\tan^2\theta+1}{\tan^2\theta}=\dfrac{\sin^2\theta+\cos^2\theta}{\sin^2\theta}=\dfrac{1}{\sin^2\theta}\)

\(\quad\dfrac{1}{\tan^2\theta}=\left(\dfrac{e^{2kx}+1}{2^{2kx}-1}\right)^2-1\)

\(\qquad=\dfrac{(e^{2kx}+1)^2-(e^{2kx}-1)^2}{(e^{2kx}-1)^2}\)

\(\qquad=\dfrac{4e^{2kx}}{(e^{2kx}-1)^2}\quad\)より

\(\quad\tan\theta=\dfrac{e^{2kx}-1}{2e^{kx}}=\dfrac{e^{kx}-e^{-kx}}{2}=\sinh kx\)

\(\quad f'(x)=z=\tan\theta=\sinh kx\quad\)なので

\(\quad f(x)=\displaystyle\int\dfrac{e^{kx}-e^{-kx}}{2}dx=\dfrac{e^{kx}+e^{-kx}}{2k}+C_2\)

境界条件が\(\quad f(0)=0\quad\)なので\(\quad C_2=-k\quad\)より

\(\quad f(x)=\dfrac{e^{kx}+e^{-kx}}{2k}-k=\dfrac{\cosh kx-1}{k}\)

1. ミクロの世界の量子状態

1-1. 箱の中の粒子

一次元の粒子の波動方程式を考える。\(\,x\,\)軸上に、\(\,x=0\,\)から\(\,x=L\,\)までの長さが\(\,L\,\)の区間(\(\,0\ge\,x\ge\,L\,\))を考え、この区間の内部でだけ粒子が運動できるモノとする。このような状況の粒子を、1次元の箱の中の粒子と呼ぶ。
箱の外に出られないように、箱の外の位置エネルギーを\(\,\infty\,\)とする。また、箱の中では粒子に力が働かず自由に運動できるように、位置エネルギーを\(\,0\,\)とする。
解くべき方程式は、定常状態の波動方程式 \(\,\hat{H}\phi(x)=E\phi(x)\,\) であり、粒子の質量を\(\,m\,\)とするとハミルトニアン\(\,\hat{H}\,\)は以下のようになる。

\(\quad \hat{H}=-\dfrac{\hbar^2}{2m}\Delta+U(x)\quad\cdots\,\)(1.1)

変数は\(\,x\,\)のみの1次式なので、ラプラシアン\(\,\Delta\,\)は\(\,x\,\)の2次微分 \(\,d^2/dx^2\,\)に等しい。箱の外では位置エネルギーが無限大だから、\(\,\hat{H}\phi=E\phi\,\)が成り立つためには、\(\phi=0\,\)にならなければならない。
すなわち、箱の外の確率密度は\(\,0\,\)になり、粒子は箱の外には存在しない。次に箱の中では、位置エネルギーを\(\,U(x)=0\,\)としたので、波動方程式は次のような2次微分を含む方程式になる。

\(\quad -\dfrac{\hbar^2}{2m}\ \dfrac{d^2\phi}{dx^2}=E\phi\quad\cdots\,\)(1.2)

この式では\(\,\phi\,\)と\(\,x\,\)以外はすべて定数であり、次の形の基本的な微分方程式となっている

\(\quad \dfrac{d^2\phi}{dx^2}=-k^2\phi\quad\cdots\,\)(1.3)

ここで右辺の\(\,k\,\)は、求めるべきエネルギー\(\,E\,\)を含み、次式で与えられる定数である

\(\quad k=\left(\dfrac{2mE}{\hbar^2}\right)^{1/2}\quad\cdots\,\)(1.4)

式 (1.3)の一般解は指数関数または三角関数を含む次式の形に表される

\(\quad \phi(x)=ae^{ikx}+be^{-ikx}=A\sin kx+B\cos kx\quad\cdots\,\)(1.5)

波動関数の重要な性質として、連続性の条件があり、境界条件を考慮しなければならない。
箱の外では\(\,\phi=0\,\)であるから、箱の中(\(\,0\ge x\ge L\,)\)の\(\,\phi\,\)も箱の両端で値が\(\,0\,\)にならなくてはいけない。
そこで境界条件として \(\,\phi(0)=\phi(L)=0\,\) を使うと、\(\,x=0\,\)で\(\,\cos\,\)関数の部分は\(\,0\,\)にならないので、\(\sin\,\)の関数部分だけが残る。
また、\(\,x=L\,\)で\(\,\sin kx\,\)の値が\(\,0\,\)になるには \(\,kL=n\pi\,\)となり、その結果、波動関数とエネルギーがそれぞれ次のように決まる。

\(\quad\phi(x)=A\sin(n\pi x/L)\quad\cdots\,\)(1.6)

\(\quad E=\dfrac{(n\pi\hbar/L)^2}{2m}\quad\cdots\,\)(1.7)

ここで整数\(\,n\,\)は\(\,1,2,3\,\)などの自然数であり、これは量子数と呼ばれる。

(1.7) 式から、箱の中の粒子に許されるエネルギーは、図のように飛び飛びのエネルギー準位になる。

この結果は有限な空間に閉じ込められた粒子のエネルギーが飛び飛びになって量子化されることを示している。
ここで注目するべきことは、エネルギーが最低の基底状態、すなわち\(\,n=1\,\)の状態でも、エネルギーが\(\,0\,\)にならないことである。このエネルギーをゼロ点エネルギーという。

(1.7) 式で与えられる波動関数は図の中央に示したように、エネルギー準位を区別する量子数\(\,n\,\)が増えるにつれ上下の振動が激しくなる。波動関数の値が\(\,0\,\)になる位置を節(ふし)という。

節になる位置は、\(\phi\,\)でも\(\,\phi^2\,\)も同じであり、図の右に示した\(\,\phi^2\,\)の図から明らかなように、節の位置に粒子が観測される確率は\(\,0\,\)である。

1-2. 調和振動子

波動方程式の簡単な例として、左図のようなバネ(バネ定数を\(\,k\,\))につるされた質量\(\,m\,\)の錘(おもり)の振動を考える。位置エネルギーは、\(\,U(x)=\frac{1}{2}kx^2\,\) であり、釣り合った位置から\(\,x\,\)だけ錘がずれると、図の右に示した放物線に沿って上がって行き、常に中心に引き戻される力が働いて振動する。このような振動を調和振動子という。

この場合のハミルトニアンの運動エネルギーは箱の中の粒子の場合と同じになり、位置エネルギーは\(\,U(x)=\frac{1}{2}kx^2\,\) であるから、この調和振動子の波動方程式は次のように表される。

\(\left(-\dfrac{\hbar^2}{2m}\dfrac{d^2}{dx^2}+\dfrac{1}{2}kx^2\right)\Psi=E\Psi\quad\cdots\,\)(1.8)

この微分方程式は、次のように波動関数を仮定すると解くことができる。

\(\quad\Psi(x)=e^{-ax^2}\,\displaystyle\sum_{i=0}b_ix^i\quad\cdots\,\)(1.9)

解いた結果、次式に従うエネルギー準位が得られる。

\(\quad E_n=\left(n+\dfrac{1}{2}\right)h\nu\qquad(n=0,1,2,\cdots)\quad\cdots\,\)(1.10)

ここで、\(\,n=0,1,2,\cdots\,\) は振動の量子数(振動量子数)であり、\(h\,\)はプランク定数で、\(\nu\,\)はこの振動の固有振動数である。

調和振動子に許されるエネルギー準位は図の左に示すように等間隔になる。その間隔は、振動のエネルギー量子\(\,h\nu\,\)であり,固有振動数\(\,\nu\,\)はバネの力の定数が\(\,k\,\)で質量が\(\,m\,\)の古典力学の振動子と全く同じで、次式によって与えられる。

\(\quad\nu=\dfrac{1}{2\pi}\sqrt{\dfrac{k}{m}}\quad\cdots\,\)(1.11)

量子論の振動子には、古典力学と明らかに違った注目すべき特徴がある。エネルギーが最低の状態(基底状態)でも、そのエネルギーは\(\,0\,\)にならず、\(\,\dfrac{1}{2}h\nu\,\)というエネルギーをもって振動する。このエネルギーをゼロ点エネルギーといい、その時の振動をゼロ点振動という。

図の左に示したように、調和振動子の波動関数は中心の周りで振動子、遠くに行くと減衰する。値が\(\,0\,\)になる節の数は量子数\(\,n\,\)に等しく、エネルギーが高いほど多くなる。

1-3. 2粒子系の運動方程式

2個の粒子を含む場合すなわち2粒子系の波動方程式は、以下のように取り扱うと1個の粒子の問題になり簡単になる。図の左に示したように、質量\(\,m_1\,\)と\(\,m_2\,\)の2個の球が距離$r$で結ばれているとして、その運動を考える。
ここで重心を座標原点に固定すると、\(\,r_1\,\)と\(\,r_2\,\)の和が\(\,r\,\)になり、2個の粒子の相対運動は長さが\(\,r\,\)の棒の先端につけた1つの球の運動の問題に単純化される。そのとき先端の粒子の質量\(\,\mu\,\)は、2個の粒子の質量の逆数の和の逆数になる。

\(\quad\)換算質量\(\quad\mu=\dfrac{1}{\frac{1}{m_1}+\frac{1}{m_2}}\)

こうして2個の粒子の相対運動は、換算質量をもつ1つの粒子の運動として取り扱うことが出来る。したがって、その波動方程式は次のように、1つの粒子の方程式と同じになる。

\(\quad\left(-\dfrac{\hbar^2}{2\mu}\Delta+U\right)\Psi=E\Psi\quad\cdots\,\)(1.13)

このような問題では、原点からの距離\(\,r\,\)が特別に重要なので、直交座標系\(\,(x,y,z)\,\)の代わりに、図の右に示す3次元の極座標\(\,(r,\theta,\varphi\,)\,\)を変数にとる。すると体積要素\(\,dxdydz\,\)とラプラシアン\(\,\Delta\,\)は、以下のように表される。

\(\quad\)体積要素\(\quad dxdydz=r^2\sin\theta dr d\theta d\varphi\quad\cdots\,\)(1.14)

\(\quad\)ラプラシアン\(\quad \Delta=\dfrac{1}{r^2}\dfrac{\partial}{\partial r}\left(r^2\dfrac{\partial}{\partial r}\right)+\dfrac{1}{r^2}\Lambda\quad\cdots\,\)(1.15)

\(\quad\)ルジャンドリアン\(\quad \Lambda=\dfrac{1}{\sin\theta}\dfrac{\partial}{\partial\theta}\left(\sin\theta\dfrac{\partial}{\partial\theta}\right)+\dfrac{1}{\sin^2\theta}\dfrac{\partial^2}{\partial\varphi^2}\quad\)(1.16)

ここで、(1.15)式のラプラシアンに含まれる\(\,\Lambda\,\)は角度に依存する演算子であり、ルジャンドリアンと呼ばれる。

1-4. 剛体回転子

2量子系の波動方程式は、換算質量\(\,\mu\,\)を使うと原点から距離\(\,r\,\)のところに繋がれた1個の粒子の回転運動になる。距離\(\,r\,\)が一定なので\(\,r\,\)に関する微分は消えてしまい、また位置エネルギーも\(\,0\,\)とおいてよいので、その波動方程式は次のように単純な式になる。

\(\quad-\dfrac{\hbar^2}{2I}\Delta\Psi=E\Psi\quad\cdots\,\)(1.17)

ここで、次式で表される\(\,I\,\)は、重心を軸とする回転運動の慣性モーメントである。

\(\quad I=\mu r^2\quad\cdots\,\)(1.18)

(1.17)式の波動方程式は三角関数を組み合わせた式を仮定して解くことが出来、その結果、剛体回転子に許されるエネルギーは次式のように、整数値をとる回転量子数\(\,J\,\)によって飛び飛びの準位になる。

\(\quad E=\dfrac{J(J+1)\hbar^2}{2I}\qquad (J=0,1,2,\cdots)\quad\cdots\,\)(1.19)

剛体回転子の波動方程式は、数学でよく知られている球面調和関数\(\,Y_{J,m}(\theta,\varphi)\,\)というものになり、それは\(\,J\,\)のほかにもう一つの量子数\(\,m\,\)にも依存している。\(\,J\,\)が決まれば、それに従って\(\,m\,\)のとりうる範囲が決まる。\(\,m\,\)は\(\,J\,\)から\(\,J-1,J-2,\cdots,0,\cdots,-J\,\)まで、合計\(\,2J+1\,\)通りの値が可能である。

\(\quad Y_{00}(\theta,\varphi)=\left(\dfrac{1}{4\pi}\right)^{1/2}\quad\cdots\,\)(1.20)
\(\quad Y_{10}(\theta,\varphi)=\left(\dfrac{3}{4\pi}\right)^{1/2}\cos\theta\quad\cdots\,\)(1.21)
\(\quad Y_{1\pm1}(\theta,\varphi)=\mp\left(\dfrac{3}{8\pi}\right)^{1/2}\sin\theta\ e^{\pm i\varphi}\quad\cdots\,\)(1.22)
\(\quad Y_{20}(\theta,\varphi)=\left(\dfrac{5}{16\pi}\right)^{1/2}(3\cos^2\theta-1)(\quad\cdots\,\)(1.23)
\(\quad Y_{2\pm1}(\theta,\varphi)=\mp\left(\dfrac{15}{8\pi}\right)^{1/2}\sin\theta\ \cos\theta\ e^{\pm i\varphi}\quad\cdots\,\)(1.24)
\(\quad Y_{2\pm2}(\theta,\varphi)=\mp\left(\dfrac{15}{32\pi}\right)^{1/2}\sin^2\theta\ e^{\pm i2\varphi}\quad\)(1.25)

2.分子の回転

2-1. 回転運動のエネルギー

質量\(\,m\,\)、速度\(\,v\,\)の直線運動の大きさが運動量\(\,p=mv\,\)で、運動に伴うエネルギーは

\(\quad E=\dfrac{1}{2}mv^2=\dfrac{1}{2}\dfrac{p^2}{m}\quad\cdots\,\)(2.1)

である。一方、回転運動の大きさを表す量が角運動量であり、回転の軸から距離\(\,r\,\)にあって、運動量\(\,p\,\)で運動している物体の角運動量\(\,p_{\theta}\,\)は、

\(\quad p_{\theta}=r\times p\quad\cdots\,\)(2.2)

と表される。\(\,p_{\theta}\,\)は運動量を含むベクトルである。この角運動量を用いて、回転運動のエネルギーは

\(\quad E=\dfrac{1}{2}\dfrac{p_{\theta}^2}{I}\quad\cdots\,\)(2.3)

のように表される。式(2.1)と(2.3)を比較すると対応していることが分かる。

直線運動を行う物体で運動の持続性(慣性)に影響する量である質量に対応するのが、回転運動における\(\,I\,\)であり、慣性モーメントと呼ばれる。慣性モーメントは

\(\quad I=\sum_im_ir_i^2=\sum_im_i(x_i^2+y_i^2+z_i^2)\quad\)(2.4)

であるが、2原子分子では換算質量\(\,\mu\,\)と原子間の距離\(\,r\,\)を使うと、以下の簡単な式で表すことが出来る

\(\quad \mu=\dfrac{m_1\cdot m_2}{m_1+m_2}\quad\cdots\,\)(2.5)

\(\quad I=\mu r^2\quad\cdots\,\)(2.6)

式(2.3)には運動量\(\,p_{\theta}\,\)が含まれており、この式を量子の世界にもちこむには、その\(\,p_{\theta}\,\)を演算子に置き換えればよい。その結果

\(\quad p_{\theta}=\sqrt{J(J+1)}\hbar\quad\cdots\,\)(2.7)

と置き換えることが出来る。従って、式(2.3)の2原子分子の場合の量子論的な表現は式(1.19)のように

\(\quad E=\dfrac{J(J+1)\hbar^2}{2I}=BJ(J+1)\qquad (J=0,1,2,\cdots)\quad\cdots\,\)(2.8)

となる。ここで\(\,J\,\)は回転量子数であり、\(\,0\,\)を含む整数のみが許される。すなわち回転エネルギーは量子化されており、飛び飛びの準位のみが許されることになる。\(\,B\,\)は回転定数と呼ばれるもので次式で表される。

\(\quad B=\dfrac{h}{8\pi^2I}\quad\)周波数単位}\(\qquad \dfrac{h}{8\pi^2c_0I}\quad\)波数単位\(\quad\cdots\,\)(2.9)

2-2. 観測される回転スペクトル

回転運動だけに起因するスペクトルはマイクロ波領域から遠赤外線領域に現れる。図は期待状態のCO分子の吸収スペクトルである。規則的な間隔の吸収線からなる。式(2.8)より2原子分子の場合の回転エネルギー準位は下図のようになることが理論的に予測できる。

回転量子数の\(\,J=0,1,2,3,\cdots\,\)に対して、エネルギー準位は\(\, 0,2B,6B,12B,\cdots\,\)となる。従って隣接するエネルギー準位間の大きさは \(\,2B,4B,6B,\cdots\,\)である。これを横軸にエネルギーをとってプロットすると図の下に示したようにスペクトルの位置が\(\,2B\,\)の間隔で並ぶ。これは実測のスペクトルとぴったり一致する。

このようにして、光子を1つ吸収することによる(2原子分子の)回転遷移は、\(\,J\,\)の値が1つだけ変化する状態の間に起こることがわかる。吸収はエネルギーが増える方向への遷移であるが、エネルギーが減る方向の瀬には発光スペクトルとして観測され、実際同様なスペクトルを示すことがわかっている。すなわち回転遷移の選択則は\(\,\Delta J=\pm 1\,\)である。ここでは、この選択則を実測との対応から導いたが、理論的にも同じ結果が導かれる。

以上のようにして、実測スペクトルの間隔から回転定数\(\,B\,\)が決定できる。この\(\,B\,\)は(2.9)式に定義された量であるから、\(\,B\,\)がわかるということは慣性モーメント\(\,I\,\)が実験的に求められたことを意味する。これより\(\,B\to I\to r\,\)のようにして、分子の構造を与える情報\(\,r\,\)が得られることになる。

2原子分子であるNO分子についても同様な解析が可能であり、実測スペクトルから分子構造に関する情報を導くことができる。それらの結果を以下の表にまとめる。

COの場合、\(J\!=\!1\!-\!0\,\)の回転遷移では、\(E=BJ(J+1)\,\) より、\(115.9\)GHz となり、\(^{13}\)COでは、\(110.8\)GHzとなる。\(J\!=\!2\!-\!1\,\)の場合は \(\,6B-2B=4B\,\)なので、\(^{12}\)CO は\(231.7\)GHzで、\(^{13}\)COは\(\,221.5\,\)GHzになる。

2-3. 観測されたスペクトル

爆発的星形成銀河（starburst銀河）M82の中心数100pc領域における、赤外線から電波に至る波長領域のスペクトル。大別して、(1)周波数に対して、連続的に放射強度が変化している「連続波放射」と、(2)ある特定の周波数でのみ、強い放射を示す「スペクトル線放射」とに分けられる。
連続波放射は、主に波長の長い電波域におけるシンクロトロン放射(磁場の中を高速で運動している電子からの放射)と、サブミリ波～赤外線域にかけて、大きな山をつくっている成分、すなわち、星間ダスト(星間空間中に存在する、シリケイトなどの固体微粒子)からの熱放射の重ね合わせとしてよく説明でできる。
この他、ミリ波付近には自由・自由遷移放射(電離領域中を運動する自由電子からの放射)も見られる。一方、この波長帯のスペクトル線放射としては、原子からの再結合線(桃色)・微細構造線(赤)・超微細構造線、分子からの回転線(青)などがみられる。

Google Chrome を使っていると、http でのアクセスの場合、保護されていないページと表示されるので、SSL化をしようとトライしてみました。
ドメインはデフォルトのままでは、無料SSL導入が出来ない(ロリポップでは)らしいので、しょうがないからドメインを導入しました。”tamami“という名前は、1970年に桜花賞をぶっちぎりで勝った「タマミ」からとっているので、これをそのまま使うことにしました。*.tokyo だと安かったんですけど、何となく “.tech“にしてしまいました。よく値段を確認しなかったんですけど、更新費用が年額 3,980円でサーバーのレンタル料(3,000円)より高い．．．．
結局、GMO としてこういう形でコストを回収しているんだっと言うことを、やっと勉強した。😭

これで、ようやくSSL化が出来ました。内部のアドレス変換で RegEx というツールで変換したのが、失敗で \(\TeX\,\)の制御コードのバックラッシュ “\”が消えてしまいました。何か設定があるんだと思うけど、やってしまったことはしょうがないので、手動で復旧した。これで、 Google Chrome から文句は言われないかな。

【問４】

\(\quad x^5=i\quad\)を解け。

【解法】
\(\quad i^5=i\quad\)なので、与式は\(\,\,x^5-i^5=0\quad\)と変形できる。

\(\quad x^5-i^5=(x-i)(x^4+x^3i+x^2i^2+xi^3+i^4)\)
\(\qquad =(x-i)(x^4+x^3i-x^2-xi+1)\quad\)ここで、\(x\ne0\quad\)なら

\(\qquad =(x-i)\left(x^2+xi-1-\dfrac{1}{x}+\dfrac{1}{x^2}\right)=0\)

ここで、\(\,x-\dfrac{1}{x}=t\,\,\)とすると、
\(\quad t^2+2+it-1=t^2+it+1=0\quad\)より

\(\quad\left(t+\dfrac{1}{2}i\right)^2+\dfrac{5}{4}=0\quad\)なので

\(\quad t+\dfrac{1}{2}i=\pm\dfrac{\sqrt{5}}{2}i\quad\)よって\(\quad t=\dfrac{-1\pm\sqrt{5}i}{2}\)

ここで\(\,\,x=\cos\theta+i\sin\theta\quad\)とおくと

\(\quad t=\cos\theta+i\sin\theta-(\cos\theta+i\sin\theta)^{-1}\)

\(\qquad =\cos\theta+i\sin\theta-\cos\theta+i\sin\theta=(2\sin\theta)i\)

\(\quad \sin\theta=\dfrac{-1\pm\sqrt{5}}{4}\quad\)となる

\(\quad \cos^2\theta=1-\sin^2\theta=\dfrac{10\pm2\sqrt{5}}{16}\quad\)より

\(\quad \cos\theta=\pm\dfrac{\sqrt{10\pm2\sqrt{5}}}{4}\quad\)なので、解は \(\,i\,\)を頂点とする５角形の点。

\(\quad x=i\,\,\,,\,\,\pm\dfrac{\sqrt{10+2\sqrt{5}}}{4}+\dfrac{-1+\sqrt{5}}{4}i\,\,,\)

\(\qquad\qquad\quad \pm\dfrac{\sqrt{10-2\sqrt{5}}}{4}+\dfrac{-1-\sqrt{5}}{4}i\)

また、\(\,i=\cos 5\theta+i\sin 5\theta=\cos\dfrac{\pi}{2}+i\sin\dfrac{\pi}{2}\quad\)なので

\(\quad \theta=\dfrac{\pi}{10}\quad\)となり、\(\sin\dfrac{\pi}{10}=\dfrac{-1+\sqrt{5}}{4}\quad\)である。

【問５】

\(\quad \displaystyle\int_0^{\pi}e^x\sin^2x\,dx>8\quad\)を示せ

【解答】
\(\quad I=\displaystyle\int_0^{\pi}e^x\sin^2x\,dx\quad\)とする。

\(\quad\sin^2x=\dfrac{1-\cos2x}{2}\quad\)なので、\(\,I=\dfrac{1}{2}\displaystyle\int_0^{\pi}e^x(1-\cos2x)dx\)

\(\quad I=\dfrac{1}{2}\left\{\displaystyle\int_0^{\pi}e^x\,dx-\displaystyle\int_0^{\pi}e^x\cos2x\,dx\right\}\)

ここで、\(\displaystyle\int_0^{\pi}e^x\cos2x\,dx=I_2\quad\)とおく。さきに微分を計算する。

\(\quad e^x\sin2x=a\,\,,\quad e^x\cos2x=b\quad\)とおく

\(\quad(e^x\cos2x)^{\prime}=e^x\cos2x-2e^x\sin2x=b-2a\quad\cdots\,\)(1)

\(\quad(e^x\sin2x)^{\prime}=e^x\sin2x+2e^x\cos2x=a+2b\quad\cdots\,\)(2)

\(\quad b=\dfrac{(1)\!+\!(2)\!\times\!2}{5}=\dfrac{1}{5}(e^x\cos2x+2e^x\sin2x)^{\prime}=e^x\cos2x\)

\(\quad\therefore\,\,I_2=\dfrac{1}{5}\bigl[e^x\cos2x+2e^x\sin2x\bigr]_0^{\pi}=\dfrac{1}{5}(e^{\pi}-1)\)

\(\quad I=\dfrac{1}{2}\left(e^{\pi}-1-\dfrac{1}{5}(e^{\pi}-1)\right)=\dfrac{2}{5}(e^{\pi}-1)\quad\)となるので

\(\quad e^{\pi}>21\quad\)を示せばいいことになる。

ここで\(\,\,y=e^x\quad\)のグラフで\(\quad x=3\quad\)での接線を考える。

接線の式は\(\quad y=e^3(x-3)+e^3=e^3x-2e^3\quad\)となる。

\(\quad e^{\pi}\quad\)は、接線の\(\,\,x=\pi\quad\)の点より上にあるので

\(\quad e^{\pi}>e^3\pi-2e^3=e^3(\pi-2)>2.7^3\times1.1\fallingdotseq21.65\)

よって、標題は証明された。

【問６】

\(\,x>0\,\)のとき\(\quad\dfrac{x^4+x^2+1}{x^3+x}\quad\)の最小値を求めよ

【解答】
与式の分母、分子を\(\quad x^2\quad\)で割ると

\(\quad =\dfrac{x^2+1+\frac{1}{x^2}}{x+\frac{1}{x}}\quad\cdots\,\)(3)　となる。

ここで\(\quad x+\dfrac{1}{x}=t\quad\cdots\,\)(4)　とおくと、式(3)は

\(\quad \dfrac{t^2-1}{t}=t-\dfrac{1}{t}\quad\cdots\,\)(5)　となる。

ここで\(\quad x>0\quad\dfrac{1}{x}>0\quad\)なので、

\(\quad t=x+\dfrac{1}{x}\ge 2\sqrt{x\dfrac{1}{x}}=0\quad\)なので、\(\,t\ge 2\)

ここで、式(5)を\(\quad f_1(t)=t\quad f_2(t)=-\dfrac{1}{t}\quad\)と分けると

\(\quad t\ge 2\quad\)で、\(f_1(t)\,\,,\,f_2(t)\quad\)はともに単調増加

なので、最小値は\(\quad t=2\,\to\,x=1\quad\)最小値は\(\quad\dfrac{3}{2}\)

最近、Youtube で流れている問題を見ていると、すっかり忘れている事や知らないことが多い。そこで、いくつかの問題とその解法を備考録に書いていきます。

【問１】

\(\quad\sqrt{3}\sin 2x+2\sin^2x-1\,\) の\(\,0\!\le\!x\!\le\!\pi\,\)での最大値と最小値を求めよ。

【解法】
三角関数って、やっていないと直ぐ忘れる。この場合は倍角の公式を使う。
\(\quad \cos 2\alpha=1-2\sin^2\alpha\qquad\cdots\,\)(1)
これを使うと与式は
\(\quad\sqrt{3}\sin 2x-\cos 2x=2\left(\dfrac{\sqrt{3}}{2}\sin 2x+\left(-\dfrac{1}{2}\right)\cos 2x\right)\,\cdots\,\)(2)

ここで\(\,\sin\,\)の加法定理が
\(\quad\sin(\alpha+\beta)=\sin\alpha\cos\beta+\cos\alpha\sin\beta\qquad\cdots\,\)(3)

であり、\(\cos(-\pi/6)=\sqrt{3}/2\,,\quad\sin(-\pi/6)=-1/2\,\)なので、(2)式は

\(\quad 2\left(\sin 2x\cos\left(-\dfrac{\pi}{6}\right)+\cos 2x\sin\left(-\dfrac{\pi}{6}\right)\right)\qquad\cdots\,\)(4)

(4)と加法定理の(3)を見比べれば
与式は以下のように変形できる。

\(\quad 2\left(\sin\left(2x-\dfrac{\pi}{6}\right)\right)\qquad\cdots\,\)(5)

最大値は\(\,2x-\pi/6=\pi/2\,\)のときなので \(\,\,x=\pi/3\,\)のとき\(\,\,2\)
最小値は\(\,2x-\pi/6=3\pi/2\,\)のときなので \(\,\,x=5\pi/6\,\)のとき\(\,\,-2\)

【問２】

\(\quad x^5+x+1\quad\)を因数分解せよ。

【解法】
ネットでは、まず１の３乗根を使うとあります。これは以下の方程式の

\(\quad x^3-1=(x-1)(x^2+x+1)=0\qquad\cdots\,\)(6)

根は、\(x=1,\,(-1\pm\sqrt{3}\,i)/2\,\)である。

\(\omega=(-1+\sqrt{3}\,i)/2\,\)とすると

\(\quad\omega^2=(-1-\sqrt{3}\,i)/2\,,\quad\omega^3=1\,\)である。

\(\quad\omega^2+\omega+1=0\,\)である。与式に\(\,\omega\,\)を代入すると

\(\quad\omega^5+\omega+1=\omega^2+\omega+1=0\,\)なので割り切れる。

そこで、与式を\(\,x^2+x+1\,\)で割ってみる。

このように、割り切れる。商が\(\,x^3-x^2+1\,\)となるので、答えは

\(\quad x^5+x+1=(x^2+x+1)(x^3-x^2+1)\quad\)となります。

別解(?)として、以下のように展開する方式があります。

\(\quad x^5+x+1=x^5-x^2+x^2+x+1=x^2(x^3-1)+x^2+x+1\)
\(\qquad =x^2(x-1)(x^2+x+1)+x^2+x+1\)
\(\qquad =(x^2+x+1)(x^3-x^2+1)\)

【問３】

\(\quad a^2+b^2=1224\quad\)にあてはまる自然数の組を答えよ。

【解法】
この手の整数の問題は合同式を使うのがセオリーらしい。

\(\,{\large \bigstar\,}\)合同式
一般に、ある２つの整数\(\,a,\,b\,\)を自然数\(\,m\,\)で割った余りが等しいとき、\(a,\,b\,\)は\(\,m\,\)を法として合同であるといい
\(\quad a\equiv b\,(\mathrm{mod}\, m)\qquad\)と表す。

この式のことを合同式といい、\(a-b\,\)が\(\,m\,\)で割りれることを表す。
合同式では加法・減法・乗法が成り立つ。

\(a\equiv b\,(\mathrm{mod}\,m),\,\,c\equiv d\,(\mathrm{mod}\,m)\,\)のとき
\(\quad a+c\equiv b+d\,(\mathrm{mod}\,m)\)
\(\quad a-c\equiv b-d\,(\mathrm{mod}\,m)\)
\(\quad ac\equiv bd\,(\mathrm{mod}\,m)\)

除法については\(\, ac\equiv bc\,(\mathrm{mod}\,m)\,\)で
\(\quad\cdot\,c\,\)と\(\,m\,\)が互いに素のとき、
\(\qquad a\equiv b\,(\mathrm{mod}\,m)\)

\(\quad\cdot\,c\,\)と\(\,m\,\)が互いに素でないとき、\(c\,\)と\(\,m\,\)の最大公約数を\(\,d(>1)\,\)とすると
\(\qquad a\equiv b\,(\mathrm{mod}\,\dfrac{m}{d})\quad\)が成立する。

\(\,{\large \bigstar\,}\)平方剰余と平方数の性質

\(\quad a^2\equiv b\,(\mathrm{mod}\,m)\,\)となる整数\(\,a\,\)が存在するとき、\(b\,\)は\(\,m\,\)を法とする平方剰余であるという。

・平方数に関する重要な性質
\(1\,\)：平方数を\(\,3\,\)で割った余りは必ず\(\,0\,\)か\(\,1\,\)。余りが\(\,2\,\)になることはない。
\(2\,\)：平方数を\(\,4\,\)で割った余りは必ず\(\,0\,\)か\(\,1\,\)。余りが\(\,2\,\)か\(\,3\,\)になることはない。

これらをふまえて、解法は以下のようになる

\(\,1224=2^3\,3^2\,17\,\)なので
・\(\,a^2+b^2=1224=4\times 306\,\)より

\(\quad a\equiv0\,(\mathrm{mod}\,2)\,\)かつ\(\,b\equiv0\,(\mathrm{mod}\,2)\)

・\(\,a^2+b^2=1224=9\times 136\,\)より

\(\quad a\equiv0\,(\mathrm{mod}\,3)\,\)かつ\(\,b\equiv0\,(\mathrm{mod}\,3)\quad\)なので
\(\therefore\quad a\equiv0\,(\mathrm{mod}\,6)\,\)かつ\(\,b\equiv0\,(\mathrm{mod}\,6)\quad\)となるので
\(\quad a=6x,\,b=6y\,\)とおくと\(\,\,x^2+y^2=34\,\,\)になり、

\(\quad (x,y)=(5,3)\,,\,(3,5)\,\)となるので\(\,(a,b)=(30,18)\,or\,(18,30)\,\)

特殊相対論の質量とエネルギーの関係式を確認してみる

1．前準備

1-1. テイラー展開

閉区間 [\(\,a\,,\,b\,\)]で\(\,n\,\)回微分可能な関数\(\,f(x)\,\)において、\(n\to\infty\,\)にて

\(f(b)=f(a)\!+\!\dfrac{f'(a)}{1!}(b\!-\!a)\!+\!\dfrac{f”(a)}{2!}(b\!-\!a)^2\cdot\!+\!\dfrac{f^{(n)}(a)}{n!}(b\!-\!a)^n\,\)(1)

とできる。この (1)式を\(\,b=x\,\)とするのが別形である。
そして\(\,a=0\,\)とするのが　マクローリン展開( Maclaurin expansion )　である。

\(\quad f(x)=\displaystyle\sum_{n=0}^{\infty}\dfrac{f^{(n)}(0)}{n!}\,x^n\quad\cdots\,\)(2)

ここでは、ローレンツ変換で用いる関数\(\,\,f(x)=(1-x^2)^{-1/2}\,\,\)をマクローリン展開する。

1-2. \(\,\,f(x)\,\)の微分

\(\quad f(x)=\dfrac{1}{\sqrt{1-x^2}}=(1-x^2)^{-\frac{1}{2}}\quad\cdots\,\)(3)

を微分する。まず１回微分すると

\(\quad f'(x)=(-\dfrac{1}{2})(-2x)(1-x^2)^{-\frac{3}{2}}=x(1-x^2)^{-\frac{3}{2}}\quad\cdots\,\)(4)

２回微分では

\(\quad f^{\prime\prime}(x)=(1-x^2)^{-\frac{3}{2}}+(-\dfrac{3}{2})(-2x)x(1-x^2)^{-\frac{5}{2}}\)
\(\qquad=(1-x^2)^{\frac{3}{2}}+3x^2(1-x^2)^{-\frac{5}{2}}\quad\cdots\,\)(5)

さらに微分すると

\(\quad f^{\prime\prime\prime}(x)=(-\dfrac{3}{2})(-2x)(1-x^2)^{-\frac{5}{2}}+6x(1-x^2)^{-\frac{5}{2}}\)
\(\qquad +3x^2(-\dfrac{5}{2})(-2x)(1-x^2)^{-\frac{7}{2}}\)
\(\qquad =3x(1-x^2)^{-\frac{5}{2}}+6x(1-x^2)^{-\frac{5}{2}}+15x^3(1-x^2)^{-\frac{7}{2}}\)
\(\qquad =9x(1-x^2)^{-\frac{5}{2}}+15x^3(1-x^2)^{-\frac{7}{2}}\quad\cdots\,\)(6)

\(\quad f^{(4)}(x)=9(1-x^2)^{-\frac{5}{2}}+9x(-\dfrac{5}{2})(-2x)(1-x^2)^{-\frac{7}{2}}\)
\(\qquad+45x^2(1-x^2)^{-\frac{7}{2}}+15x^3(-\dfrac{7}{2})(-2x)(1-x^2)^{-\frac{9}{2}}\)
\(\qquad =9(1-x^2)^{-\frac{5}{2}}+90x^2(1-x^2)^{-\frac{7}{2}}+105x^4(1-x^2)^{-\frac{9}{2}}\quad\)(7)

\(\quad f^{(5)}(x)=225x(1-x^2)^{-\frac{7}{2}}+1050x^3(1-x^2)^{-\frac{9}{2}}\)
\(\qquad +945x^5(1-x^2)^{-\frac{11}{2}}\quad\cdots\,\)(8)

\(\quad f^{(6)}(x)=225(1-x^2)^{-\frac{7}{2}}+4725x^2(1-x^2)^{-\frac{9}{2}}\)
\(\qquad +7875x^4(1-x^2)^{-\frac{11}{2}}+10395x^6(1-x^2)^{-\frac{13}{2}}\quad\cdots\,\)(9)

1-3. \(\,\,f(x)\,\)のマクローリン展開

\(\quad f^{(n)}(x)\,\)に\(\,x=0\,\)を代入すると、

\(\quad f(0)=1\,,\quad f'(0)=0\,,\quad f^{\prime\prime}(0)=1\,,\quad f^{\prime\prime\prime}(0)=0\)
\(\quad f^{(4)}(0)=9\,,\quad f^{(5)}(0)=0\,,\quad f^{(6)}(0)=225\)

なので、それぞれの値を(2)式に代入すると

\(\quad f(x)=(1-x^2)^{-\frac{1}{2}}=1+\dfrac{1}{2}x^2+\dfrac{3}{8}x^4+\dfrac{5}{16}x^6+\cdots\,\,\)(10)

のように展開することが出来る

2．４次元時空での運動

2-1. 固有時間

慣性系\(\,S(t,x,y,z)\,\)と慣性系\(\,S'(t’,x’,y’,z’)\,\)において、
\(\quad dt=t’-t\,,\quad dx=x’-x\,,\quad dy=y’-y\,,\quad dz=z’-z\)

とすると、事象の隔たりを　世界距離( world distance :\(\,ds\,\) )といい、次の式で表される。

\(\quad ds^2=c^2dt^2-dx^2-dy^2-dz^2\quad\cdots\,\)(11)

ここで空間座標をまとめて\(\,\mathbf{x}\,\)で表すと

\(\quad ds^2=c^2dt^2-d\mathbf{x}^2\quad\cdots\,\)(12)

慣性系\(\,S’\,\)が\(\,S\,\)に対して、速度\(\,v\,\)で等速運動しているとすると、

ローレンツ変換 (Lorentz transformation )は

\(\quad\left(\begin{array}{c}\bf{x’}\\ct’\end{array}\right)=\gamma\left(\begin{array}{cc}1&\beta\\ \beta&1\end{array}\right)\left(\begin{array}{c}\bf{x}\\ct\end{array}\right)\quad\cdots\,\)(13)

ここで、\(\beta=\dfrac{v}{c}\,\)であり\(\quad \gamma=\dfrac{1}{\sqrt{1-\beta^2}}\,\)である。

この\(\,ds^2\,\)はローレンツ変換で不変である。

\(\quad ds^2=-c^2dt^2+d\mathbf{x}^2=-c^2dt’^2+d\mathbf{x}’^2\quad\cdots\,\)(14)

世界線の長さは\(\,\int\sqrt{-ds^2}\,\)として定義出来る。パラメータとして慣性系\(\,S\,\)の時間座標\(\,t\,\)そのものを選ぶことが出来る。

\(\quad\displaystyle\int\sqrt{-ds^2}=c\displaystyle\int dt\sqrt{1-\dfrac{1}{c^2}\left|\dfrac{d\mathbf{x}}{dt}\right|^2}=c\displaystyle\int d\tau\quad\cdots\,\)(15)

ここで時間定数\(\,\tau\,\)を定義する。固有時 (proper time )。

\(\quad d\tau=dt\sqrt{1-\dfrac{1}{c^2}\left|\dfrac{d\mathbf{x}}{dt}\right|^2}=dt’\sqrt{1-\dfrac{1}{c^2}\left|\dfrac{d\mathbf{x}’}{dt’}\right|^2}\quad\cdots\,\)(16)

式(15)は世界線に沿っての固有時間の長さを距離の単位で与える。これは、質点の軌道の形そのものによって定まる量なので、ローレンツ変換で不変である。

2-2. ４元速度

時間\(\,t\,\)を加えた４次元で固有時間に対する速度を考える。

\(\quad u_t=\dfrac{dct}{d\tau}\,,\,u_x=\dfrac{dx}{d\tau}\,,\,u_y=\dfrac{dy}{d\tau}\,,\,u_z=\dfrac{dz}{d\tau}\quad\cdots\,\)(17)

ここで、速度の大きさを考える場合、世界距離の計算と同様に時間成分の２乗はマイナスとなる。大きさを計算(内積をとる)する。

\(\quad-u_t^2+u_x^2+u_y^2+u_z^2=\dfrac{1}{d\tau^2}(-dct^2+dx^2+dy^2+dz^2)\)
\(\qquad =\dfrac{1}{d\tau^2}(-c^2d\tau^2)=-c^2\quad\cdots\,\)(18)

となり、４元速度の大きさは一定値になってしまう。
ここで、四元速度の表記を以下のようにする。

\(\quad u^{\mu}\equiv\dfrac{dx^{\mu}}{d\tau}\quad\cdots\,\)(19)

なお、\(u^0=u_t\,,\,u^1=u_x\,,\,u^2=u_y\,,\,u^3=u_z\,\)である。

式(16)、(17)より

\(\quad u^{\mu}=\dfrac{1}{\sqrt{1-\beta^2}}\dfrac{dx^{\mu}}{dt}=\dfrac{1}{\sqrt{1-\beta^2}}v^{\mu}\quad\cdots\,\)(20)

また、

\(\quad u^0=c\dfrac{dt}{d\tau}=\gamma\dfrac{cdt}{dt}=c\dfrac{1}{\sqrt{1-\beta^2}}\quad\cdots\,\)(21)

2-3. ４元運動量
４元速度を用いて４元運動量を\(\,m\,\)を静止質量として以下のように定義する。

\(\quad p^{\mu}=m\,u^{\mu}=m\gamma v^{\mu}\quad\cdots\,\)(22)

これにより、\(p^0\,\)は

\(\quad p^0=mu^0=mc\gamma=\dfrac{mc}{\sqrt{1-\beta^2}}\quad\cdots\,\)(23)

ここで、前準備の式(10)より

\(\quad\dfrac{1}{\sqrt{1-\beta^2}}=1+\dfrac{1}{2}\beta^2+\dfrac{3}{8}\beta^4+\dfrac{5}{16}\beta^6+\cdots\quad\cdots\,\)(24)

\(\,\beta=v/c\,\)なので、式(23)と(24)より

\(\quad p^0=mc+\dfrac{mc}{2}\dfrac{v^2}{c^2}+\dfrac{3mc}{8}\dfrac{v^4}{c^4}+\cdots\quad\cdots\,\)(25)

ここで、\(v/c\to 0\,\)の低速において、式(25)の両辺に光速\(\,c\,\)をかけると

\(\quad p^0c=mc^2+\dfrac{1}{2}mv^2+\cdots\qquad\cdots\,\)(26)

\(\,p^0c\,\)は質点のエネルギーを表しており、運動エネルギーと
静止エネルギー\(\,mc^2\,\)の和として表される。

1.プランクの法則 (Planck’s law)
物理学における黒体から輻射(放射)される電磁波の分光放射輝度、もしくはエネルギー密度の波長分布に関する公式。プランクの公式とも呼ばれる。ある温度\(\,T\,\)における黒体からの電磁輻射の分光放射輝度を全波長領域において正しく説明することができる。1900年、ドイツの物理学者マックス・プランクによって導かれた。

プランクはこの法則の導出を考える中で、輻射場の振動子のエネルギーが、あるエネルギー素量(現在ではエネルギー量子と呼ばれている)\(\,\epsilon=h\nu\,\)の整数倍になっていると仮定した。このエネルギーの量子仮説（量子化）はその後の量子力学の幕開けに大きな影響を与えている。
 プランクの放射式では、単位振動数\(\,\nu\,\)あたりの黒体輻射強度を以下のように表す。

\(\quad B_{\nu}(T)=\dfrac{2h}{c^2}\dfrac{\nu^3}{e^{h\nu/kT}-1}\quad\cdots\,\)(1)

なお、ここで\(\,c\,\)は光速、\(\,h\,\)はプランク定数、\(\,k\,\) はボルツマン定数である。このプランクの放射式に温度\(\,T\,\)に3000K\(\sim\)7000Kを代入してプロットしたのが下図である。

2.全黒体輻射強度 (Total blackbody radiation intensity)
このプランクの放射式を全振動数で積分したものを、全黒体輻射強度\(\,B(T)\,\) という。(1)式で、\(\,x=h\nu/kT\,\)とおくと、

\(\quad B_x(T)=\dfrac{2k^3T^3}{c^2h^2}\dfrac{x^3}{e^x-1}\quad\cdots\,\)(2)

となり、全振動数で積分すると

\(\quad \displaystyle\int_0^{\infty}B_{\nu}(T)d\nu=\displaystyle\int_0^{\infty}B_x(T)\,\dfrac{kT}{h}\,dx\quad\cdots\,\)(3)

\(\quad B(T)=\dfrac{2k^4T^4}{c^2h^3}\displaystyle\int_0^{\infty}\dfrac{x^3}{e^x-1}\,dx\quad\cdots\,\)(4)

(4)式の積分は\(\,\pi^4/15\,\)となるので、
全黒体輻射強度(放射輝度)[単位 Wm\(^{-2}\)sr\(^{-2}\)K\(^{-4}\)] は

\(\quad B(T)=\dfrac{2\pi^4k^4T^4}{15c^2h^3}=\dfrac{\sigma}{\pi}T^4\quad\cdots\,\)(5)
となる。

3.\(\,\,\,\,x^3/e^x-1\,\)の\(\,\,\)積分の計算

\(\quad\dfrac{x^3}{e^x-1}=\dfrac{x^3}{e^x}\dfrac{1}{1-e^{-x}}\quad\cdots\,\)(6)

(6)式の最後の項は公比\(\,e^{-x}\,\)の等比数列の和とすると

\(\quad\dfrac{1}{1-e^{-x}}=1+e^{-x}+e^{-2x}+\cdots\)
となるので

\(\quad\dfrac{x^3}{e^x-1}=x^3(e^{-x}+e^{-2x}+e^{-3x}+\cdots)\quad\cdots\,\)(7)

ここで\(\,x^pe^{-ax}\,\)の形の関数を積分するために、ガンマ関数を利用する。
ガンマ関数の定義から

\(\quad\varGamma(q)=\displaystyle\int_0^{\infty}y^{q-1}e^{-y}dy=(q-1)\,!\quad\cdots\,\)(8)

ここで、\(\,p=q-1\,\)とすると
\(\quad p\,!=\displaystyle\int_0^{\infty}y^pe^{-y}dy\quad\cdots\,\)(9)

\(\,\,y=ax\,\)とすると\(\,dy=a\,dx\,\)で、(9)式に代入すると

\(\quad p\,!=\displaystyle\int_0^{\infty}(ax)^pae^{-ax}dx=a^{p+1}\displaystyle\int_0^{\infty}x^pe^{-ax}dx\quad\cdots\,\)(10)

よって

\(\quad\displaystyle\int_0^{\infty}x^pe^{-ax}dx=\dfrac{p\,!}{a^{p+1}}\quad\cdots\,\)(11)

これを使って(7)式を積分すると

\(\quad\displaystyle\int_0^{\infty}\!\dfrac{x^3}{e^x-1}dx=\displaystyle\int_0^{\infty}\!x^3(e^{-x}+e^{-2x}+\cdots)dx=\dfrac{3!}{1^4}+\dfrac{3!}{2^4}+\cdots\quad\)

ここでゼータ関数より

\(\quad\zeta(4)=\displaystyle\sum_{n=1}^{\infty}\dfrac{1}{n^4}=\dfrac{1}{1^4}+\dfrac{1}{2^4}+\cdots=\dfrac{\pi^4}{90}\quad\cdots\,\)(12)

なので

\(\quad\displaystyle\int_0^{\infty}\dfrac{x^3}{e^x-1}dx=\dfrac{3!}{1^4}+\dfrac{3!}{2^4}+\cdots=3\,!\,\zeta(4)=\dfrac{\pi^4}{15}\quad\cdots\,\)(13)

4.放射エネルギー
全立体角に対する積分は
\(\quad\displaystyle\int\cos\theta\,d\Omega=\displaystyle\int_0^{2\pi}\displaystyle\int_0^{\frac{\pi}{2}}\cos\theta\,\sin\theta\,d\theta\,d\phi=\pi\quad\cdots\,\)(14)

よって黒体の単位面積当たりの放射エネルギー\(\,I\,\)は

\(\quad I=\pi\displaystyle\int B_{\nu}d\nu=\dfrac{2\pi^5k^4}{15c^2h^3}T^4=\sigma T^4\quad\cdots\,\)(15)

ここで、\(\,\sigma\,\)はステファン・ボルツマン係数である

\(\quad\sigma=\dfrac{2\pi^5k^4}{15c^2h^3}=5.67\times 10^{-8}\)W/m\(^2\)K\(^4\quad\cdots\,\)(16)

1. あらゆる物体は光を放っている
 電磁波というのは電荷を持った粒子が振動することで発生する。そして電荷を持った粒子は、どこからか飛んできた電磁波を受けると、その影響で振動する。そこらじゅうにある原子は内部に電荷を持っているのだから、何らかの方法で電磁波を放出し、また受け止めていると考えられる。

例えば太陽の光だって電磁波の一種だが、それを受け止めた物体は熱くなる。電磁波の形で遠くから届いたエネルギーが熱に変わるからだ。逆のことも起きている。キャンプファイヤーが終わった後で火を消すと辺りは真っ暗になるが、その真っ暗闇の中で、さっきまで焼かれていた石がぼんやりと赤黒く光っていることがある。そこらにある普通の石が、内部から光を放っているように見えて、少し不気味ではあるが、何だか綺麗だ。

他にも色んな例がある。ストーブの中で焼けた鉄板が赤く光っている。焚き火のときに炭が赤く光っている。ハロゲンヒーターが赤熱している。その光を浴びると、すぐに暖かく感じる。これは間にある空気を伝わって暖められるわけではない。光によって直接暖められるせいだ。このように、熱せられた物体から光や電磁波が出る現象を「熱放射」と呼ぶ。

では物体は何度くらいから光を出し始めるのだろう。実は目に見えないだけで、低い温度の物体でも電磁波を出しているらしい。体温の高い人が近くに来ると直接触らなくても暖かさを感じることがあるだろう。それはその人の体から、可視光線よりも低いエネルギーの光である赤外線が放射されているのだ。

2. ステファン・ボルツマンの法則
では物体は一体どれほどのエネルギーを熱放射の形で放出しているのだろうか？それを調べたのはヨーゼフ・ステファンという 19 世紀の科学者であり、実はあのボルツマンの師である。彼は他の学者たちが過去に行った幾つかの実験結果をまとめた上で、自らも白金の針金に電流を流して赤熱させる実験を行ったりして、熱放射のエネルギー\(\,K\,\)と物体の温度\(\,T\,\)との間に次のような関係があるのではないかという提案を行った。 (1879)

\(\quad K=\sigma T^4\quad\cdots\,\)(1)

この式は「ステファン・ボルツマンの法則」と呼ばれている。なぜボルツマンの名前が入っているかというと、後に、弟子のボルツマンが、当時の最新理論である電磁気学と熱力学を使って、この法則の根拠を説明することに成功したからである。 (1884)

にもかかわらず、この結果はしばらくの間認められることはなかった。多くの科学者がこの実験結果を確認しようとしたが、追試がうまくいかなかったのである。今ではその理由ははっきり分かっているのだが、この法則には少しばかり適用条件があって、当時はそのことがまだ良く分かっていなかったのである。ステファンの実験がうまく行ったのは偶然の幸運であろうと言われている。

3. 理論的導出
 電磁波はエネルギーだけでなく運動量をも持つのだということを、電磁気学で学んだ。そこで、こんなことを考えてみよう。容器の壁が温度\(\,T\,\)の物質で出来ていて、そこから放出された電磁波で容器の中が一様に満たされており、そのエネルギー密度が\(\,u\,\)になっていたとする。

 この仮定は次のような考えを根拠にしている。もし容器の壁から内側へ向けて電磁波が放出される一方だとすると、容器内部の電磁場のエネルギーは無限に増える一方となるだろう。しかしそうならないのは、容器の壁は電磁場を放出すると同時に吸収も行っていて、一定の状態で安定しているせいだと考えられる。容器の中の電磁場と容器の壁とは熱平衡に達しているというわけだ。電磁場に対しても、温度という概念が使えるということでもある。というわけで、\(\,u\,\)は\(\,T\,\)の関数だと言えるだろう。
さて、そのときの運動量密度\(\,w\,\)と\(\,u\,\)の間には\(\,u=cw\,\)という関係がある。

4. 電磁波のエネルギー密度
 マックスウェルの方程式から次の波動方程式が導かれる
\(\quad\nabla^2\mathbf{B}-\dfrac{1}{c^2}\dfrac{\partial^2\mathbf{B}}{\partial t^2}=0\qquad\nabla^2\mathbf{E}-\dfrac{1}{c^2}\dfrac{\partial^2\mathbf{E}}{\partial t^2}=0\quad\cdots\,\)(2)

空間に領域Dをとり、その内部で電磁場が荷電粒子系に及ぼす単位体積当たりの力は
\(\quad\mathbf{f}=\rho\mathbf{E}+\mathbf{j}\times\mathbf{B}\quad\cdots\,\)(3)
であるので荷電粒子系の運動方程式は
\(\quad\dfrac{d\mathbf{P}_{\text{粒子}}}{dt}=\displaystyle\int_{\mathrm{D}}\mathbf{f}\,dv\quad\cdots\,\)(4)
となる。\(\mathbf{P}_{\text{粒子}}\,\)は粒子系の全運動量である。

ここで電場についてのガウスの法則と拡張されたアンペールの法則を使って\(\,\mathbf{f}\,\)を表すと
\(\mathbf{f}=\epsilon_0(\nabla\cdot\mathbf{E})\mathbf{E}-\epsilon_0\dfrac{\partial\mathbf{E}}{\partial t}\times\mathbf{B}+\dfrac{1}{\mu_0}(\nabla\times\mathbf{B})\times\mathbf{B}\quad\cdots\,\)(5)

これよりD内の電磁場の運動量は単位体積当たり以下の値を持つ
\(\quad\mathbf{p}_{\text{場}}\equiv\epsilon_0\mathbf{E}\times\mathbf{B}\quad\cdots\,\)(6)

ここで電束密度\(\,\mathbf{D}=\epsilon_0\mathbf{E}\)、磁束密度\(\,\mathbf{B}=\mu_0\mathbf{H}\,\)であり、真空中を伝搬する電磁場の電場\(\,\mathbf{E}\)と磁束密度\(\,\mathbf{B}\)の大きさの間には\(\,|\mathbf{E}|=c|\mathbf{B}|\)なる関係が成立する。

電場の大きさを\(\,E=|\mathbf{E}|\)とすると電磁場のエネルギー密度は
\(\quad u=\dfrac{1}{2}(\mathbf{E}\cdot\mathbf{D}+\mathbf{B}\cdot\mathbf{H})=\dfrac{1}{2}\left(\epsilon_0E^2+\dfrac{E^2}{c^2\mu_0}\right)\)
\(\qquad =\dfrac{1}{2}(\epsilon_0E^2+\epsilon_0E^2)=\epsilon_0E^2\quad\cdots\,\)(7)

一方、運動量密度\(\,w\,\)は
\(\quad w=|\mathbf{D}\times\mathbf{B}|=\epsilon_0E\cdot\dfrac{E}{c}=\dfrac{\epsilon_0E}{c}=\dfrac{u}{c}\quad\cdots\,\)(8)

5. 電磁波による圧力

 ここで、容器の内部の電磁波を考える。電磁波は容器の壁に当たるたびに壁に運動量を与えるのだから、容器の壁を圧力\(\,p\,\)で押すことだろう。まずは、その\(\,p\,\)と\(\,u\,\)の関係を導く。

説明が簡単に済むように、一辺の長さが\(\,L\,\)の立方体の容器を考えてみる。そしてどの壁を考えても同じことではあるのだが、極座標の積分計算を分かりやすくする都合で、\(\,z\,\)軸の正の方向にある壁に与える圧力を考える。

容器内の全ての電磁波が一つの方向を向いているとする。その全運動量\(\,P\,\)は運動量密度\(\,w\,\)と体積\(\,V=L^3\,\)を掛けて、\(P=wL^3\,\)と表せる。

全ての電磁波の向かう一つの方向というのが、\(\,z\,\)軸に対して\(\,\theta\,\)の角度を持った方向だとすると、壁に当たって反射するたびに、合計で\(\,2P\cos\theta\,\)の運動量を与えることになる。

その頻度であるが、電磁波の速度の\(\,z\,\)軸方向成分が\(\,c\,\cos\theta\,\)であり、往復\(\,2L\,\)の距離を進むたびに再び同じ壁に当たるのだから、1秒間に\(\,c\,\cos\theta/2L\,\)回の衝突を行うことになる。

力というのは1 秒あたりの運動量変化のことであり、それを面積\(\,L^2\,\)で割れば圧力\(\,p\,\)が求まることになる。

\(\quad p=2P\cos\theta\times\dfrac{c\,\cos\theta}{2L}\div L^2\)
\(\qquad =2wL^3\cos\theta\times\dfrac{c\cos\theta}{2L^3}\)
\(\qquad =wc\,\cos^2\theta=u\,\cos^2\theta\quad\cdots\,\)(9)

という形になる。しかしこの結果は、全ての電磁波が一つの方向を向いていると仮定した上での話であり、実際にはあらゆる方向を向いているはずだから、全方向について平等に平均をしてやらないといけないだろう。その為には立体角の考えを使えば良さそうだ。自分が半径1の球の中心にいて、周囲をぐるっと見回していると想像してみよう。どの方向も等しいのだから、この球の表面積にわたって今の結果を積分してやればいい。平均を取るためにはそれを表面積\(\,4\pi\,\)で割ってやればいい。

\(\quad p=\dfrac{u}{4\pi}\displaystyle\int_0^{\pi}\cos^2\theta\sin\theta\,d\theta\displaystyle\int_0^{2\pi}d\phi\)
\(\qquad =\dfrac{u}{4\pi}\left[-\dfrac{1}{3}\cos^3\theta\right]_0^{\pi}\times 2\pi\)
\(\qquad =\dfrac{u}{3}\quad\cdots\)(10)

6. 熱力学的考察
 可逆過程での内部エネルギー(\(\,U\,\))変化において以下の関係式がある。
\(\quad dU=T\,dS-p\,dV\quad\) これを変形して
\(\quad\left(\dfrac{\partial U}{\partial V}\right)_T=T\left(\dfrac{\partial p}{\partial T}\right)_V-p\quad\cdots\,\)(11)

という式を作ることが出来る。これに\(\,T\,\)一定の条件下で\(\,V\,\)を変化させ、マクスウェルの関係式を当てはめる。電磁波の場合、この左辺はエネルギー密度\(\,u\,\)になる。

圧力\(\,p\,\)に式(10)を代入すると式(11)は
\(\quad u=\dfrac{1}{3}\,T\,\dfrac{du}{dT}-\dfrac{u}{3}\\[4pt]\therefore\,4u=T\,\dfrac{du}{dT}\qquad\dfrac{du}{u}=4\,\dfrac{dT}{T}\\[4pt]\quad\log u=4\log T+C\qquad\therefore\,\,u=aT^4\quad\cdots\text{(12)}\)となり、形の上ではステファン・ボルツマンの法則と同じものが導かれたことになる。

しかし\(\,u\,\)というのは、物体と電磁場とが熱平衡にあるときの電磁場のエネルギー密度を示しているのであって、物体から放出されるエネルギー\(\,K\,\)とは少し違う概念である。

7. 概念の整理
\(\,K\,\)の意味を「物体の単位表面積から単位時間あたりに放出されるエネルギー」だと言えるようにまとめることが出来れば使いやすい法則になりそうである。
微小面積\(\,dS\,\)から放出される電磁波をイメージする。
その微小面の法線方向から\(\,\theta\,\)だけずれた方向へ向かって微小時間\(\,dt\,\)に出て行った電磁波は、図のような筒状体積の内部にあることになるだろう。その体積は\(\,c\,\cos\theta\,dtdS\,\)である。

しかしこの筒状体積とエネルギー密度\(\,u\,\)をかけるだけで、そのエネルギー量が求まるほど単純ではない。なぜなら、このエネルギー密度\(\,u\,\)というのは、あらゆる方向に向かう電磁波が均等に重なり合った合計のエネルギー密度を表しているからである。

今知りたいのは、\(\,\theta\,\)方向へ流れる電磁波のエネルギーだけなのだが、その割合は全方向に向かう電磁波のエネルギーの中でどれくらいあると言えるだろうか。いや、厳密に\(\,\theta\,\)方向に向かう電磁波だけに限定してしまうとそれは無に等しいのである。だからある程度の微小幅を許容して考えないといけない。

こういう時は立体角の考えを使う。半径1の球殻を考え、その表面積の比で表すのである。角度\(\,d\theta\,\)と\(\,d\phi\,\)によって張られる球殻上の微小表面積は\(\,sin\theta\,d\theta\,d\phi\,\)と表されるから、全球殻の面積\(\,4\pi\,\)と比べた割合は\(\,1/4\pi\,\sin\theta\,d\theta\,d\phi\,\)となるだろう。

先ほど考えた筒状領域を、厳密に一方向として考えるのではなく、角度\(\,d\theta\,\)と\(\,d\phi\,\)の微小角度の幅を許してやれば、その範囲を向いている電磁波のエネルギーは

\(\quad\dfrac{1}{4\pi}uc\,\cos\theta\sin\theta\,d\theta\,d\phi\,dt\,dS\quad\cdots\,\)(13)

と表せることになる。この式(13)を微小面より上にある全ての方向について積分してやれば、微小面\(\,dS\,\)から微小時間\(\,dt\,\)内にあらゆる方向へ出て行くエネルギーを表せる。

\(\quad\dfrac{1}{4\pi}uc\,dt\,dS\displaystyle\int_0^{\pi/2}\cos\theta\sin\theta\,d\theta\displaystyle\int_0^{2\pi}d\phi\\[3pt]\qquad=\dfrac{1}{2}uc\,dt\,dS\displaystyle\int_0^{\pi/2}\dfrac{1}{2}\sin 2\theta\,d\theta\\[3pt]\quad =\dfrac{1}{2}uc\,dt\,dS\left[-\dfrac{1}{4}\cos 2\theta\right]_0^{\pi/2}\\[3pt]\quad =\dfrac{1}{4}uc\,dt\,dS\quad\cdots\,\text{(14)}\)

意味を考えれば今計算したのは\(\,K\,dt\,dS\,\)のことでもある。つまり\(\,K\,\)というのは、

\(\quad K=\dfrac{c}{4}u\quad\cdots\,\)(15)
に相当することが言えるのである。ここに先ほどの\(\,u=aT^4\,\)を代入すれば、

\(\quad K=\dfrac{ca}{4}\,T^4=\,\sigma T^4\quad\cdots\,\)(16)

という形が成り立っていることが言えるのである。普通はステファン・ボルツマンの法則と言えば、ここで使ったような意味での\(\,K\,\)と温度\(\,T\,\)の関係のことであり、「ステファン・ボルツマンの定数」として良く知られた\(\,\sigma\,\)もその意味に合うような数値が書かれていることが多い。

放送大学での天文での研究として、恒星について報告書（修論）を出したんですが、結果が中途半端でした。連星での質量のやりとりをH\(\alpha\)線のスペクトルを観測して調べようというものです。近接連星なので片方（質量の小さい）から大量の水素が流れ込んでいます。水素の温度は1万\(^{\circ}\)Kを超えており、もう片方の星に落ちるときはそのエネルギーで、H\(\alpha\)線のスペクトルが輝線または吸収線として現れます。
表面近くからスペクトル線が観測されると予想していたら、両方の共通重心の位置からしか観測されません。これを少しでも理解するためには、物理の知識が必要です。
高温なので当然、プラズマになっているので形としては「電磁流体解析」となるようです。まずは、流体力学を勉強しないとということで、入門書から始めます。
テキストとして、日本機械学会のJSMEシリーズの「流体力学」用います。

レーン=エムデン方程式は以下の形で表される

\(\quad\dfrac{1}{\xi^2}\dfrac{d}{d\xi}\left(\xi^2\dfrac{d\theta}{d\xi}\right)=-\theta^N\qquad\cdots\,\)(1)

境界条件は、\(\xi\to 0\,\)のとき\(\,\theta=1\,,\,\dfrac{d\theta}{d\xi}=0\,\)である。

［１］中心近傍での扱い
この方程式は\(\xi=0\,\)が確定特異点になっている。この場合、特異点近傍でローラン級数展開を考える。\(\xi=0\,\)の周りの形式解を以下のようにする。

\(\quad\theta(\xi)=1+a_2\xi^2+a_4\xi^4+\cdots\qquad\)(2)

式(2)を式(1)の左辺に代入する。まず、\(\xi\,\)で微分すると

\(\quad\dfrac{d\theta}{d\xi}=2a_2\xi+4a_4\xi^3+\cdots\qquad\)(3)

さらに\(\xi^2\,\)をかけて微分すると

\(\quad\dfrac{d}{d\xi}\left(\xi^2\dfrac{d\theta}{d\xi}\right)=6a_2\xi^2+20a_4\xi^4+\cdots\qquad\)(4)

\(\therefore\,\,\)左辺\(\,\,=6a_2+20a_4\xi^2+\cdots\qquad\)(5)

右辺は、まず\(B_2=a_2\xi^2\,,\,B_4=a_4\xi^4\,,\cdots\,\)として、

\(\quad(1+B_2+B_4+\cdots)^N\,\)の展開を考える。

最初の\(\,1\,\)の項は\(\,1^NB_2^0B_4^0\cdots\,\)であるので、

係数は\(\,\dfrac{N!}{N!0!0!\cdots}\,\)となる。

次の項は\(1^{N-1}B_2^1B_4^0\cdots\,\)で、

係数は\(\,\dfrac{N!}{(N-1)!1!0!\cdots}\,\)となり、

次の\(1^{N-1}B_2^0B_4^1B_6^0\cdots\,\)の係数は

\(\,\dfrac{N!}{(N-1)!0!1!0!\cdots}\,\)となる。因って

\(\quad(1+B_2+B_4+\cdots)^N=1+NB_2+NB_4+\cdots\qquad\)(6)

右辺\(\,=-\theta^N=-1-Na_2\xi^2-Na_4\xi^4-\cdots\,\)となるので

\(\,\xi^m\,\)の係数を比較すると

\(\quad 6a_2=-1\,,\qquad 20a_2=-Na_4\qquad\)となり

\(\quad a_2=-\dfrac{1}{6}\,,\qquad a_4=\dfrac{N}{120}\quad\)を得る。

よって中心近傍での解は

\(\quad \theta(\xi)\simeq 1-\dfrac{1}{6}\xi^2+\dfrac{N}{120}\xi^4\quad\cdots\,\)(7)

と表すことが出来る。

［２］数値計算
この方程式は(その2)で求めた解析解以外の解は数値計算に因らざる得ない。
\(\bigstar\,\)常微分方程式の解法
\(\quad\dfrac{dy_i}{dx}=f_i[x,y_i(x)]\quad i=1,2,\dots\quad\)(8)

を初期値問題として数値的に解く。
まず、\(x=x_0\,\)において初期値\(\,y_i(x_0)=y_{i0}\,\)が与えられたとき、\(x=x_0+h\,\)での\(\,y_i\,\)の値\(\,y_i(x_0+h)\,\)をもとめ、次のステップでその\(\,(x,y_i)\,\)を新たな初期値として処理を繰り返す。
\(x_0\,\)と\(\,x_0+h\,\)の間の中間の\(\,x\,\)における導関数、\(f(x,y)\,\)の関数値を何回か計算し、それらの適当な平均を用いることによって、高次のTaylor展開を達成するのが Runge-Kutta 法である。

\(\bigstar\,\)Runge-Kutta 法
\(\quad k_1=hf(x_0,y_0)\qquad\qquad\qquad\cdots\,\)(9)
\(\quad k_2=hf(x_0+\frac{h}{2},y_0+\frac{k_1}{2})\qquad\cdots\,\)(10)
\(\quad k_3=hf(x_0+\frac{h}{2},y_0+\frac{k_2}{2})\qquad\cdots\,\)(11)
\(\quad k_4=hf(x_0+h,y_0+k_3)\qquad\cdots\,\)(12)
\(\quad y(x_0+h)=y(x_0)+\frac{1}{6}[K_1+2(K_2+k_3)+k_4]\,\cdots\,\)(13)

まず点(\(\,x_i,y_i\,\))での微係数から\(\,y\,\)の増分\(\,k_1\,\)を推定する(式9)。次に(\(\,x_i,y_i\,\))と今得た点(\(\,x_i+\Delta x,y_i+k_1\,)\)の中間点での微係数から、これに平行な破線で示す直線によって増分\(\,k_2\,\)の推定を行う(式10)。これより第３の推定点\((\,x_i+\Delta x,y_i+k_3\,)\)を得る。第４の推定点まで求めたところで、これら４つの増分推定を\(\,1:2:2:1\,\)の重みを付けて平均化し、これを最終的な増分推定として\(\,y_{i+1}\,\)を求める(式13)。この式は Taylor展開の４次の精度を持つ。

\(\bigstar\,\) 初期値とプログラミング
解くべき式は(1)式だが、\(\xi=x\,,\,\theta=y\,\)と書き換えて見やすくする。

\(\quad\dfrac{1}{x^2}\dfrac{d}{dx}\left(x^2\dfrac{dy}{dx}\right)=-y^N\quad x\!=\!0\quad y\!=\!0\,,\,\dfrac{dy}{dx}=0\,\cdots\,\)(14)
この常微分方程式を連立１階微分方程式に変換する。

\(y1=y\,,\quad y2=x^2\dfrac{dy}{dx}\,\)とおくと
\(\quad \dfrac{dy1}{dx}=\dfrac{y2}{x^2}\qquad\cdots\,\)(15)
\(\quad \dfrac{dy2}{dx}=-x^2y1^N\qquad\cdots\,\)(16)

初期値としては\(\,xi=x=0\,\)でなく、演算の刻み\(\,\xi=x=h\,\)とし、\(y1\,\)は式(7)を使い
\(\quad y1=1-\dfrac{1}{6}h^2+\dfrac{N}{120}h^4\quad\cdots\,\)(17)

\(y2\,\)は、\(\xi^2\dfrac{d\theta}{d\xi}=y2\quad\)としたので
\(\quad y2=h^2\left(-\dfrac{1}{3}h+\dfrac{N}{30}h^3\right)=-\dfrac{1}{3}h^3+\dfrac{N}{120}h^4\quad\cdots\,\)(18)
で設定する。
計算は中心\(\,\,\xi=x=h\,\,\)から始めて、密度係数がゼロ\(\,\,\,\theta=y1=0\,\,\)となる表面までを積分する。また\(\,\,h=1.0\times10^{-6}\,\,\)とした。

［３］数値計算による結果
\(N\,\)を\(0\sim 4.5\,\)まで計算した結果をグラフと表に示す。

この表は、\(N=n\,\)でガス球の表面（圧力ゼロ\(\,\theta=0\,\)）の時の\(\,\xi_n\,\)とそこでの\(\,-1/\xi\,(d\theta/d\xi)\,\)で、中心密度\(\rho_c\,\)を計算するためのパラメータである。

\(\bigstar\,\)結果の見方
\(\xi\,\)は \(\,R=\alpha\xi=\left[\dfrac{(N+1)K}{4\pi G}\rho_c^{\frac{1}{N}-1}\right]^{1/2}\xi\,\cdots\)(21)　である。

また中心から\(\,\xi\,\)より内側に含まれる質量\(\,m(\xi)\,\)は
\(\quad m(\xi)=\displaystyle\int_0^{\alpha\xi}4\pi\rho r^2dr\)
\(\qquad\quad =4\pi\alpha^3\rho_c\displaystyle\int_0^{\xi}\xi^2\theta^Nd\xi\)
\(\qquad\quad =-4\pi\alpha^3\rho_c\displaystyle\int_0^{\xi}\dfrac{d}{d\xi}\left(\xi^2\dfrac{d\theta}{d\xi}\right)d\xi\)
\(\qquad\quad =-4\pi\alpha^3\rho_c\xi^2\dfrac{d\theta}{d\xi}\quad\cdots\,\)(22)
で与えられる。したがって星の全質量\(\,M\,\)は
\(\quad M=-4\pi\alpha^3\rho_c\left(\xi^2\dfrac{d\theta}{d\xi}\right)_{\xi=\xi_n}\)
\(\qquad =-4\pi\left[\dfrac{(N+1)K}{4\pi G}\rho_c^{\frac{1}{N}-1}\right]^{3/2}\left(\xi^2\dfrac{d\theta}{d\xi}\right)_{\xi=\xi_n}\,\cdots\,\)(23)
ここで星の平均密度\(\,\bar{\rho}\,\)を考える。
\(\quad\bar{\rho}=\dfrac{M}{4\pi R^3/3}=\dfrac{3M}{4\pi R^3}\)
\(\qquad =-\dfrac{3}{4\pi\alpha^3\xi_n^3}\cdot4\pi\alpha^3\rho_c\left(\xi^2\dfrac{d\theta}{d\xi}\right)_{\xi=\xi_n}\)
\(\qquad =-\dfrac{3}{\xi_n}\rho_c\left(\dfrac{d\theta}{d\xi}\right)_{\xi=\xi_n}\quad\cdots\,\)(24)

\(\bigstar\,\)太陽の場合
これらの関係を太陽の場合で計算してみる。太陽の半径\(\,R=6.960\times10^8\)m 、質量\(\,M=1.9884\times10^{30}\)kg と、恒星の場合に適応できるポリトロープ指数\(\,N=3\,\,(\Gamma=\frac{4}{3}\,)\,\)の計算結果を用いる。
\(\quad\rho_c=\dfrac{M}{4\pi R^3}\left(-\dfrac{1}{\xi}\dfrac{d\theta}{d/xi}\right)_{\xi=\xi_n}^{-1}\quad\cdots\,\)(25)

\(\qquad=\dfrac{1.9884\times10^{30}}{4\pi\cdot (6.960\times10^8)^3\cdot6.152\times10^{-3}}\)

\(\qquad =7.629\times10^4\,\mathrm{kg}\cdot\mathrm{m}^{-3}\)
理科年表によると\(\,\rho_c=1.56\times10^5\,\mathrm{kg}\cdot\mathrm{m}^{-3}\,\)*となっており、計算値は約半分である。まあ、天文のざっくりモデルでの計算としては、まあ合っているといえる。
*[参考]理科年表の太陽の構造の部分

その１で導出したレーン=エムデン方程式は
\(\quad\dfrac{1}{\xi^2}\dfrac{d}{d\xi}\left(\xi^2\dfrac{d\theta}{d\xi}\right)=-\theta^N\quad\cdots\,\)(1)
で、\(\,\theta\,\)は無次元化した密度、\(\,\xi\,\)は無次元化した半径である。中心での境界条件は、\(\xi\to 0\,\)のとき\(\,\theta=1\,,\,\dfrac{d\theta}{d\xi}=0\,\)である。

［１］\(N=0\,\)の時の解析解
式(1)を一般の\(\,N\,\)について解析解は得られないが、\(N=0\,,\,1\,,\,5\,\)の時は得ることが出来る。
なお、以降は見やすくするために\(\,\xi\,,\,\theta\,\)を\(\,x\,,\,y\,\)と記述する。式(1)は
\(\quad\dfrac{1}{x^2}\dfrac{d}{dx}\left(x^2\dfrac{dy}{dx}\right)=-x^N\quad\cdots\,\)(2)
となる。ここに\(\,N=0\,\)を代入すると
\(\quad\dfrac{1}{x^2}\dfrac{d}{dx}\left(x^2\dfrac{dy}{dx}\right)=-1\quad\cdots\,\)(3)
両辺に\(\,x^2\,\)をかけると
\(\quad\dfrac{d}{dx}\left(x^2\dfrac{dy}{dx}\right)=-x^2\quad\cdots\,\)(4)

この式の両辺を\(\,x\,\)で積分すると
\(\quad x^2\dfrac{dy}{dx}=-\dfrac{1}{3}x^3+C_1\,\cdots\,\)(5)
となる。\(C_1\,\)は積分定数。両辺を\(\,x^2\,\)で割ると
\(\quad\dfrac{dy}{dx}=-\dfrac{1}{3}x+\dfrac{C_1}{x^2}\quad\cdots\,\)(6)

再度\(\,x\,\)で積分すると
\(\quad y=-\dfrac{1}{6}x^2-\dfrac{C_1}{x}+C_2\quad\cdots\,\)(7)
ここで、\(\,x=0\,\)のとき\(\,y=1\,\)の条件を満たすためには、\(C_1=0\,\)で\(\,C_2=1\,\)とする必要がある。よって、\(N=0\,\)の時は
\(\quad y=1-\dfrac{1}{6}x^2\quad\cdots\,\)(8)
となる。

［２］\(N=1\,\)の時の解析解
式(2)に\(\,N=1\,\)を代入すると
\(\quad\dfrac{1}{x^2}\dfrac{d}{dx}\left(x^2\dfrac{dy}{dx}\right)+y=0\quad\cdots\,\)(9)
ここで、\(z=xy\,\)とすると
\(\quad\dfrac{dy}{dx}=\dfrac{d}{dx}\left(\dfrac{z}{x}\right)=\dfrac{1}{x}\dfrac{dz}{dx}-\dfrac{z}{x^2}\quad\cdots\,\)(10)

式(10)の両辺に\(\,x^2\,\)をかけると
\(\quad x^2\dfrac{dy}{dx}=x\dfrac{dz}{dx}-z\quad\cdots\,\)(11)
となり、これを\(\,x\,\)で微分すると
\(\quad\dfrac{d}{dx}\left(x^2\dfrac{dy}{dx}\right)=\dfrac{d}{dx}\left(x\dfrac{dz}{dx}-z\right)=x\dfrac{d^2z}{dx^2}\quad\cdots\,\)(12)
式(12)を式(9)に代入すると
\(\quad\dfrac{1}{x}\dfrac{d^2z}{dx^2}+y=\dfrac{d^2z}{dx^2}+z=0\quad\cdots\,\)(13)
という2階微分方程式になる。

この方程式の解は
\(\quad z=A\cos x+B\sin x\)
で、\(z=xy\,\)と変数を元に戻すと
\(\quad y=\dfrac{A\cos x}{x}+\dfrac{B\sin x}{x}\quad\cdots\,\)(14)

ここで、\(x\to 0\,\)の時、\(y\to 1\,\)なので、\(A=0\,,B=1\,\)となり
\(\quad y=\dfrac{\sin x}{x}\quad\cdots\,\)(15)
となる。

［３］\(N=5\,\)の時の解析解
式(2)に\(\,N=5\,\)を代入すると
\(\quad\dfrac{1}{x^2}\dfrac{d}{dx}\left(x^2\dfrac{dy}{dx}\right)=-y^5\quad\cdots\,\)(16)

ここで、特殊な変数変換を行う。
\(\quad x=e^{-t}\quad\cdots\,\)(17)
\(\quad y=\sqrt{\dfrac{1}{2}e^t}\,Z\quad\cdots\,\)(18)
式(17)より
\(\quad\dfrac{dx}{dt}=-x\quad\cdots\,\)(17)

また、\(Z\,\)の\(\,x\,\)微分は
\(\quad\dfrac{dZ}{dx}=\dfrac{dZ}{dt}\dfrac{dt}{dx}=-\dfrac{1}{x}\dfrac{dZ}{dt}\quad\cdots\,\)(18)

微分演算子には以下の関係がある
\(\quad\dfrac{d}{dx}=-\dfrac{1}{x}\dfrac{d}{dt}\quad\cdots\,\)(19)

まず、式(18)を\(\,x\,\)で微分すると
\(\quad\dfrac{dy}{dx}=-\dfrac{1}{2\sqrt{2}}x^{-\frac{3}{2}}Z+\dfrac{1}{\sqrt{2}}x^{-\frac{1}{2}}\dfrac{dZ}{dx}\)
\(\qquad=-\dfrac{1}{2\sqrt{2}}x^{-\frac{3}{2}}Z-\dfrac{1}{\sqrt{2}}x^{-\frac{3}{2}}\dfrac{dZ}{dt}\quad\cdots\,\)(20)

式(20)の両辺に\(\,x^2\,\)をかけると
\(\quad x^2\dfrac{dy}{dx}=-\dfrac{1}{2\sqrt{2}}x^{\frac{1}{2}}Z-\dfrac{1}{\sqrt{2}}x^{\frac{1}{2}}\dfrac{dZ}{dt}\quad\cdots\,\)(21)

この式(21)を\(\,x\,\)で微分する
\(\quad\dfrac{d}{dx}\left(x^2\dfrac{dy}{dx}\right)=-\dfrac{1}{4\sqrt{2}}x^{-\frac{1}{2}}Z-\dfrac{1}{2\sqrt{2}}x^{\frac{1}{2}}\dfrac{dZ}{dt}\)
\(\qquad\qquad-\dfrac{1}{2\sqrt{2}}x^{-\frac{1}{2}}\dfrac{dZ}{dt}-\dfrac{1}{\sqrt{2}}x^{\frac{1}{2}}\dfrac{d}{dx}\left(\dfrac{dZ}{dt}\right)\)
\(\qquad=-\dfrac{1}{4\sqrt{2}}x^{-\frac{1}{2}}Z+\dfrac{1}{\sqrt{2}}x^{-\frac{1}{2}}\dfrac{d^2Z}{dt^2}\quad\cdots\,\)(22)

ここで
\(\quad -x^2y^5=-x^2\dfrac{1}{4\sqrt{2}}x^{-\frac{5}{2}}Z^5=-\dfrac{1}{4\sqrt{2}}x^{-\frac{1}{2}}Z^5\quad\cdots\,\)(23)
なので式(22)(23)を(16)に代入すると
\(\quad-\dfrac{1}{4}Z+\dfrac{d^2Z}{dt}=-\dfrac{1}{4}Z^5\quad\)より
\(\quad\dfrac{d^2Z}{dt^2}=\dfrac{1}{4}Z(1-Z^4)\quad\cdots\,\)(24)
という形に書き直すことができる。式(24)の両辺に\(\dfrac{dZ}{dt}\,\)をかけると
\(\quad\left(\dfrac{dZ}{dt}\right)\left(\dfrac{d^2Z}{dt^2}\right)=\dfrac{1}{4}Z\left(\dfrac{dZ}{dt}\right)-\dfrac{1}{4}Z^5\left(\dfrac{dZ}{dt}\right)\,\cdots\,\)(25)

ここで左辺は以下のように変形でき
\(\quad\left(\dfrac{dZ}{dt}\right)\left(\dfrac{d^2Z}{dt^2}\right)=\dfrac{1}{2}\dfrac{d}{dt}\left[\left(\dfrac{dZ}{dt}\right)^2\right]\quad\cdots\,\)(26)

右辺の計算も
\(\quad Z\left(\dfrac{dZ}{dt}\right)=\dfrac{1}{2}\dfrac{d}{dt}\left(Z^2\right)\quad\cdots\,\)(27)
と変形できるので、式(25)は
\(\quad\dfrac{d}{dt}\left[\left(\dfrac{dZ}{dt}\right)^2\right]=\dfrac{1}{4}\dfrac{d}{dt}\left(Z^2\right)-\dfrac{1}{12}\dfrac{d}{dt}\left(Z^6\right)\quad\cdots\,\)(28)

\(\,t\,\)で積分すると
\(\quad\left(\dfrac{dZ}{dt}\right)^2=\dfrac{Z^2}{4}-\dfrac{Z^6}{12}+D\quad\cdots\,\)(29)
となる。\(\,D\,\)は積分定数である。\(x\to 0\,\)のとき、\(Z=0\,dZ/dt=0\,\)なので、\(D=0\,\)である。よって次の式が得られる。

\(\quad\dfrac{dZ}{Z\sqrt{1-\frac{1}{3}Z^4}}=\pm\dfrac{1}{2}dt\quad\cdots\,\)(30)
左辺の正負の符号は境界条件より負を選択する。
\(\quad\dfrac{dZ}{Z\sqrt{1-\frac{1}{3}Z^4}}=-\dfrac{1}{2}dt\quad\cdots\,\)(31)

この方程式を解くために、さらなる変数変換を行う。
\(\quad\dfrac{1}{3}Z^4=\sin^2\phi\quad\cdots\,\)(32)
この式の微分形は以下のようになる。
\(\quad\dfrac{4}{3}Z^3\,dZ=2\sin\phi\,\cos\phi\,d\phi\quad\cdots\,\)(33)

式(32)(33)を(31)の左辺に代入すると
\(\quad\dfrac{dZ}{Z\sqrt{1-\frac{1}{3}Z^4}}=\dfrac{1}{Z\sqrt{1-\sin^2\phi}}\dfrac{3}{4Z^3}2\sin\phi\,\cos\phi\,d\phi\)
\(\qquad=\dfrac{3}{2}\dfrac{1}{Z^4}\dfrac{\sin\phi\,\cos\phi\,d\phi}{\sqrt{1-\sin^2\phi}}=\dfrac{1}{2\sin^2\phi}\dfrac{\sin\phi\,\cos\phi}{\cos\phi}d\phi\)
\(\qquad=\dfrac{1}{2}\dfrac{1}{\sin\phi}\,d\phi=\dfrac{1}{2}\mathrm{cosec}\,\phi\,d\phi\quad\cdots\,\)(34)
より(31)の方程式は以下の形に変形できる。

\(\quad\mathrm{cosec}\,\phi=-dt\quad\cdots\,\)(35)
両辺を積分すると
\(\quad\ln\left|\tan\dfrac{\phi}{2}\right|=-t+C’\quad\cdots\,\)(36)
より
\(\quad\tan\dfrac{\phi}{2}=Ce^{-t}\quad\cdots\,\)(37)

ここで、式(32)から
\(\quad Z^4=3\sin^2\phi=\dfrac{12\tan^2\frac{\phi}{2}}{\left(1+\tan^2\frac{\phi}{2}\right)^2}\quad\cdots\,\)(38)
この式(38)に式(37)を代入すると
\(\quad Z^4=\dfrac{12C^2e^{-2t}}{(1+C^2e^{-2t})^2}\)
\(\quad Z=\pm\left[\dfrac{12C^2e^{-2t}}{(1+C^2e^{-2t})^2}\right]^{\frac{1}{4}}\quad\cdots\,\)(39)

ここで、\(Z\,\)は式(18)より正で、式(17)(18)を代入すると
\(\quad y=\left[\dfrac{3C^2}{(1+C^2e^{-2t})^2}\right]^{\frac{1}{4}}=\left[\dfrac{3C^2}{(1+C^2x^2)^2}\right]^{\frac{1}{4}}\quad\cdots\,\)(40)

\(x=0\,\)で\(\,y=0\,\)という条件なので、
\(\quad 1=(3C^2)^{\frac{1}{4}}\quad\)より\(\quad C^2=\dfrac{1}{3}\quad\)が得られ
\(\quad y=\left[\dfrac{1}{(1+\frac{1}{3}x^2)^2}\right]^{\frac{1}{4}}=\dfrac{1}{\sqrt{1+\frac{1}{3}x^2}}\quad\cdots\,\)(41)
となる。

やってみたかった事

・高校では『天文研』に所属して、太陽の黒点の観測や夏の奥多摩での観測合宿などやってそれなりに、天文してました。完成はしませんでしたが、放課後交代で反射鏡を磨いていました。その後、大学と生業は『化学』という天文とは多少距離がある分野で長らく過ごしていました。
・やってみたかったと言うほどではないですが、いつかは．．．的な思いがあったのは、確かです。

放送大学で

・一応、定年（６０歳）後を見据えて５０代から、放送大学の専科履修生として勉強していました。定年後に放送大学の「自然と環境コース」に全科履修生として入学しました。(2013年度)
・それまで、専科履修生で取得した単位があったので、3年からの編入としました。昔の大学の単位を認定してくれる制度はあったようですが、面倒だと考えて語学も含めて全ての単位を放送大学で取得することにしました。その夏には卒業研究の申請を提出しています。
・それまで学習センターは試験や面接授業でしか利用していませんでしたが、いろいろなサークルなどがあり、先輩諸兄が指導してくれる環境がありがたかったです。卒論はどうせなら\(\,\mathrm{\TeX}\,\)を使いたいと思っていたところ、サークルにて\(\,\mathrm{\TeX}\,\)の講座を開催してくれたので、渡りに船と勉強させてもらいました。
・自然科学系の卒研に必要な科目と言ったら、まず数学です。一応、工学部だったのですが、昔取った杵柄が全く役に立ちませんでした。その数学を親切に教えてくださる先輩がいて、初歩から(再?)勉強させてもらっています。

天文での卒研

・卒業研究は天文を選択しました。報告書（卒論）はもちろん\(\,\mathrm{\TeX}\,\)で記述しました。放送大学のWAKABAの(資料)の「卒業研究資料_2014年度_06_自然と環境」に掲載されています。テーマはKeplar望遠鏡の観測データを使用した「連星の黒点挙動」です。

\(\,\mathrm{\TeX}\,\)について

・そんな感じで、備考録では\(\,\mathrm{\TeX}\,\)を使える環境をということで、ロリポップのレンタルサーバーの上で WordPress を使っています。
参考までに MathJaxの表示設定はこちらです。
<script type=”text/javascript” async=”” src=”https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.1/MathJax.js?config=TeX-AMS_CHTML”></script>
この方式でなくともプラグインして［mathjax］(括弧は半角)と文頭で宣言すればよいらしいですが．．．今のところ、うまく行かないので、文頭にスクリプトを書いています。

ポリトロープについて

・恒星は天文の勉強の基本です。ポリトロープは恒星の入門コースとなっています。これから導かれるレーン=エムデン方程式は、恒星物理学の基本です。約6年前勉強した内容を忘れないようにということで、ここに残します。
・振り子は最近までちゃんとした計算を知りませんでした。そんなんで、勉強した楕円積分とともに、備考録の最初に残しておきます。

［１］Lane Emden 方程式とは

・宇宙物理学や流体力学において、レーン＝エムデン方程式（Lane-Emden equation）は、球対称な密度分布を示す力学平衡にある自己重力流体を記述する微分方程式である。名称は宇宙物理学者のジョナサン・ホーマー・レーンとロバート・エムデンに由来する。
・Jonathan Homer Lane(1819-1880,ニューヨーク州のゲネシー生まれ)エール大学でオルムステッド（Denison Olmsted）に師事し、数学と自然科学を学ぶ。1857 年にワシントン市内に特許事務所を開き、特許関係の顧問や特許弁護士などの業務を始めた。1860 年に南北戦争が勃発すると、ペンシルバニア州北部の兄弟の住む村に移り、冷却用ガス膨張装置を用いた極低温装置の開発や天文学にも興味をもち、ひとりで文献に当たりながら太陽内部構造の研究に思いを深めていた。
・太陽の進化についてはジョン・ハーシェルやケルビン卿の収縮説に興味をもっていたが、内部構造を数学的に扱っていない点がレーンには不満であった。そこで彼は太陽をガス球と仮定し、重力平衡と熱力学過程を基本にして内部構造に関する基本方程式を書き始めた。
「太陽の理論的温度」(Lane,J.H.1870,Amer.J.Sci.&Arts,50,57-74) 自己重力のもとで球状に保たれるガス体の全体的平衡状態を取り扱った歴史的論文である。彼は太陽内部を完全気体と仮定したが、これは大胆な発想であった。完全気体を考えると\(\,p=R\rho T\)、静力学的平衡は\(\,dp=-\dfrac{Gm(r)}{r^2}\rho\,dr\,\)、内部は断熱的に対流平衡となるので、内部構造は\(\,1-\left(\dfrac{\rho}{\rho_c}\right)^{\gamma-1}=\displaystyle\int_0^x\dfrac{\mu(x)}{x^2}dx\,\)にて与えられる。

［２］アウグスト・リターのガス球論

・August Ritter(1826-1908)1853年にゲッチンゲン大学卒業。アーヘン工業大学(現RWTH-Aachen)の構造力学と一般力学の担当教授。「大気の高さに関する考察とガス状天体の内部構造について」ガス球の平衡・振動・安定性の問題に取り組み、自己重力にあるポリトロープガス球の平衡状態についての先駆的研究。
・ポリトロープとは圧力と密度に\(\,pV^n=K\,\) または\(\,p=K\rho^{1+1/n}\,\)の関係を満たす過程のこと。
\(\,T/T_c=y\,,\quad r/a=x\,\)とおいて、重力平衡式と組み合わせると以下の基本微分方程式が得られる。
\(\quad\dfrac{d^2y}{dx^2}+\dfrac{2}{x}\dfrac{dy}{dx}+y^n=0\,\quad\cdots\,\)(1)
この式は ポリトロープガス球の内部構造に関する基本方程式として、エムデンやエディントンに大きな影響を与えた。

［３］エムデンのガス球論

・ロベルト・エムデン(Jacob Robert Emden,1862-1940)スイスの天体物理学および気象学の専門家。1907 年にミュンヘン工科大学の教授。「ガス球論」出版。リターの基本方程式もそのままの形で再度導かれており、通常、エムデンの微分方程式と呼ばれている。リターらと異なり、エムデンはガス球半径を有限の場合と無限の場合を併行して考察している。エムデンの基本的考え方はレーンやリターと変わっていないが、理論が整備され、応用例が多いことからガス球理論の集大成と呼ばれ、4 版を重ねて、星の内部構造論の基本的テキストとなって広く知られるようになった。なお、エムデンの甥、マーチン・シュワルツシルドは 1958 に著書「星の構造と進化｣を著わしているが、この書で歴史的なマイルストーンの研究者として、順にエムデン、エディントン、チャンドラセカールの 3 人の名前を挙げている。

［４］星の力学的安定

・連続の式:
　半径\(\,r\,\)、ガス密度は半径の関数とし、\(\,\rho(r)\,\) と表すと、その半径より内側に含まれる質量\(\,M_r\,\)は
\(\quad M_r=\displaystyle\int_0^r4\pi r’^2\rho(r)dr’\quad\cdots\,\)(2)
と書くことができる。
この式の両辺を\(\,r\,\)で微分して微分形式で書き表すと
\(\quad\dfrac{dM_r}{dr}=4\pi r^2\rho(r)\quad\cdots\,\)(3)
となる。この式は流体力学の連続の式(continuity equation)の一種である。

・静水圧平衡:
　星は自分自身の重力による収縮しようとする力とガス圧による対抗する力がバランスして、力学平衡を保っている。
\(\,r\,\)を中心からの距離、\(\,P\,\)をそこでの圧力、\(\,G\,\)を重力定数として
\(\quad\dfrac{dP}{dr}=-\dfrac{GM_r}{r^2}\rho\quad\cdots\,\)(4)
と表せる。
連続の式(3)、静水圧平衡の式(4)を解くことで、天体の内部構造を明らかにすることができる。 しかしこの段階では、独立変数が半径\(\,r\,\)で、従属変数は密度\(\,\rho\)、圧力\(\,P\)、質量\(\,M_r\,\)の３つで方程式は２つしかない。ためもうひとつ変数を結びつける式が必要となる。

・ポリトロープ関係
ここで、前出のポリトロープ関係と呼ばれる関係式を導入する。
\(\quad P=K\rho^{1+\frac{1}{N}}\quad\cdots\,\)(5)
\(\,K\,,\,\,N\,\)は定数のパラメータであり、\(\,N\,\)のことをポリトロープ指数と呼ぶ。

［５］レーン=エムデン方程式の導出

天体内部の構造を決めるためには、式(3)、(4)、(5)を解けば良いことになる。まず式(4)から
\(\quad\dfrac{r^2}{\rho}\dfrac{dP}{dr}=-GM_r\quad\cdots\,\)(6)
となり、両辺を\(\,r\,\)で微分すると
\(\quad\dfrac{d}{dr}\left(\dfrac{r^2}{\rho}\dfrac{dP}{dr}\right)=-G\dfrac{dM_r}{dr}\quad\cdots\,\)(7)
となる。これを式(3)に代入すると
\(\quad\dfrac{d}{dr}\left(\dfrac{r^2}{\rho}\dfrac{dP}{dr}\right)=-4\pi Gr^2\rho\quad\cdots\,\)(8)
となる。また、式(5)を\(\,r\,\)で微分すると
\(\quad\dfrac{dP}{dr}=K\left(1+\dfrac{1}{N}\right)\rho^{\frac{1}{N}}\dfrac{d\rho}{dr}\quad\cdots\,\)(9)
が得られる。式(9)を式(8)に代入すると
\(\quad K\left(1+\dfrac{1}{N}\right)\dfrac{d}{dr}\left(r^2\rho^{\frac{1}{N}-1}\dfrac{d\rho}{dr}\right)=-4\pi Gr^2\rho\quad\cdots\,\)(10)

ここで密度\(\,\rho(r)\,\)を星の中心密度\(\,\rho_c\,\)として
\(\quad\rho=\rho_c\theta^N\quad\cdots\,\)(11)
の形で表すこととする。\(\theta\,\)は無次元化した密度。
両辺を\(\,r\,\)で微分すると
\(\quad\dfrac{d\rho}{dr}=\dfrac{d}{dr}\left(\rho_c\theta^N\right)=\rho_cN\theta^{N-1}\dfrac{d\theta}{dr}\quad\cdots\,\)(12)
また、式(11)から以下の関係が導かれる
\(\quad\rho^{\frac{1}{N}-1}=\rho_c^{\frac{1}{N}-1}\,\theta^{1-N}\quad\cdots\,\)(13)

式(10)に式(12)(13)を代入すると
\(\quad K\left(1+\dfrac{1}{N}\right)\dfrac{d}{dr}\left(r^2\rho_c^{\frac{1}{N}-1}\theta^{1-N}\rho_cN\theta^{N-1}\dfrac{d\theta}{dr}\right)\!=\!-4\pi Gr^2\rho_c\theta^N\quad\)
\(\quad\dfrac{K(N+1)\rho_c^{\frac{1}{N}-1}}{4\pi G}\dfrac{1}{r^2}\dfrac{d}{dr}\left(r^2\dfrac{d\theta}{dr}\right)=-\theta^N\quad\cdots\,\)(14)

また式(5)を星の中心にあてはめると
\(\quad P_c=K\rho_c^{1+\frac{1}{N}}\quad\cdots\,\)(15)
この式(15)を式(14)に代入して\(\,K\,\)を消去すると
\(\quad\dfrac{(N+1)P_c}{4\pi G\rho_c^2}\dfrac{1}{r^2}\dfrac{d}{dr}\left(r^2\dfrac{d\theta}{dr}\right)=-\theta^N\quad\cdots\,\)(16)

左辺の係数部分を以下のように定義すると
\(\quad\alpha\equiv\left[\dfrac{(N+1)K}{4\pi G}\rho_c^{\frac{1}{N}-1}\right]^{1/2}\quad\cdots\,\)(17)
\(\quad\dfrac{\alpha^2}{r^2}\dfrac{d}{dr}\left(r^2\dfrac{d\theta}{dr}\right)=-\theta^N\quad\quad\cdots\,\)(18)

ここで、半径を以下の形に無次元化する。
\(\quad\xi=\dfrac{r}{\alpha}\quad\cdots\,\)(19)
これより無次元化した密度\(\,\theta\,\)の\(\,r\,\)の微分は以下のようになる
\(\quad\dfrac{d\theta}{dr}=\dfrac{d\theta}{d\xi}\dfrac{d\xi}{dr}=\dfrac{1}{\alpha}\dfrac{d\theta}{d\xi}\quad\cdots\,\)(20)

微分演算子は以下の形になる
\(\quad\dfrac{d}{dr}=\dfrac{1}{\alpha}\dfrac{d}{d\xi}\quad\cdots\,\)(21)

式(18)にこれらを用いると
\(\quad\dfrac{1}{\xi^2}\dfrac{d}{d\xi}\left(\xi^2\dfrac{d\theta}{d\xi}\right)=-\theta^N\quad\cdots\,\)(22)
この式は無次元化した密度\(\,\theta\)、無次元化した半径\(\,\xi\,\)を結びつけるもので、
レーン=エムデン方程式と呼ばれている。

天文においてポリトロープとは、圧力と密度がポリトロピック関係を満たし力学平衡にある球対称な流体である。
一般に、圧力\(\,P\,\)と密度\(\,\rho\,\)の間に
\(\quad P=K\rho^{1+\frac{1}{N}}\quad\cdots\,\) (1)
の関係があるとき、この関係式をポリトロピックと呼ぶ。ここで\(\,N\,\)はポリトロピック指数、\(\,K\,\)は比例定数である。
熱力学では \(\,1+\dfrac{1}{N}=\gamma\,\) であるが、天文では後述の Lane-Emden 方程式の解法のために、この形の表記を行う。
また、断熱圧縮率(\(\,\Gamma_1\,\))を以下のような形で表す。
\(\quad\Gamma_1\equiv\,\left(\dfrac{\partial\ln P}{\partial\ln\rho}\right)_S=\left(\dfrac{\partial\ln P}{\partial\ln\rho}\right)_T\gamma\quad\cdots\,\) (2)
理想気体の場合\(\, (\partial\ln P/\partial\ln\rho)_T=1\,\)となるので、
\(\quad\Gamma_1=\gamma\equiv\,\dfrac{C_P}{C_V}\quad\cdots\,\) (3)
ポリトロピック指数と星の種類
・\(N=0\,\)：密度一定のガス
・\(N=0.5\sim 1\,\)：中性子星
・\(N=\frac{3}{2}\,\)：非相対論的ガス球　白色矮星、ガス惑星　RGの核
・\(N=3\sim3.5\,\)：相対論的ガス　放射優勢の太陽などの恒星
・\(N=5\,\)：半径が無限大の星
・\(N=\infty\,\)：等温のガス球

［1］状態方程式
粒子がエネルギー\(\,\epsilon\,\)、化学ポテンシャル\(\,\mu\,\)に存在する確率は分布関数
\(\quad f(\epsilon)=\dfrac{1}{\exp[(\epsilon-\mu)/k_BT]\pm1}\quad\cdots\)(4)
符号\(\,\pm\,\)は+がフェルミ粒子、-がボース粒子に対応する。
高温あるいは低密度で粒子の存在期待値が小さい場合、量子力学の効果が無視出来て
\(\quad f(\epsilon)=e^{\mu/k_BT}\,e^{-\epsilon/k_BT}\,\)と近似出来る。
マックスウェル・ボルツマン分布が成り立つ条件は\(\,e^{\mu/k_BT}\ll1\,(\mu<0\,)\,\)であり、理想気体の状態方程式が数密度が以下を充たす場合である。\(A_m\,\)を粒子当たりの平均分子量とすると
\(\quad n=\dfrac{\rho}{A_mM_u}\ll\dfrac{g(2\pi mk_BT)^{3/2}}{h^3}\quad\cdots\,\)(5)

［2］光子ガス
光子はスピン1のボース粒子でボース-アインシュタイン統計に従う。エネルギーを\(\epsilon\)とすると、\(\epsilon=pc\,,\quad v=c\,\)で統計的重みは\(\,g=2\)、光子数は保存されないので\(\mu=0\)である。これより放射の圧力とエネルギー密度は
\(\quad U_r=3P_r=aT^4\quad\cdots\,\) (6)　ステファンボルツマンの式
放射の圧力(\(P_r\))の寄与の比を \(\,1-\beta\equiv P_r/P\,\)と表し、気体の比熱比を\(\,\gamma\,\)とすると、断熱指数は
\(\quad\Gamma_1=\dfrac{4}{3}+\dfrac{\beta(4-3\beta)(3\gamma-4)}{3\beta+36(\gamma-1)(1-\beta)}\quad\cdots\,\)(7)
と表せる。ガス圧が優勢の\(\,\beta=1\,\)の時は\(\,\Gamma_1=\gamma\,\)で、
放射優勢の極限\(\,\beta\to 0\,\)で\(\,\Gamma_1=4/3\,\)となる。

［3］電子縮退気体
化学ポテンシャル\(\,\mu\,\)と熱エネルギー\(\,k_BT\,\)との比として縮退パラメーター\(\,\Phi\,\)を定義する
\(\quad\Phi=\mu/k_BT\quad\cdots\,\)(8)
(5)式より縮退による理想気体からのズレは\(\,\Phi\simeq 0\,\)で現れるが、低温・高密度で\(\,\Phi\to\infty\,\)の極限で分布関数が階段関数で近似できる場合を完全縮退という。このときの化学ポテンシャルの最大値をフェルミエネルギー\(\,\epsilon_F\,\)、運動量をフェルミ運動量\(\,p_F\,\)という。
ここでフェルミモーメント\(\,p_F/m_ec=x_F\,\)とおくと、電子の数密度、圧力、運動エネルギー密度はそれぞれ
\(\quad n_e=\dfrac{8\pi}{h^3}\displaystyle\int_0^{p_F}\!p^2dp=\dfrac{8\pi}{3}\left(\dfrac{mc}{h}\right)^3x_F^3=Bx_F^3\quad\cdots\,\)(9)
\(\quad P_e=\dfrac{8\pi c}{3h^3}\displaystyle\int_0^{p_F}\!\dfrac{p^4dp}{\sqrt{p^2+m_e^2c^2}}=8A\!\displaystyle\int_0^{x_F}\!\dfrac{x^4\,dx}{\sqrt{1+x^2}}\cdots\,\)(10)
\(\quad U_e=\dfrac{8\pi}{h^3}\displaystyle\int_0^{p_F}(\sqrt{p^2c^2+m_e^2c^4}-m_ec^2)p^2dp\)
\(\qquad=\dfrac{8\pi m_e^4c^5}{h^3}\displaystyle\int_0^{x_F}\!(\sqrt{1+x^2}-1)x^2dx=Ag(x_F)\quad\cdots\,\)(11)
ここで、\(\,A\,,B\,\)はそれぞれ
\(\quad A\equiv\dfrac{\pi m_e^4c^5}{3h^3}=6.002\times10^{22}\)dyn cm\(^{-2}\)
\(\quad B\equiv\dfrac{8\pi}{3}\left(\dfrac{mc}{h}\right)^3=9.739\times10^5\)g cm\(^{-3}\)
また
\(\quad f(x)=8\displaystyle\int_0^x\!\dfrac{y^4\,dy}{\sqrt{1+y^2}}=\!x(2x^2-3)\sqrt{1+x^2}+3\sinh^{-1}x\)
\(\quad g(x)=24\displaystyle\int_0^x(\sqrt{1+y^2}-1)y^2\,dy\)
\(\qquad=3\{2x(x^2+1)^{3/2}-x\sqrt{x^2+1}-\sinh^{-1}x\}-8x^3\)
\(\qquad=8x^3(\sqrt{1+x^2}-1)-f(x)\)
のように定義されます。
\(\quad U_e=Ag(x_F)\,,\quad P_e=Af(x_F)\quad\)また、
\(\quad P=(\gamma-1)U\quad\) なので、電子ガスにおいて
\(\quad \gamma=\dfrac{f(x_F)}{g(x_F)}+1\quad\cdots\,\)(12)
と表すことが出来る。

図のように、\(p_F/m_ec=x_F\ll 1\,\)では、\(\gamma\to 5/3\,\)となり、
相対論的極限の \(x_F\gg 1\,\)では、\(\gamma\to 4/3\,\)となる。

ポリトロープとは、熱力学で使うときと天文で使うときで多少感じが違う。まずは、熱力から

● ポリトロープ変化

圧力\( P\)と比容積\( V\)が以下の関係で結ばれる変化をポリトロープ変化という（この変化は可逆変化）
\(\quad PV^n=C\quad n\): ポリトロープ指数
\(\qquad n=0\qquad P=C\quad\,\,\): 等圧変化
\(\qquad n=1\qquad PV=C\quad\): 等温変化
\(\qquad n=\gamma\qquad PV^{\gamma}=C\quad\): 断熱変化
\(\qquad n=\infty\quad\,\, V=C\qquad\): 等積変化

【1】気体の比熱

熱力学の第1法則は
\(\,U\,\)を内部エネルギー、\(\,Q\,\)を熱量、\(\,W\,\)を仕事とすると
\(\quad dU=d’W+d’Q\quad \cdots\)(1)
と表せる。ここで \( d’\,\)は準静的過程を表す。
\(\quad d’W=-PdV\quad \cdots\)(2)
\(\quad dU=-PdV+d’Q\quad \cdots\)(3)
体積が一定として、\(dT\,\)で割ると
\(\quad C_V=\dfrac{d’Q}{dT}=\left(\dfrac{\partial U}{\partial T}\right)_V\quad \cdots\)(4)
同様に、圧力が一定とすれば
\(\quad C_P=\left(\dfrac{\partial U}{\partial T}\right)_P+P\left(\dfrac{\partial V}{\partial T}\right)_P\quad \cdots\)(5)

また、(4)式を \(T\)で積分すると\(\,\,U=C_VT\quad \cdots\)(6)
1モルでは \(\,\,PV=\mathrm{R}T\quad \cdots\)(7)\(\,\,\)で
\(\quad P\left(\dfrac{\partial V}{\partial T}\right)_P=\mathrm{R}\quad \cdots\)(8)　なので、(5)式は
\(\quad C_P=C_V+\mathrm{R}\quad \cdots \)(9)

【2】断熱変化

\(dU\,\)の微分を考える
\(\quad dU=\left(\dfrac{\partial U}{\partial T}\right)_VdT+\left(\dfrac{\partial U}{\partial V}\right)_TdV\quad \cdots\)(10)
これを(3)式に代入すると
\(\quad d’Q=\left(\dfrac{\partial U}{\partial T}\right)_VdT+\left[P+\left(\dfrac{\partial U}{\partial V}\right)_T\right]dV\quad \cdots\)(11)
\(P=\)一定として(11)式を\(dT\,\)で割ると
\(\quad C_P=C_V+\left[P+\left(\dfrac{\partial U}{\partial V}\right)_T\right]\left(\dfrac{\partial V}{\partial T}\right)_P\quad \cdots\)(12)
断熱変化では \(d’Q=0\,\)なので(11)式は
\(\quad \left(\dfrac{\partial U}{\partial T}\right)_VdT+\left[P+\left(\dfrac{\partial U}{\partial V}\right)_T\right]dV=0\quad \cdots\)(13)
理想気体では\(\,(U\,\)は温度だけの関数なので \(\left(\dfrac{\partial U}{\partial V}\right)_T=0\)
(13)式に(6)式と(7)式を代入すると
\(\quad C_VdT+\dfrac{\mathrm{R}T}{V}dV=0\quad \cdots\)(14)
比熱比を\(\,\gamma=\dfrac{C_P}{C_V}\quad \cdots\)(15)\(\quad\)として(14)式を\(C_VT\,\)で割ると
\(\quad \dfrac{dT}{T}+(\gamma-1)\dfrac{dV}{V}=0\quad \cdots\)(16)\(\quad\)なので
\(\quad \ln T+(\gamma-1)\ln V= \text{一定} \cdots\)(17)\(\quad\)よって
\(\quad TV^{\gamma-1}= \text{一定} \cdots\)(18)
式(7)より\(\quad T=\dfrac{PV}{\mathrm{R}}\quad\cdots\)(19)を代入すると
\(\quad PV^{\gamma}=\text{一定}\quad \cdots\)(20)「ポアソンの法則」が得られる。
常温での\(\,C_V\,\)は並進運動と回転運動の自由度の\(\frac{\mathrm{R}}{2}\,\)倍となる。
単原子分子では並進運動の自由度が３なので、\(\,C_V=\frac{3}{2}\mathrm{R}\)。
２原子分子では\(\theta\,,\,\phi\,\)の回転が加わり、自由度が５になるので、\(\,C_V=\frac{5}{2}\mathrm{R}\)となる。
３原子分子(直線形)では、並進の自由度と回転の自由度は、２原子分子と同様だが振動の自由度が4となる。この振動は常温での比熱に寄与しないので、\(\,C_V\,\)は同じとなる。
よって、単原子分子は \(\gamma=\frac{5}{3}\,\)となり、２原子分子と直線形の３原子分子は \(\,\gamma=\frac{7}{5}\,\)となる。

［１］楕円の周の長さを求める

まず \(\,x=a\sin\theta\,,\,y=b\cos\theta\,\) で表される楕円の周の長さを求める．計算は第１象限の弧だけ求めて4倍することとする。

曲線の微少長さ\(\,ds\,\)はピタゴラスの定理から以下のようになる。

\(\quad ds=\sqrt{dy^2+dx^2}\quad\cdots\,\) (1)

弧の長さ(\(L\))は、これを全体に渡って積分することで求めることが出来る。

\(\quad x=a\sin\theta\,\) から\(\quad dx=a\cos\theta\,d\theta\,\)
\(\quad y=b\cos\theta\,\) から\(\quad dy=-b\sin\theta\,d\theta\,\)なので

\(\quad L=\displaystyle\int\sqrt{dy^2+dx^2}\,\)

\(\qquad =4\displaystyle\int_0^{\frac{\pi}{2}}\sqrt{b^2\sin^2\theta+a^2\cos^2\theta}\,d\theta\)

\(\qquad =4a\displaystyle\int_0^{\frac{\pi}{2}}\sqrt{1-\sin^2\theta+\frac{b^2}{a^2}\sin^2\theta}\,d\theta\)

ここで\(\,a\gt b\,\)として、\(\,q=\dfrac{a^2-b^2}{a^2}\,\)とおくと

\(\quad L=4a\displaystyle\int_0^{\frac{\pi}{2}}\sqrt{1-q\sin^2\theta}\,d\theta\,\)となる

［２］級数展開

\(\sqrt{1-q\sin^2\theta}\,\)はそのままでは積分できないので、級数展開する。

例えば、\(\,(1+x)^m\,\)をマクローリン展開すると

\(\quad(1+x)^m=1+m\dfrac{x}{1!}+m(m-1)\dfrac{x^2}{2!}+\,m(m-1)(m-2)\dfrac{x^3}{3!}\)

\(\qquad \cdots+m\cdots(m-n+1)\dfrac{x^n}{n!}+\)

ここで、例えば　\(m=-\dfrac{1}{2}\)を代入すると

\(\quad(1+x)^{-\frac{1}{2}}=1+\left(-\dfrac{1}{2}\right)x+\left(-\dfrac{1}{2}\right)\left(-\dfrac{3}{2}\right)\dfrac{x^2}{2!}\)

\(\qquad+\left(-\dfrac{1}{2}\right)\left(-\dfrac{3}{2}\right)\left(-\dfrac{5}{2}\right)\dfrac{x^3}{3!}+\cdots\)

\(\qquad =1-\dfrac{1}{2}\dfrac{x}{1!}+\dfrac{1\cdot 3}{2^2}\dfrac{x^2}{2!}-\dfrac{1\cdot 3\cdot 5}{2^3}\dfrac{x^3}{3!}+\)

◎ここでは、\((1-x)^m=1-m\dfrac{x}{1!}+m(m-1)\dfrac{x^2}{2!}-\)で

\(\quad m=\dfrac{1}{2}\,\qquad x=q\sin^2\theta\,\)とすると

\((1-q\sin^2\theta)^{\frac{1}{2}}=1-\dfrac{1}{2}q\sin\theta-\dfrac{1}{2^2\cdot 2!}q^2\sin^4\theta\)

\(\qquad-\dfrac{1\cdot 3}{2^3\cdot 3!}q^3\sin^6\theta\cdots-\dfrac{1\cdot 3\cdot 5\cdot(2n-3)}{2^n\cdot n!}q^n\sin^{2n}\theta\,\cdots\)

［３］二重階乗

二重階乗\(\,n!!\,\)は自然数\(\,n\,\)に対して一つ飛ばしに積を取る。

\(\quad(2n)!!=(2n)(2n-2)(2n-4)\cdots(4)(2)=2^nn!\)

\(\quad(2n+1)!!=(2n+1)(2n-1)\cdots(3)(1)=\dfrac{(2n+1)!}{(2n)!!}\,\)

なので、一般項の分母 \(\,2^nn!\,\)は\(\,(2n)!!\,\)となる。

分子は\(\quad 1\cdot 3\cdot 5\cdot(2n-3)=\dfrac{(2n-1)!!}{2n-1}\,\) と表せる

［４］\(\,\sin^n\theta\,\)の積分

[2]で出てきた\(\,\,\displaystyle\int_0^{\frac{\pi}{2}}\sin^{2n}\theta\,\,\)はそのままでは計算できない。

まず、部分積分を再確認。

\(\,\,\displaystyle\int f(x)g'(x)dx=f(x)g(x)-\displaystyle\int f'(x)g(x)dx\,\)である。

\(\,I_n=\displaystyle\int_0^{\frac{\pi}{2}}\sin^nx\,dx\,\,\)とすると、
\( I_0=\dfrac{\pi}{2}\,,\quad I_1=\left[-\cos x\right]_0^{\frac{\pi}{2}}=1\,\)となる

\(\,I_n=\displaystyle\int_0^{\frac{\pi}{2}}\sin^{n-1}x\sin x\,dx\,\)で

\(\,f(x)=\sin^{n-1}x\quad g(x)=-\cos x\,\)とおくと部分積分より

\(I_n=\left[-\sin^{n-1}x\cos x\right]_0^{\frac{\pi}{2}}+(n-1)\displaystyle\int_0^{\frac{\pi}{2}}\sin^{n-2}x\cos^2x\,dx\)

\(\,\,=0+(n-1)\displaystyle\int_0^{\frac{\pi}{2}}\sin^{n-2}x(1-\sin^2x)dx\)

\(\,\,=(n-1)\displaystyle\int_0^{\frac{\pi}{2}}\sin^{n-2}x\,dx-(n-1)\displaystyle\int_0^{\frac{\pi}{2}}\sin^nx\,dx\)

\(n\displaystyle\int_0^{\frac{\pi}{2}}\sin^nx\,dx=(n-1)\displaystyle\int_0^{\frac{\pi}{2}}\sin^{n-2}x\,dx\,\)なので

\(\displaystyle\int_0^{\frac{\pi}{2}}\sin^nx\,dx=\dfrac{n-1}{n}\displaystyle\int_0^{\frac{\pi}{2}}\sin^{n-2}x\,dx\quad\)となるので

\(I_n=\dfrac{n-1}{n}I_{n-2}\quad\)という漸化式が求められる。より

\(I_n=\displaystyle\int_0^{\frac{\pi}{2}}\sin^nx\,dx=\dfrac{(n-1))(n-3)\cdots}{n(n-2)(n-4)\cdots}\,\)

となり、(n,)が奇数の場合は

\(I_n=\dfrac{(n-1)(n-3)\cdots 4\cdot 2}{n(n-2)\cdots 5\cdot 3}I_1\,\)となり偶数では

\(I_n=\dfrac{(n-1)(n-3)\cdots 3\cdot 1}{n(n-2)\cdots 4\cdot 2}I_0\,\)となる。よって

\(\displaystyle\int_0^{\frac{\pi}{2}}\sin^{2n}\theta\,d\theta=\dfrac{(2n-1)!!}{(2n)!!}\dfrac{\pi}{2}\,\)となる。

[５］弧の長さの計算

これで、準備が終わったので \(\,L\,\)の長さを計算する。

\(L=4a\displaystyle\int_0^{\frac{\pi}{2}}\sqrt{1-q\sin^2\theta}\,d\theta\qquad q=(a^2-b^2)/a^2\)

\(\quad=4a\displaystyle\int_0^{\frac{\pi}{2}}\left(1-\displaystyle\sum_{n=1}^{\infty}\dfrac{(2n-1)!!}{(2n)!!}\dfrac{q^n}{2n-1}\sin^{2n}\theta\right)d\theta\)

\(\,=4a\displaystyle\int_0^{\frac{\pi}{2}}d\theta-\displaystyle\sum_{n=1}^{\infty}\dfrac{(2n-1)!!}{(2n)!!}\dfrac{q^n}{2n-1}4a\displaystyle\int_0^{\frac{\pi}{2}}\sin^{2n}\theta\,d\theta\)

\(\,=2a\pi\left(1-\displaystyle\sum_{n=1}^{\infty}\left[\dfrac{(2n-1)!!}{(2n)!!}\right]^2\dfrac{q^n}{2n-1}\right)\)

これで、計算完了。結構面倒くさい！！

［６］楕円積分

楕円の弧の長さの計算に使った積分は「第2標準形の完全楕円積分」というらしい。

\(E(m)=\displaystyle\int_0^{\frac{\pi}{2}}\sqrt{1-m\sin^2\theta}\,d\theta=\dfrac{\pi}{2}\left(1-\displaystyle\sum_{n=1}^{\infty}\left[\dfrac{(2n-1)!!}{(2n)!!}\right]^2\dfrac{m^n}{2n-1}\right)\)

振り子の計算で使った「第1標準形の完全楕円積分」は\(\,K(m)\,\)の形で表される積分で

\(K(m)=\displaystyle\int_0^{\frac{\pi}{2}}\dfrac{d\theta}{\sqrt{1-m\sin^2\theta}}\quad\)の様に書かれます。

楕円の弧の計算と同様に級数展開します。

\((1-m\sin^2\theta)^{-\frac{1}{2}}=1+\dfrac{1}{2}m\sin^2\theta+\dfrac{3}{8}m^2\sin^4\theta\)

\(\qquad+\dfrac{15}{48}m^3\sin^6\theta+\cdots=\displaystyle\sum_{n=1}^{\infty}\dfrac{(2n-1)!!}{(2n)!!}m^n\sin^{2n}\theta\)

この級数も \(\,m\lt 1\,\)の範囲で一様収束するので、項別積分が可能です。

\(\displaystyle\int_0^{\frac{\pi}{2}}((1-m\sin^2\theta)^{-\frac{1}{2}}d\theta=\displaystyle\int_0^{\frac{\pi}{2}}\displaystyle\sum_{n=1}^{\infty}\dfrac{(2n-1)!!}{(2n)!!}m^n\sin^{2n}\theta\)

\(\qquad=\displaystyle\sum_{n=1}^{\infty}\dfrac{(2n-1)!!}{(2n)!!}m^n\displaystyle\int_0^{\frac{\pi}{2}}\sin^{2n}\theta\,d\theta\)

\(\qquad=\dfrac{\pi}{2}\left(\displaystyle\sum_{n=1}^{\infty}\left[\dfrac{(2n-1)!!}{(2n)!!}\right]^2\,m^n\right)\)

エネルギー保存の法則から振り子の周期 \(T\) の式を求める。
振り子の支点より水平方向に \(\ x\ \)軸をとり、鉛直方向に \(\ y\ \)軸をとると、エネルギー保存則から、次式が成立する

\( \dfrac{1}{2}mv^2=mg\Delta h\quad\) (21)

ここで、\(\ v\ \)は重りが円周方向に移動する速度で、\(\Delta h\ \)は重りの落下距離である。ここで

\(\Delta h=l\cos\theta-l\cos\theta_0=l\ (\cos\theta-\cos\theta_0)\quad\) (22)
なので

\( \dfrac{1}{2}mv^2=mgl(\cos\theta-\cos\theta_0)\)
\( \therefore \ v=\sqrt{2gl(\cos\theta-\cos\theta_0)}\quad\) (23)

さらに

\( v=\dfrac{ds}{dt}=l\dfrac{d\theta}{dt}\quad\) (24)

なので、以下の式が得られる。

\( l\dfrac{d\theta}{dt}=\sqrt{2gl(\cos\theta-\cos\theta_0)}\quad\)(25)

\( dt=\dfrac{1}{\sqrt{2}\omega}\dfrac{d\theta}{\sqrt{\cos\theta-\cos\theta_0}}\quad\) (26)

\(\quad\because\,\,\omega:=\sqrt{\dfrac{g}{l}}\)

この(26)式を積分すると周期を求めることが出来る。積分範囲を4分の1周に相当する \(\ t=T/4\ \)とすると、右辺の積分範囲は\(\ \theta=0\sim\theta_0\ \)となるので、

\( \displaystyle\int_0^{\frac{T}{4}}dt=\dfrac{1}{\sqrt{2}\omega}\displaystyle\int_0^{\theta_0}\dfrac{d\theta}{\sqrt{\cos\theta-\cos\theta_0}}\quad\) (27)

\(\therefore\,\, T=4\dfrac{1}{\sqrt{2}\omega}\displaystyle\int_0^{\theta_0}\dfrac{d\theta}{\sqrt{\cos\theta-\cos\theta_0}}\quad\) (28)

この積分を行うために、半角公式を使って以下の変換を行う

\( \cos\theta=1-2\sin^2\frac{\theta}{2}\quad\cos\theta_0=1-2\sin^2\frac{\theta_0}{2}\)

ここで \(\,\,k:=\sin\frac{\theta_0}{2}\quad \sin\phi:=\sin\frac{\theta}{2}/k\,\,\)とすると

\( \cos\theta-\cos\theta_0=1-2\sin^2\frac{\theta}{2}-\left(1-2\sin^2\frac{\theta_0}{2}\right)\)

\(\quad=2\left(\sin^2\frac{\theta_0}{2}-\sin^2\frac{\theta}{2}\right)=2k^2\cos^2\phi\,\,\)となる

\(k\sin\phi=\sin\frac{\theta}{2}\,\,\)を両辺微分して\(\,\,k\cos\phi\,d\phi=\frac{1}{2}\cos\frac{\theta}{2}\,d\theta\)

\(d\theta=\dfrac{2k\cos\phi\,d\phi}{\cos\frac{\theta}{2}}=\dfrac{2k\cos\phi\,d\phi}{\sqrt{1-\sin^2\frac{\theta}{2}}}=\dfrac{2k\cos\phi\,d\phi}{\sqrt{1-k^2\sin^2\phi}}\)

なので

\(\dfrac{d\theta}{\sqrt{\cos\theta-\cos\theta_0}}=\dfrac{1}{\sqrt{2}k\cos\phi}\dfrac{2k\cos\phi\,d\phi}{\sqrt{1-k^2\sin^2\phi}}\)

\(\quad=\dfrac{\sqrt{2}\,d\phi}{\sqrt{1-k^2\sin^2\phi}}\,\,\)となり、第１種楕円積分の形になる。

［１］運動方程式

振り子には、図のような力が働く。ここで 重力は\( \ W=mg\ \)で糸の張力は\(\ \ T=W\cos\theta\ \) である。

\(W\ \)と\(\ T\ \)の合力\(\ f\ \)の向きは図のように半径\(\ l\ \)の円周の接線方向となり、その大きさは

\(f=W\sin\theta=mg\sin\theta\)　(1)

である。この合力\(\ f\ \)により、重りは円弧を描くように往復運動をする。ここで、運動を円軌道として考えるために、振り子の釣り合いの位置 O から、円弧に沿って\(\ s\ \)軸をとる。この合力\(\ f\ \)は重りを\(\ s\ \)軸の負の方向に運動させる力となるため、重りの運動方程式は

\( m\dfrac{d^2s}{dt^2}=-f=-mg\sin\theta\ \)　(2)

すなわち

\( \dfrac{d^2s}{dt^2}=-g\sin\theta\ \)　(3)

と表すことが出来る。　円弧の長さは \(\ s=l\ \theta\ \)で表せられ、２回微分は

\(\ \dfrac{d^2s}{dt^2}=l\dfrac{d^2\theta}{dt^2}\ \)　(4)

となるので、運動方程式は重りの角度\(\ \theta\ \)に対して

\( \dfrac{d^2\theta}{dt^2}=-\dfrac{g}{l}\sin\theta\ \)　(5)

と表すことが出来る。

［２］運動方程式の近似解

　(5)式の運動方程式は、三角関数を含んでいるため、このままでは階を求めることが難しい。但し、円弧の角度\(\ \theta\ \)が微少の場合、三角形の高さ\(\ x=\sin\theta\ \)は円弧の長さ\(\ s\ \)にほぼ等しくなり
\(\quad x\ \approx\ s\ \)とおけるため

\( \sin\theta=\dfrac{x}{l}\ \approx\ \dfrac{s}{l}=\dfrac{l\theta}{l}=\theta\ \)　(6)

となり、運動方程式は

\( \dfrac{d^2\theta}{dt^2}=-\omega^2\theta\qquad \omega:=\sqrt{\dfrac{g}{l}}\quad\) (7)

と表すことが出来る。これは単振動で良く近似出来る。一般解は

\( \theta=A\sin\omega t+B\cos\omega t \)　(8)

ここで、 \( t=0\ \)で\(\ \theta=\theta_0\ \)として、(8)式に代入すると

\( \theta_0=A\sin\omega\times 0+B\cos\omega\times 0=B\quad\therefore B=\theta_0\)

つぎに、(8)式を \(t\) で微分して、\(\ t=0\ \)で\(\ \dfrac{d\theta}{dt}=0\ \) とすると

\( \dfrac{d\theta}{dt}=A\omega\cos\omega t-B\omega\sin\omega t\)

\(\qquad =A\omega\cos\omega\times 0-B\omega\sin\omega\times 0=A\omega=0 \)

よって、\(\ A=0\ \)となり、(8)式は\(\quad\theta=\theta_0\cos\omega t\ \) (9)

と決定される。このように解が周期関数で表される振動を単振動といい、単振動で近似出来る振り子を単振り子といいます。

\(\ \omega=\sqrt{\cfrac{g}{l}}\ \) は角振動数 [rad/s] と呼ばれる。

角振動数\(\ \omega\ \)を 2\(\pi\) で割った値は、重りが1秒間に往復する回数を表し、振動数 \(\ f\ \)[Hz}と呼ばれています。

\( f=\dfrac{\omega}{2\pi}=\dfrac{1}{2\pi}\sqrt{\dfrac{g}{l}}\ \,\,\) Hz\(\quad\)(10)

重りが一往復するのに必要な時間は、周期 \(T\) [s] といい、振動数の逆数になる。

\(T=\dfrac{1}{f}=\dfrac{2\pi}{\omega}=2\pi\sqrt{\dfrac{l}{g}}\quad\)(11)

TeX test

$$ x^2=c$$
$$\nabla=\dfrac{\partial}{\partial x} $$