2026-07-30

ライブ前に考えたい騒音と音の減衰：物理学の視点から

物理学

この記事の目的
粘性( 古典力学 ) による減衰
内部振動による緩和
人間の聴覚特性
音の減衰のグラフ
まとめ
参考文献

この記事の目的

今年は近年では珍しく梅雨らしい梅雨を迎え、エルニーニョ現象もあり涼しい夏が久々に来るかもという私の予想（希望ともいう）に反して、梅雨明けから猛暑になり今年も暑い夏を迎えそうな感がある。夏になると盛んに行われるのが野外音楽ライブだろう。私自身は音楽に明るくないので、どこで誰がどういうライブをするのかというのは良くわからない。そんな私でも去年の7月末にとあるミュージシャンのライブが騒音で話題になったのを耳にすることがあった。このときに主に重低音に言及している人が多かったのと、物理的に考えると重低音の方が遠くに届くという記事をみかけたのをふと思い出したので、それを数式を用いて少し物理学的に考えたい。

最後の章に行く途中でかなり数式が出てくるので、難しい数式を見たくない人は「音の減衰のグラフ」の部分だけを見てもらえればと思う。

粘性( 古典力学 ) による減衰

この章ではまず、音波が空気の粘性によってどのように弱まるかを見る。この部分だけを考えると、減衰はおおよそ周波数の二乗に比例するため、高い音ほど空気中で失われやすい。
物理として考えるにあたって、音というのは気体の微小振動であるので、まずは気体や液体（合わせて流体と言う)を記述するNavier-Stokes方程式 [1] から減衰のある場合に音波を記述する波動方程式はどんな形になるのかを考える。

$\displaystyle \frac{\partial }{\partial t}\bf{v} + (\bf{v}\cdot\nabla)\bf{v} = -\frac{1}{\rho}\nabla p + \frac{\mu}{\rho}\nabla^2 \bf{v}+\frac{\chi+\frac{1}{3}\mu}{\rho}\nabla(\nabla\cdot \bf{v})$

これが Navier-Stokes方程式である。 $\bf{v}$ は流体の速度、 $\rho$ は密度、 $p$ は圧力を表している。
また、流体は次の連続の式を満たす。

$\displaystyle \frac{\partial \rho}{\partial t} + \nabla\cdot (\rho\bf{v})= 0$

音は、速度が0で静止している流体に微小な動きが加わったときを考える。静止した状態では速度は $\bf{0}$ , 密度は $\rho_0$ , 圧力は $p_0$ である。この静止した状態からの摂動を考える。すると速度は $\bf{u}$ , 密度は $\rho_0 + \rho'$ , 圧力は $p_0 + p'$ であり、 $|{\bf u}| \ll 1$ , $\rho' \ll \rho_0$ , $p' \ll p_0$ を満たすとする。
Navier-Stokes方程式と連続の式に $\bf{u}$ , $\rho_0 + \rho'$ , $p_0 + p'$ を代入して、微小量の一次までとって線形化すると以下のようになる。

$\frac{\partial }{\partial t}\bf{u} = -\frac{1}{\rho_0}\nabla p' + \frac{\mu}{\rho_0}\nabla^2 \bf{u}+\frac{\chi+\frac{1}{3}\mu}{\rho_0}\nabla(\nabla\cdot \bf{u})$

$\frac{\partial \rho'}{\partial t} + \rho_0\nabla\cdot (\bf{u})= 0$

ここで、波動方程式を導出するために $\rho'$ の時間に関する二階微分を考えると、

$\frac{\partial^2 \rho'}{\partial t^2} = -\rho_0\nabla\cdot \left(\frac{\partial \bf{u}}{\partial t}\right) \\ =-\rho_0\nabla\cdot \left\{ -\frac{1}{\rho_0}\nabla p' + \frac{\mu}{\rho_0}\nabla^2 \bf{u}+\frac{\chi+\frac{1}{3}\mu}{\rho_0}\nabla(\nabla\cdot \bf{u}) \right\} \\ =\nabla^2 p'-\left(\chi+\frac{4}{3}\mu\right)\nabla^2(\nabla\cdot\bf{u}) \\ =\nabla^2 p'+\frac{\left(\chi+\frac{4}{3}\mu\right)}{\rho_0}\nabla^2\left(\frac{\partial \rho'}{\partial t}\right)$

ここで、状態方程式 $p' = \left(\frac{\partial p}{\partial \rho}\right)_s\rho'$ を用いると、最終的に次の波動方程式を得る。

$\left(\frac{\partial^2 }{\partial t^2} - c_s^2 \nabla^2\right) \rho'=\frac{\left(\chi+\frac{4}{3}\mu\right)}{\rho_0}\nabla^2\left(\frac{\partial \rho'}{\partial t}\right)$

ここで、 $c_s^2 = \left(\frac{\partial p}{\partial \rho}\right)_s$ とした。この $c_s$ が音波の速度になる。この微分方程式において粘性係数が0であれば、右辺が0になり一般的な波動方程式と一致する。ここで、密度 $\rho'$ が振動すると仮定して ( $\rho' \propto \exp({\rm i}{\bf k}\cdot{\bf x}-{\rm i}\omega t)$ )、微分方程式に代入することで、次のような波数と角周波数の間の分散関係を得る。

$-\omega^2+c_s^2k^2 = {\rm i} \frac{\chi + \frac{4}{3}\mu}{\rho_0} k^2\omega$

これを平方完成して $\omega$ について解くと、

$\omega = \pm \sqrt{c_s^2k^2-\left(\frac{\chi+\frac{4}{3}\mu}{2\rho_0}k^2\right)^2}-{\rm i}\frac{\chi+\frac{4}{3}\mu}{2\rho_0}k^2$

となり、さらに弱粘性極限を考えると、

$\omega \simeq \pm c_sk - {\rm i}\frac{\chi+\frac{4}{3}\mu}{2\rho_0}k^2$

となる。

以下では音波の形として球形波を考える。これは、点音源から音波が均等に放出されているという単純化をしている。(そのため指向性がある場合は必ずしも以下の議論は成立しない場合もある。）

弱粘性極限では、微分方程式の解の形はほぼ変わらないと考えて良いので、

$p' \simeq p'_0 \frac{\exp\left({\rm i}kr-{\rm i}\omega t\right)}{r} \\ \simeq p'_0 \frac{\exp\left({\rm i}kr \mp {\rm i}kc_s t \right)\exp\left(-\frac{\chi+\frac{4}{3}\mu}{2\rho_0}k^2t\right)}{r}$

$\exp\left(-\frac{\chi+\frac{4}{3}\mu}{2\rho_0}k^2t\right)$ が粘性による音波の減衰を表す。ここで、時間は音の進む距離を音速で割ったものなので $t = r/c_s$ となることを使うと減衰は次のように書ける。

$p' \simeq p'_0 \frac{\exp\left({\rm i}kr \mp {\rm i}kc_s t \right)\exp\left(-\frac{\chi+\frac{4}{3}\mu}{2\rho_0c_s}k^2r\right)}{r} \\ = p'_0 \frac{\exp\left({\rm i}kr \mp {\rm i}kc_s t \right)\exp\left(-\alpha_{\rm fluid}r\right)}{r} \\ \alpha_{\rm fluid} = \frac{\chi+\frac{4}{3}\mu}{2\rho_0c_s^3}\omega^2$

ここでは $k = \omega/c_s$ となることを用いた。

以上では計算を簡略化するため粘性のみを考慮したが、実際には熱伝導によるエネルギー損失も同程度に寄与する。熱伝導のエネルギー損失は、上記の Navier-Stokes 方程式と連続の式に加えてエネルギー輸送方程式を考えてこれら3つの式を連立させることで導くことができる。これらを合わせた古典的減衰係数は、実効的な粘性係数を用いて

$\alpha_{\rm cl} = \frac{\omega^2}{2\rho_0c_s^3}\left[\left(\chi+\frac{4}{3}\mu\right)+\kappa\left(\frac{1}{c_v} - \frac{1}{c_p} \right)\right] = \frac{\omega^2}{2\rho_0c_s^3}\left[\left(\chi+\frac{4}{3}\mu\right)+\frac{\kappa}{c_p}\left(\gamma - 1 \right)\right] \\ = \frac{2\pi^2 f^2}{\rho_0c_s^3}\left[\left(\chi+\frac{4}{3}\mu\right)+\frac{\kappa}{c_p}\left(\gamma - 1 \right)\right]$

とあらわすことができる[2]。ここで、最後の式変形では $\omega = 2\pi f$ を用いた。

この式の前半は流体の粘性による寄与であり、式の後半が熱伝導の効果である。また、 $\kappa$ , $c_v$ , $c_p$ , $\gamma$ はそれぞれ熱伝導率, 定積比熱, 定圧比熱, 比熱比である。これらは古典力学によって導かれるもので、粒子同士が弾性衝突してエネルギーが分散されている過程によるものである。
ここで、 $\alpha_{\rm cl}$ の大きさがどのくらいか具体的に数値を見ておく。300K での空気の粘性係数 $\mu$ は $1.845 \times 10^{-5}$ $\rm Pa \ s$ であり [3]、空気の体積粘性係数 $\chi$ は以下の式で代表値として $1.794 \times 10^{-5}$ $\rm Pa \ s$ とする [4]。

$\chi(\text{air}) = \chi_0(\text{air}) + cT \\ \chi_{0}(\text{air}) = (-5.316 \pm 0.245) \times 10^{-5} {\rm Pa \ s} \\ c = (0.0237 \pm 0.0007) \times 10^{-5} {\rm Pa \ s \ K^{-1}}$

この組み合わせで $\chi+\frac{4}{3}\mu$ は $4.254 \times 10^{-5}$ $\rm Pa \ s$ 程度の大きさになる。空気の密度は $1.177$ $\rm kg \ m^{-3}$ 程度とする [5]。音速は $c_s \simeq 348 \ {\rm m/s}$ とし、比熱比 $\gamma = 1.402$ 、定圧比熱は $c_p = 1006 \ {\rm J/ kg \ K}$ を代表値として選び、熱伝導係数は $\kappa = 0.0262 \ {\rm W / m \ K}$ を代表値とする。これらの値を $\alpha_{\rm cl}$ の式に代入して計算すると、周波数 $f$ に対して $\alpha_{\rm cl} \simeq 2.109 \times 10^{-11} \cdot f^2 \ {\rm Np / m}$ となる。

内部振動による緩和

ここまでは古典的な流体力学を使った音波の減衰を見てきた。この減衰は微視的には空気の分子同士が弾性衝突していくことで、エネルギーが分散していく効果に由来している。
分子の内部にバネのような振動の自由度があると話が少し複雑になる。音波では空気が周期的に圧縮・膨張し、それに伴って温度も周期的に変化する。しかし、分子内部の振動エネルギーが温度変化に瞬時に追随できないと、エネルギーの出し入れに遅れが生じる。この遅れが、音波のエネルギーを分子内部自由度へ逃がすため、追加の吸収として現れる。
ここでは、分子は二原子が角振動数が $\omega_0$ のバネで結合している模型を考える。したがって、この振動は量子化されていて、分子同士の衝突によって振動エネルギーが変化する効果について考える。この調和振動子は外部の熱浴と ${\mathcal H}_{\rm int} = k\hat{x}\hat{b}$ というハミルトニアンで結合しているとする。ここで、 $\hat{x}$ はバネの変位の演算子、 $\hat{b}$ は熱浴の演算子である。
このとき、系が $n$ 準位にある確率を $P_n$ とすると、その遷移を表すマスター方程式は次のように書くことができる [6]。

$\frac{\rm d}{{\rm d} t} P_n = -k_{n, n+1}P_n + k_{n+1, n}P_{n+1} -k_{n, n-1}P_n + k_{n-1, n}P_{n-1}$

ここで、 $k_{i, j}$ は $i$ 番目の準位から $j$ 番目の準位への遷移を表す係数で、フェルミの黄金律から $|\langle i | k\hat{x}\hat{b}|j\rangle|^2$ に比例し、生成消滅演算子を用いると $\hat{x} = \sqrt{\frac{\hbar}{2m_0\omega_0}}(\hat{a} + \hat{a}^\dagger)$ より $k_{n, n-1} = nk_{1, 0}$ , $k_{n + 1, n} = (n + 1) k_{1, 0}$ を得る。
さらに、平衡系を考えると $n$ 番目の準位から $n - 1$ 番目の準位への遷移と $n - 1$ 番目の準位から $n$ 番目の準位への遷移とが等しくなるとする。つまり $k_{n, n-1}P_n = k_{n-1, n}P_{n-1}$ が平衡系で成り立つ（詳細釣り合い条件という）と仮定する。ここで、 $P_n^{\rm eq}$ を平衡系での確率とするとこれはボルツマン分布なので $P_n^{\rm eq}\propto \exp\left(-\beta \hbar\omega_0(n + \frac{1}{2}) \right)$ となる。これを使うと次を得る。

$k_{n-1, n} = k_{n, n-1} \frac{P_n^{\rm eq}}{P_{n-1}^{\rm eq}} = k_{n, n-1} e^{-\beta\hbar\omega_0}$

このとき、基底状態エネルギーの基準位置を取り直してエネルギーの期待値は

$E_v = \hbar\omega_0 \sum_{n = 0}^\infty nP_n$

となり、このエネルギーを時間について微分しマスター方程式を用いると最終的に次を得る。

$\frac{\rm d}{{\rm d}t}E_v =\frac{1}{\tau}(E_v^{\rm eq}-E_v) \\\\E_v^{\rm eq} = \frac{\hbar\omega_0}{e^{\beta\hbar\omega_0} -1} \quad , \quad \frac{1}{\tau} = k_{1, 0}(1-e^{-\beta\hbar\omega_0})$

つまりバネのエネルギーは熱浴との相互作用によって逆温度 $\beta$ での振動エネルギーへと緩和していく。
音波では温度が $T = T_0 + T_1 e^{-{\rm i}\omega t}$ となるので ( $|T_1|\ll T_0$ ) 、 $E_v^{\rm eq}$ を $T_1$ の一次の項まで展開し、 $E_v$ を $E_v^{\rm eq}(T_0)$ からのずれで書くと

$E_v^{\rm eq}(T) = E_v^{\rm eq}(T_0) + \frac{\partial E_v^{\rm eq}}{\partial T}T_1 e^{-{\rm i}\omega t} \\E_v = E_v^{\rm eq}(T_0) + \delta E e^{-{\rm i}\omega t}$

となり、これをエネルギーの時間発展の式に代入して整理すると次のようになる。

$(-{\rm i}\omega + \tau^{-1})\delta E = \tau^{-1} \frac{\partial E_v^{\rm eq}}{\partial T}T_1 \\\delta E = \frac{1}{(1 -{\rm i}\omega\tau )} \frac{\partial E_v^{\rm eq}}{\partial T}T_1 = \frac{C_v^{vib}}{(1 -{\rm i}\omega\tau )} T_1$

ここで、 $C_v^{vib} = \frac{\partial E_v^{\rm eq}}{\partial T}$ は振動の比熱である。よって、実効的な振動の比熱は複素数になることが分かる。これに、分子の並進と回転の比熱 $C_v^{rot, tr}$ を加えたものが全体の比熱となる。

$C_v(\omega) = C_v^{rot, tr} + \frac{C_v^{vib}}{1-{\rm i}\omega\tau}$

音は断熱的な振動なので熱力学第一法則から、 $dU = C_v(\omega)dT = (p/\rho^2)d\rho$ が成り立つ。
ここで、理想気体の状態方程式 $p = \rho RT$ において平衡状態 $(p_0, \rho_0, T_0)$ とそこからの摂動 $(p_0 + p', \rho_0 + \rho', T_0 + T')$ を考えると次のように線形化できる。

$\frac{p'}{p_0} = \frac{T'}{T_0} + \frac{\rho'}{\rho_0}$

熱力学第一法則より断熱過程において $C_v(\omega) T' = (p_0/\rho_0^2)\rho'$ が成り立つので、これを代入して整理すると

$p' = RT_0\left(1 + \frac{R}{C_v(\omega)}\right)\rho'$

理想気体における音速は $c_s^2 = (1+R/C_v)RT_0$ で与えられるので、この式から音速は $c^2(\omega) =RT_0\left(1 + \frac{R}{C_v(\omega)}\right)$ であることが分かる。ここで、 $C_v^{vib} \ll C_v^{rot, tr}$ であることを用いると、 $C_v(\omega)$ の逆数は

$\frac{1}{C_v(\omega)} = \frac{1}{C_v^{rot, tr}}\left(1+\frac{C_v^{vib}/C_v^{rot, tr}}{1-{\rm i}\omega\tau}\right)^{-1} \\ \simeq \frac{1}{C_v^{rot, tr}}\left(1-\frac{C_v^{vib}/C_v^{rot, tr}}{1-{\rm i}\omega\tau}\right)$

となり、これを用いて音速を書き直すと次のようになる。

$c^2_s(\omega) = c_\infty^2\left(1-\frac{m}{1-{\rm i}\omega\tau}\right) \\c_s(\omega) = c_\infty\left(1-\frac{m}{2}\frac{1}{1-{\rm i}\omega\tau}\right)$

ただし、 $c_\infty$ , $c_0$ はそれぞれ $\omega = \infty$ , $\omega = 0$ での音速で、 m は無次元定数で次のように定義される

$c_\infty^2 = RT_0\left(1 + \frac{R}{C_v^{rot, tr}}\right) \quad , \quad c_0^2 = RT_0\left(1 + \frac{R}{C_v^{rot, tr} + C_v^{vib}}\right) \quad , \quad m = \frac{c_\infty^2 - c_0^2}{c_\infty^2}$

$C_v^{vib} \ll C_v^{rot, tr}$ より $m \ll 1$ である。波数は $k = \omega / c_s(\omega)$ で計算されるので

$k = \frac{\omega}{c_s(\omega)} = \frac{\omega}{c_\infty}\left(1 + \frac{m}{2}\frac{1}{1-{\rm i}\omega\tau}\right) = \frac{\omega}{c_\infty}\left(1 + \frac{m}{2}\frac{1+{\rm i}\omega\tau}{1+(\omega\tau)^2}\right)$

よって緩和による減衰項は $\alpha_{vib} = \frac{m}{2c_\infty\tau}\frac{(\omega\tau)^2}{1 + (\omega\tau)^2}$ となる。
空気の分子は $N_2$ と $O_2$ であると近似すると、緩和の効果は $N_2$ と $O_2$ とのそれぞれの緩和による減衰の和になる。
角周波数を使っているが、 $\omega = 2\pi f$ から周波数 $f$ を使い $\tau = (2\pi f_r)^{-1}$ と置いて書き直すと

$\alpha_{\rm relaxation} = \frac{m}{2c_\infty\tau}\frac{(\omega\tau)^2}{1 + (\omega\tau)^2} \\ = \frac{2\pi f_r m}{2c_\infty}\frac{(f/f_r)^2}{1 + (f/f_r)^2} \\ =\frac{\pi m}{c_\infty}\frac{f^2f_r}{f_r^2 + f^2}$

という形になる。これは、 $f \ll f_r$ の低周波数極限では古典減衰のように $\alpha_{\rm relaxation} \propto f^2$ となり、 $f_r \ll f$ の高周波数側では $\alpha_{\rm relaxation} \propto f_r$ もしくは $\alpha_{\rm relaxation} \propto \tau^{-1}$ となることがわかる。したがって、緩和吸収は低周波側では $f^2$ に比例して増加し、緩和周波数を越えると単一緩和モードとしては飽和する。一方で、古典減衰は $f^2$ に比例して増え続けるため、全体としては高周波ほど空気吸収を受けやすい傾向になる。
最終的には $\alpha_{\rm tot} = \alpha_{\rm cl} + \alpha_{O_2} +\alpha_{N_2}$ と近似でき、その表式は次のように与えられることが知られている[7]。

$\alpha_{\rm tot} = f^2 \left\{ 1.84\times10^{-11}\left(\frac{p_m}{p_0}\right)^{-1}\left(\frac{T}{T_0}\right)^{\frac{1}{2}} + \left(\frac{T}{T_0}\right)^{-\frac{5}{2}}\times\left[0.01275e^{-2239/T}\left(\frac{f_{rO}}{f_{rO}^2+f^2} \right) +0.1068e^{-3352/T}\left(\frac{f_{rN}}{f_{rN}^2+f^2} \right)\right] \right\} \\f_{rO} = \frac{p_m}{p_0}\left\{24 + \left[\frac{(4.04\times10^4h)(0.02 + h)}{0.391 + h}\right]\right\} \\f_{rN} =\frac{p_m}{p_0}\left(\frac{T}{T_0}\right)^{-\frac12}\left(9 + 280h\exp\left\{-4.170\left[\left(\frac{T}{T_0}\right)^{-\frac{1}{3}} - 1 \right]\right\}\right)$

ここで $h$ は水蒸気のモル濃度を表しており、相対湿度 $RH$ ( ％ )と飽和水蒸気圧 $p_{\rm sat}$ に対して $h = RH \frac{p_{\rm sat}}{p_0}$ で与えられ、 $\frac{p_{\rm sat}}{p_0} = 10^C$ , $C = -6.8346 \times \left(\frac{273.16}{T}\right)^{1.261} + 4.6151$ となっている。また、 $p_m$ は気体の圧力である。
この式では $\alpha_{\rm tot}$ の前半は $\alpha_{\rm cl}$ を表しており、後半は酸素と窒素の緩和による減衰を表している。ここで先ほどの式と違うのは熱伝導や粘性係数などは気体の圧力や温度の関数であることを反映している。さきほど、300K で $\alpha_{cl} \simeq 2.109 \times 10^{-11} \cdot f^2 \ {\rm Np / m}$ と計算したが、 $1.84\times10^{-11}$ とオーダーが同じで近い値になっている。
また、緩和による減衰の温度依存性の形はLandau-Teller 理論などを使って与えられる。

人間の聴覚特性

音の強度を表す指標としてデシベルというものがある。音圧レベルは、音圧の実効値を基準音圧と比較して対数表示した量で、音波の圧力を用いた定義は $dB = 20\log_{10}(p_{\rm rms}/p_{\rm ref})$ である。ここで、 $p_{\rm ref}$ は基準となる音波の圧力である。空気中の音圧レベルでは、基準音圧として $p_{\rm ref} = 20 \mu {\rm Pa}$ が用いられる。これは 1 kHz 付近の最小可聴閾に近い値（つまり人が聞こえる最小値）である。
では実際にどのようになるかを見ていきたい。まず、 $\alpha$ で減衰する球面波の圧力は次のようになる。

$p = p_{\rm amp}\frac{e^{{\rm i}(kr -\omega t)}e^{-\alpha r}}{r}$

ここで、 $p_{\rm amp}$ は音圧の振幅である。デシベルの定義に、この $\alpha$ で減衰する球面波の音波の圧力の式を入れると

$dB = 20 \log_{10}\left( \ \left| p/p_{\rm ref}\right| \ \right) \\ = 20 \log_{10}(p_0/p_{{\rm ref}, 0}) -20 \log_{10}(r/r_{\rm ref}) -20\log_{10}(e)\alpha (r - r_{\rm ref})$

となる。この式の第二項の $-20\log_{10}(r/r_{\rm ref})$ は球面波の形状に由来する幾何学的な減衰をあらわしている。 $-20\log_{10}(e)\alpha r$ は空気の粘性・熱伝導・緩和吸収による減衰を表している。
ちなみに文献によっては自然対数 ln による音圧レベル Np と $\log_{10}$ による音圧レベル dB の変換の関係で $\alpha$ に最初から $20\log_{10}(e) \simeq 8.686$ をかけているものもある。例えば上記の $\alpha_{\rm tot}$ は $\alpha_{\rm tot} = 8.686 \times f^2 \cdots$ のように修正されることになる。

さて、騒音問題としてとらえる場合には、人間の生物学的な特性が必要になる。人間の聴覚は周波数によって感度の差があり、これをグラフとして表したものを等ラウドネス曲線という。そして、この等ラウドネス曲線において人間の耳の感度は数kHz にピークがあり、重低音および10kHz以上の高音側では聴覚の感度が低いことがわかる。そのため、高周波数と低周波数で同じ音圧レベル（デシベル）でも聞こえやすさが違うのでそこを考慮する必要がある。

等ラウドネス曲線の、とくに低い音圧レベルで聞こえにくくなるという人間の聴覚の周波数依存性を反映するために周波数重み付けして表したものをA特性音圧レベルという。
そして周波数重み付けをほぼ行わないものは Z特性といい、A特性より低周波をあまり減点しないが完全なフラットではないものを C特性という。
その補正項は次で与えられる[8, 9]。ここで $f$ は振動数である。

$A(f) = 20\log_{10}(R_A(f)) - 20\log_{10}(R_A(1000)) \\R_A(f) = \left[\frac{f_4^2 \ f^4}{(f^2 + f_1^2)\sqrt{(f^2 + f_2^2)(f^2 + f_3^2)}(f^2 + f_4^2)}\right] \\\\f_1 = 20.6 \ , \ f_2 = 107.7 \ , \ f_3 = 737.9 \ , \ f_4 = 12194 \\ -20\log_{10}R_A(1000) \simeq 2.0$

この項を補正として入れたものをdBA と書いたりする。この $A(f)$ をプロットした図は次のようになっている。

この曲線をみると、1 kHz 付近ではほぼ 0 dB に正規化され、2〜5 kHz 付近でわずかに正の補正を持ち、低音側では著しく下がることが分かる。

なお、A特性はあくまで人間の聴覚感度を反映した重み付けであり、低周波音の物理的なエネルギーそのものが減るわけではない。低周波音は人の耳には聞こえにくくても、建物の窓や壁を振動させたり、体感的な圧迫感として問題になることがあるため、A特性で「小さく評価される」ことと「実際に届いていない」ことは同義ではない点を先に断っておきたい。

音の減衰のグラフ

さて、もともとはライブの音の減衰の話だったので、そちらに話を戻したい。ライブの音量は100 dB 程度らしい。上記の式でいうと、 $20 \log_{10}(p_0/p_{{\rm ref}, 0})$ である。この減衰の式がどう影響するのかを見ていきたい。基準となる距離が 1m でこの時の音圧レベルが 100dB とすると、上記の計算（古典減衰・分子緩和）をすべて合算すると、以下のグラフのような距離減衰になる。ただし、これは「1mの地点ですべての周波数が同じ100dBで出ている」という単純化をしたものであり、実際の音響システムでも周波数ごとの出力特性は一様とは限らない。したがって以下の議論は「空気中の伝搬による減衰の違い」だけを純粋に取り出したものであり、実際の音源特性を反映したものではない点に注意されたい。

この図は各周波数における横軸は対数距離、縦軸はデシベルとしたときのグラフである。ただし、基準となる温度と圧力はそれぞれ $T_0 = 293.15 {\rm K}$ ( 20℃ ) , $p_0 = 1013 {\rm hPa}$ (大気圧) とし、考える温度は夏なので $T = 308.15 {\rm K}$ (35℃ ) で圧力は大気圧 $p_m = 1013 {\rm hPa}$ で、相対湿度 $RH$ は 75％であるとした。
さきほどの分類でいうと、全く補正をしていない Z 特性音圧のグラフである。音楽ライブは 100dB らしいので 1m という基準地点では100dB になっていて、囁き声程度になるのが 20dB らしいのでそこを下限にしている。ただし、実際に聞こえるかどうかは周波数、背景騒音、個人差、測定条件によって変わる。
楽器や声の周波数のおおまかな目安として、サブベースは 30 ～ 60 Hz、ベース、バスドラム、男性声の低い成分は 100～500 Hz、人の声、ギター、ピアノの主成分は 500 Hz〜2 kHz、声の輪郭、アタック感、音の抜けは 2〜4 kHz、シンバル、トライアングル、小鳥の声、空気感は 4kHz ～ 20 kHz くらいのようなので [10]、そこから振動数をピックアップした。
この図から 15kHz では 200m 程度で 20dB になっていて、2kHz では約1km で20dB になることが分かる。一方で、45Hz という低周波数の重低音では 10km にならないと20dB 程度に落ちないということが分かる。これは、古典減衰が周波数の二乗に比例し、緩和吸収も低周波側では二乗に近く増加するため、高周波ほど空気吸収を受けやすいからである。

では、聴覚についてA特性音圧レベルを使って補正するとどうなるのか？
それをプロットしたのが次の図である。

図の横軸は対数距離で、縦軸は音のA特性音圧 dB(A) である。この図から20 dB において一番距離が届くのが 300Hz の音で距離は 3km 程度であり、そこから 500 Hz の約2.5 km、1000 Hz の約 1.6 km と 1500Hz の約 1.4 km、2kHz と 100Hz の約 1 km へと続いていく。A特性が入ると低音ほど影響を受け、45 Hz は音源近くで 100dB あっても、A特性では 68dB 程度に評価される、この単純なモデルでは 250 m 程度で 20 dB まで落ちることが分かる。
最初に話をしたライブで低音が話題になったように、20dB という囁き声程度であるなら300Hz という中低音が一番遠くまで届くという結果になった。これは楽器や声でいうと、ベース、バスドラム、男性声の低い成分にあたる。苦情で言われる「重低音」が、30〜60 Hz のサブベースだけでなく、100〜300 Hz 程度の低音・中低音を含んでいた可能性もある。

ただし、この図はあくまでも聴覚を通して認識できる音に関するものであり、実際の物理的な音の振動についてはZ特性のグラフのままであり、実際の物理的な音圧変動については、Z特性のグラフに示したように、低音側では空気吸収による追加減衰が小さいことに注意が必要である。
つまり、A特性では低周波が大きく割り引かれるが、Z特性で見た物理的な音圧変動そのものが消えるわけではない。低周波音では、振動が長距離まで届くことで、耳で明瞭に聞こえる音だけでなく、建物の共振や振動感として問題になることもある。

ちなみに、この結果はライブだけでなく他の現象にも応用できる。ハイキングや山登りでやまびこを試すとき、この図だけで厳密な限界距離は決められないが、数百 m 程度を越えると明瞭なやまびことして聞き取りにくくなる理由は理解できる。

それ以上離れると、ヤッホーという音が減衰して聞こえなくなってしまう。(もちろん、この距離は山の地形や植生などに影響を受けて実際は小さい場合もあるし、100dB 以上の大声を出せるなら少しは距離が伸ばせるかもしれない)

ここでは古典減衰・緩和による減衰・A特性音圧レベルを考慮したが、実際には風の向きと速度・周囲の建物とそれに伴う回折・音響機器の指向性等の性能やその配置など音の聞こえ方に影響を与えるような要素が複数存在することに注意してほしい。

まとめ

この記事では、音波の距離減衰を、古典的な粘性・熱伝導による減衰、酸素・窒素分子の緩和による減衰、人間の聴覚を近似的に反映するA特性音圧レベルの三つに分けて考えた。Z特性、つまり物理的な音圧レベルで見ると、低周波ほど空気吸収を受けにくく、遠方まで残りやすい。一方で、A特性を入れると、極端な低周波は人間の聴覚感度の低さによって大きく割り引かれる。その結果、今回の単純化した条件では、20 dB(A) 付近まで残る距離は 300 Hz 付近で最大になった。したがって、野外ライブで問題になる「重低音」は、単に低ければ低いほど聞こえるというより、物理的な減衰の少なさと人間の聴覚感度のバランスが合わさった結果だと考えられる。ただし、これは気体中の音波伝搬という側面だけを切り出した結果である点には注意したい。実際の苦情事例では、低域スピーカーの指向性が高域に比べて広い（会場外に音が回り込みやすい）ことや、低周波音が建物の壁や窓を振動させて伝わる固体伝搬・共振の効果も無視できないと考えられる。

参考文献

[1] 巽友正. 流体力学. 培風館
[2] L. D. Landau and E. M. Lifschitz, Fluid Mechanics 2nd ed.
[3] The Engineering Toolbox, https://www.engineeringtoolbox.com/air-absolute-kinematic-viscosity-d_601.html
[4] J. Shang, T. Wu, H. Wang, C. Yang, C. Ye, R. Hu, J. Tao, and X. He, IEEE Access, vol. 7, pp. 136439. ( https://www.researchgate.net/publication/335896009_Measurement_of_Temperature-Dependent_Bulk_Viscosities_of_Nitrogen_Oxygen_and_Air_From_Spontaneous_Rayleigh-Brillouin_Scattering )
[5] Science.Tools, 空気の物性値, https://cattech-lab.com/science-tools/properties-air/
[6] W. G. Vincenti and C. H. Kruger, Jr, Introduction to Physical Gas Dynamics, Krieger Pub Co.
[7] S. L. Garrett, Understanding Acoustics An Experimentalist's View of Sound and Vibration, Springer.
[8] A. N. Rimell, N. J. Mansfield and G. S. Paddan , Ind Health. 2015;53(1):21-7 ( https://pmc.ncbi.nlm.nih.gov/articles/PMC4331191 )
[9] F. Lai, Z. Huang and F. Guo, Materials 2021, 14, 4356. ( https://www.mdpi.com/1996-1944/14/16/4356 )
[10] soundzone, 何となく知っている「周波数」、きちんと説明できますか？, https://www.soundzone.jp/staffblog/10128/

2026-07-14

アナレンマ:太陽が描く８の字軌道

物理学

アナレンマとは？

定時刻に同一地点から太陽の中心を毎日記録・観測すると、その軌跡は８の字を描くことが知られている。この現象をアナレンマ (Analemma) という。
それを具体的に映像化したのが以下の YouTube 動画 ("Analemma (3-year time-lapse)" URL : https://www.youtube.com/watch?v=Deli5COMJhs ) だ。

この動画はサンディエゴの天文台で記録されたもののようだ。
アナレンマはもともとは日時計やその土台を指す言葉だったようだが、今ではこの太陽の８の字の軌跡のことを指すようになっている。
この記事では、赤経・赤緯や均時差の公式を直接使うのではなく、地球の自転を ${SO(3)}$ の回転として扱い、太陽方向を地球とともに回転する像平面へ射影することで、アナレンマの形を再現する。

設定

さて、このアナレンマは暦などで使われる赤道座標系を用いて計算するのが一般的であるようだが、私は「定点カメラで見える8の字は、自転・公転・射影だけで説明できないだろうか」という疑問を持った。そこで、この記事では別の視点から幾何学的に考えて定点カメラの撮影をイメージして、地球とともに回転する像平面への射影だけで、8の字の軌跡を再現できるかを試みたい。次のような設定で考える。

まず、地球が太陽の周りを角振動数 $\Omega = 2\pi / (365 \times 24) \ ({\rm rad/h})$ の円軌道（本当は楕円軌道だが第一近似として離心率を無視し、 $e = 0$ の円軌道とみなす）で公転している。
同時に、地軸が黄道面法線に対して角度 $\alpha = 23.4\pi/180 \ ({\rm rad})$ で傾き、地軸まわりに角振動数 ${\omega_e = \frac{2\pi}{24} + \Omega \simeq \frac{2\pi}{23.934} }$ (rad/h) で自転しているとする。地球の自転について、1日は太陽が南中してから次に南中するまで（太陽日）だが、地球の自転そのものは公転分を差し引いた恒星日ベースで考える必要がある。
ここで地軸方向の単位ベクトルを ${{\bf n} = (\sin\alpha, 0, \cos\alpha)}$ とする。ただし、黄道面に直交する方向を z 軸としている。

地球の自転と公転の数学的扱い

この地軸まわりの回転をまずは考えると、地球の自転は $SO(3)$ の生成子を用いて

\begin{align*}
\exp \left(-{\rm i}\omega_e t {\bf J}\cdot {\bf n}\right)
\end{align*}

と書くことができる。まずは $\exp$ の中身を計算すると

\begin{align*}
-{\rm i}\omega_e t {\bf J}\cdot {\bf n} & = -{\rm i}\omega_e t \left(J_x\sin\alpha + J_z\cos\alpha \right) \\
& = \omega_et\begin{pmatrix}
0 & -\cos\alpha & 0 \\
\cos\alpha & 0 & -\sin\alpha \\
0 & \sin\alpha & 0
\end{pmatrix}
\end{align*}

となる。ここで、行列を $K$ とおく。

\begin{align*}
K = \begin{pmatrix}
0 & -\cos\alpha & 0 \\
\cos\alpha & 0 & -\sin\alpha \\
0 & \sin\alpha & 0
\end{pmatrix}
\end{align*}

このとき $K$ の特性方程式は $\lambda^3 + \lambda = 0$ であるので、ケーリー・ハミルトンの定理より $K^3 + K = 0$ となる。この性質をつかって整理すると、以下のロドリゲスの回転公式と呼ばれる公式が得られる。

\begin{align*}
\exp(\omega_e t K) = 1 + \sin(\omega_e t)K + (1-\cos(\omega_e t))K^2
\end{align*}

これを使って実際に回転行列を計算すると

\begin{align*}
\exp(\omega_e t K) =
\begin{pmatrix}
\sin^2\alpha + \cos^2\alpha\cos(\omega_e t) & -\cos\alpha\sin(\omega_e t) & \sin\alpha\cos\alpha(1-\cos(\omega_e t)) \\
\cos\alpha\sin(\omega_e t) & \cos(\omega_e t) & -\sin\alpha\sin(\omega_e t) \\
\cos\alpha\sin\alpha(1-\cos(\omega_e t)) & \sin\alpha\sin(\omega_e t) & \cos^2\alpha + \sin^2\alpha\cos(\omega_e t)
\end{pmatrix}
\end{align*}

ちなみに、この記事ではロドリゲスの公式を使って計算しているが、ユニタリ変換を用いて $U^\dagger K U = {\rm diag}(-{\rm i}, 0, {\rm i})$ と対角化してから、 $U^\dagger \exp(\omega_e t K) U = {\rm diag}(e^{-{\rm i}\omega_et}, 1, e^{{\rm i}\omega_et})$ から　 $\exp(\omega_e t K) = U {\rm diag}(e^{-{\rm i}\omega_et}, 1, e^{{\rm i}\omega_et})U^\dagger$ と求める方法もある。
また、 $R_y(\alpha)$ , $R_z(\omega_e t)$ をそれぞれ y軸まわりに角度 $\alpha$ だけ回転させる行列、z軸まわりに角度 $\omega_e t$ だけ回転させる行列とすると　 $R_y(\alpha)R_z(\omega_e t)R_y^{-1}(\alpha)$ は $\exp(\omega_e tK)$ と同じ行列になるので、このやり方でも良い。

さて、簡単のために緯度 $\theta$ の地点の単位法線を ${\bf n}_\theta = (\cos(\theta-\alpha), 0, \sin(\theta-\alpha))$ とする。これに地球の半径をかけると地球の中心から緯度 $\theta$ の地点へのベクトルになる。
本当は経度 $\lambda$ も入れて $(\cos(\theta-\alpha)\cos\lambda, \cos(\theta-\alpha)\sin\lambda, \sin(\theta -\alpha))$ とするべきだがこの時に $\lambda = 0$ という特殊な点を考えている。しかしながら、アナレンマでは自転を考えており経度の違いは、地球の自転位相の違い、すなわち観測時刻 $t_0$ のずれとして吸収できる。そのため、アナレンマの形そのものを見る目的では、経度を $\lambda=0$ に固定しても本質的な一般性は失われない。このベクトル ${\bf n}_\theta$ に回転行列を作用させると次を得る。

\begin{align*}
{\bf n}_\theta(t) =
\begin{pmatrix}
\sin\theta\sin\alpha + \cos(\omega_e t)\cos\theta\cos\alpha \\
\sin(\omega_e t)\cos\theta \\
\sin\theta\cos\alpha - \cos(\omega_e t)\cos\theta\sin\alpha
\end{pmatrix}
\end{align*}

ここで、太陽に対する地球の軌道を半径が $a$ の円で近似すると、地球の軌道は次のようにあらわされる。

\begin{align*}
{\bf r}_{\rm e}(t)= a
\begin{pmatrix}
\cos(\Omega t) \\
\sin(\Omega t) \\
0
\end{pmatrix}
\end{align*}

ここで、黄道面に垂直な方向をz成分に取っている。これより, 緯度 $\theta$ の地点の位置は

\begin{align*}
{\bf r}(t) & = {\bf r}_{\rm e}(t) + R_\oplus {\bf n}_\theta(t) \\
& = a
\begin{pmatrix}
\cos(\Omega t) \\
\sin(\Omega t) \\
0
\end{pmatrix}
+ R_\oplus
\begin{pmatrix}
\sin\theta\sin\alpha + \cos(\omega_e t)\cos\theta\cos\alpha \\
\sin(\omega_e t)\cos\theta \\
\sin\theta\cos\alpha - \cos(\omega_e t)\cos\theta\sin\alpha
\end{pmatrix}
\end{align*}

となる。ここで $R_\oplus$ は地球の半径である。

射影平面の設定

アナレンマは同一時刻に同一地点から見た太陽の位置なので、基準時刻 $t_0$ における太陽方向 ${\bf s}(t_0)$ に垂直な平面をカメラの像平面として固定し、その平面を地球の自転とともに回転させる。
すなわち、ある時間 $t_0$ における地球から見た太陽方向の単位ベクトル ${\bf s}(t_0)$ を考え、それに直交する平面をなす2つの単位ベクトルを考える。これは、カメラを固定して1年間撮影し続けるような状況になっている。この2つのベクトルを地面に固定して地球と同様に回転して、時間 $t_n = 24n + t_0$ (h) においてその回転したあとの平面を射影した図形を考える。すると、その軌道はアナレンマになっているはずである。これはちょうどカメラのスクリーンを射影平面としているようなものである。
このような考えのもと、まずは地球から見た太陽方向の単位ベクトル ${\bf s}(t)$ を考える。

\begin{align*}
{\bf s}(t) & = \frac{-{\bf r}(t)}{|{\bf r}(t)|} = \frac{ -({\bf r}_{\rm e}(t) + R_\oplus {\bf n}_\theta(t))}{| {\bf r}_{\rm e}(t) + R_\oplus {\bf n}_\theta(t)|} \simeq \frac{-{\bf r}_{\rm e}(t)}{|{\bf r}_{\rm e}(t)|} =
\begin{pmatrix}
-\cos(\Omega t) \\
-\sin(\Omega t) \\
0
\end{pmatrix}
\end{align*}

ここで地球の軌道の半径 $a$ は1天文単位で約 $1.5 \times 10^{11}$ m (正確には 149597870700 m ) で、地球の半径 $R_\oplus$ は約 $6.4 \times 10^6$ m ( 6378137m) なので $R_\oplus {\bf n}_\theta(t)$ を無視した。
さて、アナレンマの同一時刻 $t_0$ での平面の基準として ${\bf c} = {\bf s}(t_0)$ とおく。このベクトルに直交するものとして、黄道面に直交する ${\bf e}_v = (0, 0, 1)$ がある。この2つのベクトルに直交する単位ベクトルは ${\bf e}_h = {\bf c}\times{\bf e}_v = (-\sin(\Omega t_0), \cos(\Omega t_0), 0)$ となる（下図参照）。
ここでは簡単のため、像平面の縦方向を黄道面法線方向に取っている。実際のカメラの縦方向に近づけるには、基準時刻の天頂方向や北向き方向を像平面へ射影して縦軸を定める必要がある。
この ${\bf e}_v$ , ${\bf e}_h$ を $\exp\left(\omega_e t {\bf K}\right)$ で回転させていったときに、それらのなす平面に $\bf{s}(t)$ を射影させるとアナレンマが描画される。

射影に使うベクトルの概略。緑のベクトルが $\bf c$ , ${\bf e}_v$ , ${\bf e}_h$ である。 $\bf c$ は太陽方向、 ${\bf e}_v$ は黄道面の法線方向、 ${\bf e}_h$ は表から裏の方向になっている。

射影に使うベクトルの概略。緑のベクトルが $\bf c$ , ${\bf e}_v$ , ${\bf e}_h$ である。 $\bf c$ は太陽方向、 ${\bf e}_v$ は黄道面の法線方向、 ${\bf e}_h$ は表から裏の方向になっている。

射影すべき平面は地球とともに回転する。このとき時刻 $t$ での ${\bf e}_v$ , ${\bf e}_h$ , $\bf c$ をそれぞれ ${\bf e}_v(t)$ , ${\bf e}_h(t)$ , ${\bf c}(t)$ とすると次のようになる。

\begin{align*}
{\bf e}_v(t) & = \exp(\omega_e(t-t_0){\bf K}){\bf e}_v = \exp(\omega_e(t-t_0){\bf K})
\begin{pmatrix}
0 \\
0 \\
1
\end{pmatrix}
\\
& =
\begin{pmatrix}
\sin\alpha\cos\alpha\left(1-\cos\omega_e(t-t_0)\right) \\
-\sin\alpha\sin\omega_e(t-t_0) \\
\cos^2\alpha+\sin^2\alpha\cos\omega_e(t-t_0)
\end{pmatrix}\\
\\
{\bf e}_h(t) & = \exp(\omega_e(t-t_0){\bf K}){\bf e}_h = \exp(\omega_e(t-t_0){\bf K})
\begin{pmatrix}
-\sin\Omega t_0 \\
\cos\Omega t_0 \\
0
\end{pmatrix}
\\
& =
\begin{pmatrix}
-(\sin^2\alpha+\cos^2\alpha\cos\omega_e(t-t_0))\sin\Omega t_0 - \cos\alpha\sin\omega_e(t-t_0)\cos\Omega t_0 \\
-\cos\alpha\sin\omega_e(t-t_0)\sin\Omega t_0 + \cos\omega_e(t-t_0)\cos\Omega t_0 \\
-\cos\alpha\sin\alpha(1-\cos\omega_e(t-t_0))\sin\Omega t_0 + \sin\alpha\sin\omega_e(t-t_0)\cos\Omega t_0
\end{pmatrix}
\\
{\bf c}(t) & = \exp(\omega_e(t-t_0){\bf K}){\bf c} = \exp(\omega_e(t-t_0){\bf K})
\begin{pmatrix}
-\cos\Omega t_0 \\
-\sin\Omega t_0 \\
0
\end{pmatrix}
\\
& =
\begin{pmatrix}
-(\sin^2\alpha + \cos^2\alpha\cos\omega_e(t-t_0))\cos\Omega t_0 + \cos\alpha\sin\omega_e(t-t_0)\sin\Omega t_0 \\
-\cos\alpha\sin\omega_e(t-t_0)\cos\Omega t_0 - \cos\omega_e(t-t_0)\sin\Omega t_0 \\
-\cos\alpha\sin\alpha(1-\cos\omega_e(t-t_0))\cos\Omega t_0 - \sin\alpha\sin\omega_e(t-t_0)\sin\Omega t_0
\end{pmatrix}
\end{align*}

ここで $\exp(\omega_e(t-t_0){\bf K})$ を作用させているのは $t = t_0$ で元のベクトルに一致させるためである。
原点から太陽方向 ${\bf s}(t)$ へ伸ばした直線 ${\bf X} = \lambda {\bf s}(t)$ が、 ${\bf c}(t)\cdot{\bf X}=1$ で定義される像平面と交わる点を考えると、交点は ${\bf X} ={\bf s}(t)/({\bf c}(t) \cdot {\bf s}(t))$ となる。これを ${\bf e}_h$ , ${\bf e}_v$ 方向に分解すると、アナレンマの x, y 成分は次のように計算される。

\begin{align*}
x(t) & = \frac{{\bf e}_h(t)\cdot {\bf s}(t)}{{\bf c}(t) \cdot {\bf s}(t)} \\
\\
y(t) & =\frac{{\bf e}_v(t)\cdot {\bf s}(t)}{{\bf c}(t) \cdot {\bf s}(t)} \\
\end{align*}

これに先ほど計算したベクトルの値を代入すると良い。そして、嬉しいことにここまでの計算は全て手で計算することができる。つまり、高緯度では日照条件のチェックが大変だが、円軌道を仮定して（ $\sin\alpha$ , $\cos\alpha$ , $\sin\theta$ や $\cos\theta$ などの値を知って）いればアナレンマのプロットまで手計算ですることができるということである。

円軌道のアナレンマ

このようにして得られた射影平面上の点 $(x(t), y(t))$ を位置 ${\bf r}(t) = {\bf r}_{\rm e}(t) + R_\oplus {\bf n}_\theta(t)$ において太陽が地表を照らしているという条件のもとで $t_n = 24n + t_0$ ( $n = 0, 1, \cdots, 364$ )についてプロットすると次のような形になる。

このときに $t = 0$ (h) から $t = 24$ (h) までの間で日照条件 ${\bf s}(t)\cdot {\bf n}_\theta(t) > 0$ の下で、この内積が最大になる時間を南中時間とし、 $t_0$ を南中時間になるように選んでいる。
これらの図から、だいたい日本の緯度くらいの北緯35度での太陽の同一時刻・同一地点での軌道が８の字の形になることが分かる。上下で円の大きさが微妙に違っているが、これは射影の分母 ${\bf c}(t)\cdot{\bf s}(t)$ が季節によって変化し、同じ角度差でも像平面上での倍率が変わるためだと考えられる。実際、赤道上ではアナレンマを射影する平面を地球の接平面にして南中時刻を基準にプロットすると以下のように円の大きさがほぼ等しいものが出てくる。（ただし、緯度が高くなるとやはり歪みが出る）

北緯90度では8の字の片方（北半分）しか描画されていないが、これは半年間太陽が地平線の下に隠れる「極夜」が発生するためである。計算上は8の字の全体が存在するが、観測可能な範囲に限定するとこのような形になる。よくアナレンマの軌道の構成要素として地球の軌道が楕円軌道であることと、地軸の傾きがあることがあげられるが、8の字の形になることは楕円軌道を仮定せずとも定性的には再現できる。ただし、動画のようなアナレンマの形になるには楕円軌道を考える必要がある。

楕円軌道のアナレンマ

楕円軌道の場合においてもこれまでに記述してきた方法はあまり変化しない。大きく変わるのは地球の軌道が $(a\cos\Omega t, a\sin\Omega t, 0)$ から極座標系を使った $(r(t) \cos\phi(t), r(t) \sin\phi(t), 0)$ となることである。
ケプラー方程式を使うと、楕円軌道の極座標における距離 $r$ と動径 $\phi$ は次のように書ける。

\begin{align*}
\Omega t & = \psi - e\sin\psi \\
r & = a(1 - e\cos\psi) \\
\tan\frac{\phi}{2} & = \sqrt{\frac{1+e}{1-e}}\tan\frac{\psi}{2}
\end{align*}

で与えられる。ここで、 [\tex: e], $\psi$ はそれぞれ離心率と離心近点角である。なお、ここでは簡単のため近日点方向を $x$ 軸に取っている。現実の地球のアナレンマにより近づけるには、近日点黄経に相当する位相 $\varpi$ を加えて $\lambda(t)=\phi(t)+\varpi$ とする必要がある。
最初の式にニュートン法を用いることで数値的に $\psi$ を求めてそこから $\phi$ を求めることでアナレンマが計算できる。
以下はその疑似コードである。

for n = 0,...,364:
    t = 24*n + t0
    solve Kepler equation
    compute true anomaly
    compute sun direction s(t)
    compute rotating image basis c(t), e_h(t), e_v(t)
    if s(t) · n_theta(t) >= 0:
        x = (e_h(t) · s(t)) / (c(t) · s(t))
        y = (e_v(t) · s(t)) / (c(t) · s(t))

その結果をプロットしたのが以下の図である。それぞれ赤道（0度)、北緯35度、北緯90度についてプロットしている。

赤道と北緯35度の場合には円軌道の場合とは異なり上側の方が小さくなっていることが見て取れる。これは楕円軌道にすることによって太陽黄経が等速でなくなる効果である。
また、円軌道だけを考えていた場合には、北緯90度でのアナレンマは片側の円しか見えなかったが、楕円軌道にすると下半分もわずかに見えることが分かる。これは北緯90度では、太陽赤緯が正の期間だけ太陽が地平線上に現れる。したがって、可視条件を課すと、円軌道モデルではアナレンマの一方のループだけがほぼ残る。楕円軌道を入れると、太陽黄経の進み方が非一様になるため、同じ可視条件を課してもプロットされる点の分布がわずかに変わり、射影平面上では下側の領域にも点が入り込む。これは極夜の期間にも太陽が見えるという意味ではなく、可視条件を満たす時期の点の配置が射影平面上で変化するということである。

また、離心率を火星の $e \sim 0.0934$ にすると水滴のような形状のアナレンマが出てくる。これは、赤道座標系などをもとにした火星のアナレンマの結果と整合的になっている。

一方で、この記事の方法は他の離心率でも全て上手くいくのかというとそうでもない。今考えているものは離心率が小さいので上手くいってるが、 $e \ge 0.68$ など比較的大きな離心率をもつような軌道については、太陽方向が像平面に近づき、 $\bf{c}\cdot {\bf s} \to 0$ となり、射影が特異になるためであると思われる。これは天体の軌跡そのものが発散するというよりは、固定した像平面への射影が、太陽方向が像平面とほぼ平行になるところで特異になるということである。

まとめ

この記事では赤道座標系を一切使わず、 $SO(3)$ と射影だけでアナレンマを再現した。この手法の結果、8の字の軌道自体は定性的に再現することが分かった。また、この手法は地球の公転軌道を円軌道で近似すると解析的に手計算することができるという利点がある。そして、地球が太陽の周りを楕円軌道で公転する効果を入れると動画のような軌道のアナレンマになることが分かった。また、この記事の射影を使った方法は離心率が高い場合には、太陽方向が基準像平面に対して大きくずれるため射影平面との角度がきつくなるという限界が生じることが分かった。機会があれば赤道座標系などの暦の計算で使われるような座標系を用いてアナレンマの計算に挑戦したい。

参考文献

国立天文台
暦計算室暦WIKI 恒星時 (Sidereal Time)( https://eco.mtk.nao.ac.jp/koyomi/wiki/B9B1C0B1BBFE.html )
K. M. Lynch and F. C. Park, "MODERN ROBOTICS MECHANICS, PLANNING, AND CONTROL"
H. Goldstein, C. P. Poole and J. Safko, "Classical Mechanics 3rd edition"

2026-07-14

フーリエ変換の固有値問題について

この記事は物理を学んだ人がフーリエ変換の固有値について考えた記事です。個人の思考をトレースしている部分が多く、あまりスッキリしない部分もあるかもしれないですが、それも味と思っていただければと思います。

フーリエ変換の固有値

フーリエ変換を四回すると恒等変換 id になる

実はあまり知られてない（かどうかわからないけど）こととして、関数にフーリエ変換を4回作用させると元の関数に戻るというものがある。この記事では、この性質を使ってフーリエ変換の固有値問題について考えたい。
以下では関数 $f$ のフーリエ変換を ${\mathcal F}[f]$ で表す。その定義は

\begin{align*}
{\mathcal F}[f](\xi) = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^\infty f(x) e^{-{\rm i}x\xi} \ {\rm d}x
\end{align*}

である。文献によっては積分の前の係数は $1/\sqrt{2\pi}$ 以外に $1/(2\pi)$ や $1$ にする定義もあるが、対称性のために $1/\sqrt{2\pi}$ としている。さて、数学的な厳密性はさておき、フーリエ変換を4回すると元の関数に戻ることは次のように示すことができる。まずは二回フーリエ変換を作用させると

\begin{align*}
{\mathcal F}^2[f](y) & = \frac{1}{(2\pi)} \int_{-\infty}^\infty e^{-{\rm i} ky} \int_{-\infty}^\infty f(x) e^{-{\rm i} kx} \ {\rm d}x \ {\rm d}k \\
& = \frac{1}{2\pi} \int_{-\infty}^\infty \int_{-\infty}^\infty f(x) e^{-{\rm i}k(x+y)} \ {\rm d}x \ {\rm d}k \\
& = \int_{-\infty}^\infty f(x) \delta(x+y) \ {\rm d}x = f(-y)
\end{align*}

となる。ここで、ディラックのデルタ関数を

\begin{align*}
\delta(x) & = \frac{1}{(2\pi)} \int_{-\infty}^\infty e^{-{\rm i}kx} \ {\rm d}k \\
\end{align*}

とした。 ${\mathcal F}^2[f](y) = f(-y)$ なので、 $h(y) = {\mathcal F}^2[f](y)$ とおき、さらに二回フーリエ変換を作用させると、

\begin{align*}
h(y) & = {\mathcal F}^2[f](y) = f(-y) \\
{\mathcal F}^4[f](x) & = {\mathcal F}^2[h](x) = f(x)
\end{align*}

となり、フーリエ変換を四回作用させると恒等変換になる。つまり、　 ${\mathcal F}^4 = {\rm id}$ となることがわかる。

固有値と仮定する固有関数の空間

さて、固有値問題を考える前に、急減少関数というものを導入する。 $\mathbb R$ 上の関数 ]tex: f] が次の性質をみたすとき、 $f$ を急減少関数という。

\begin{align*}
(1) & \quad f\in C^\infty(\mathbb{R}) \\
(2) & \quad \forall N \in \mathbb{Z}_{\ge 0} , \\
& \ |f|_N = \sum_{\alpha + k \le N}\sup_{x\in \mathbb{R}}(1 + |x|)^k\left|\left(\frac{\rm d}{{\rm d}x}\right)^\alpha f(x)\right| < \infty
\end{align*}

このような $f$ 全体を $\mathbb R$ 上のシュワルツ空間といい、 ${\mathcal S}(\mathbb R)$ と書く。
このシュワルツ空間のありがたい性質として、多項式がかかっていても、 $|x|\to \infty$ で $0$ に収束するというものがある。この性質のおかげで部分積分をしても $|x| \to \infty$ で端の影響を考えずに済む。急減少関数のよく知られた例としては、 $e^{-x^2/2}$ がある (そして実はこれがフーリエ変換の固有関数になっている）。
また、急減少関数 $f \in {\mathcal S}(\mathbb R)$ の良い性質として、急減少関数をフーリエ変換あるいは逆フーリエ変換しても急減少関数のままであるという性質がある ( ${\mathcal F}[f] \in {\mathcal S}(\mathbb R)$ , ${\mathcal F}^{-1}[f] \in {\mathcal S}(\mathbb R)$ )。したがって、フーリエ変換の固有値問題を考える際には、この空間に固有関数が存在するとして議論すると楽だということがわかる。

フーリエ変換の固有関数を $g \in {\mathcal S}(\mathbb R)$ とし、固有値を $\lambda$ とするとき、 ${\mathcal F}^4 = {\rm id}$ から $\lambda^4 = 1$ つまり $\lambda = \pm 1, \pm {\rm i}$ であることが分かる。

\begin{align*}
{\mathcal F}[g] & = \lambda g \\
g & = {\mathcal F}^4[g] = \lambda^4 g
\end{align*}

また、この性質から固有値によって固有関数の偶奇が以下のようにして分かる。

\begin{align*}
{\mathcal F}^2[g](y) & = \lambda^2 g(y) \\
{\mathcal F}^2[g](y) & = g(-y)
\end{align*}

より固有関数 $g$ は $\lambda^2 g(y) = g(-y)$ を満たすので、 $\lambda^2 = 1$ のとき $g$ は偶関数になり、 $\lambda^2 = -1$ のとき $g$ は奇関数になる。

固有関数の満たす方程式

さて、 $\lambda g = {\mathcal F}[g]$ という固有値方程式から他の固有値 $-\lambda, \pm {\rm i} \lambda$ が出てこないかを考えたい。まず、 $-\lambda$ については、次のように計算する。

\begin{align*}
\lambda \frac{{\rm d}}{{\rm d}\xi} g(\xi) & = \frac{1}{\sqrt{2\pi}} \frac{{\rm d}}{{\rm d}\xi} \int_{-\infty}^\infty g(x) e^{-{\rm i} \xi x} {\rm d}x \\
& = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^\infty g(x) \frac{\partial}{\partial \xi} e^{-{\rm i} \xi x} {\rm d}x \\
& = \frac{-{\rm i}}{\sqrt{2\pi}} \int_{-\infty}^\infty x g(x) e^{-{\rm i} \xi x} {\rm d}x \\
& = \frac{-{\rm i}}{\sqrt{2\pi}} \int_{-\infty}^\infty x g(x) \frac{1}{-{\rm i} \xi} \frac{\partial }{\partial x} e^{-{\rm i} \xi x} {\rm d}x \\
& = \frac{1}{\sqrt{2\pi}} \frac{1}{\xi} \int_{-\infty}^\infty x g(x) \frac{\partial }{\partial x} e^{-{\rm i} \xi x} {\rm d}x \\
& = \frac{1}{\sqrt{2\pi}} \frac{1}{\xi} (-1) \int_{-\infty}^\infty \frac{{\rm d}}{{\rm d}x} (x g(x)) e^{-{\rm i} \xi x} {\rm d}x \\
& = -\frac{1}{\xi}{\mathcal F}[g + xg']
\end{align*}

これを整理すると、最終的に次を得る。

\begin{align*}
\lambda \xi \frac{{\rm d}}{{\rm d}\xi} g(\xi) & = -{\mathcal F}[g + xg'](\xi) \\
-\lambda \xi \frac{{\rm d}}{{\rm d}\xi} g(\xi) & = {\mathcal F}\left[\frac{g}{2} + xg'\right](\xi) +\frac{1}{2} \lambda g(\xi) \\
-\lambda \left(\xi \frac{{\rm d}}{{\rm d}\xi} g(\xi) + \frac{1}{2} g(\xi) \right) & = {\mathcal F}\left[\frac{g}{2} + xg'\right](\xi)
\end{align*}

よって、 $\left(x\frac{\rm d}{{\rm d}x} + \frac{1}{2} \right ) g(x)$ は固有値 $-\lambda$ の固有関数になる。
次に $\pm {\rm i} {\lambda}$ が固有値になる場合を考える。

\begin{align*}
\lambda \xi g(\xi) & = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^\infty g(x) \xi e^{-{\rm i} \xi x} {\rm d}x = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^\infty g(x) {\rm i} \frac{\partial}{\partial x} e^{-{\rm i} \xi x} {\rm d}x \\
& =\frac{{\rm i}}{\sqrt{2\pi}} \int_{-\infty}^\infty g(x) \frac{\partial}{\partial x} e^{-{\rm i} \xi x} {\rm d}x = \frac{-{\rm i}}{\sqrt{2\pi}} \int_{-\infty}^\infty \frac{\rm d}{{\rm d} x} g(x) e^{-{\rm i} \xi x} {\rm d}x =-{\rm i} {\mathcal F}[g']
\end{align*}

よって、 ${\rm i}\lambda$ に関する式として次を得る。

\begin{align*}
{\rm i} \lambda \xi g = {\mathcal F}[g']
\end{align*}

$-\lambda$ の固有関数を求めた方法とは別のやり方で $\lambda g = {\mathcal F}[g]$ を微分すると

\begin{align*}
\lambda \frac{{\rm d}}{{\rm d}\xi} g(\xi) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^\infty g(x) \frac{\partial}{\partial \xi} e^{-{\rm i} \xi x} {\rm d}x = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^\infty g(x) (-{\rm i} x) e^{-{\rm i} \xi x} {\rm d}x = -{\rm i} \mathcal{F}[xg](\xi)
\end{align*}

となり、整理すると

\begin{align*}
{\rm i} \lambda g' = {\mathcal F}[xg]
\end{align*}

$xg$ と $g'$ のフーリエ変換の二つの式を組み合わせることで、 $\pm {\rm i}\lambda$ の固有値をもつ固有関数を得る。

\begin{align*}
{\mathcal F}[xg + g'] & = {\rm i} \lambda (\xi g + g') \\
{\mathcal F}[xg - g'] & = {\rm i} \lambda g' - {\rm i}\lambda \xi g = -{\rm i} \lambda (\xi g - g')
\end{align*}

この式は面白く、この $\left(x + \frac{\rm d}{{\rm d}x}\right)$ と $\left(x - \frac{\rm d}{{\rm d}x}\right)$ によって固有値が巡回的に移り変わっていくことが分かる。
例えば固有値 $\lambda$ を持つ固有関数 $g$ に $\left(x - \frac{\rm d}{{\rm d}x}\right)$ を一回作用させるごとに固有値は $-{\rm i}$ 倍されて
$\lambda \rightarrow -{\rm i} \lambda \rightarrow -\lambda \rightarrow {\rm i}\lambda \rightarrow \lambda$
のように移り変わっていく。逆に、固有値 $\lambda$ を持つ固有関数 $g$ に $\left(x + \frac{\rm d}{{\rm d}x}\right)$ を一回作用させるごとに固有値は $\rm i$ 倍されて
$\lambda \rightarrow {\rm i} \lambda \rightarrow -\lambda \rightarrow -{\rm i}\lambda \rightarrow \lambda$
のように移り変わっていく。

二次の関係から調和振動子へ

ここまでは、一階微分をみてきたので二階微分を考える。

\begin{align*}
\lambda g & = {\mathcal F}[g] \\
\lambda \frac{\rm d^2}{{\rm d}\xi^2} g(\xi) & = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^\infty g(x)\frac{\partial^2}{\partial \xi^2}(e^{-{\rm i} \xi x}) {\rm d}x \\
& = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^\infty g(x) \ (-x^2)e^{-{\rm i} \xi x} {\rm d}x \\
& = -{\mathcal F}[x^2 g]
\end{align*}

先ほどの場合と同様に $x^2 g$ のフーリエ変換と二階微分がつながったので、次は逆を考える。

\begin{align*}
\lambda \xi^2g(\xi) & = \xi^2{\mathcal F}[g](\xi) \\
& = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^\infty g(x) \xi^2 e^{-{\rm i} \xi x} {\rm d}x \\
& = -\frac{1}{\sqrt{2\pi}}\int_{-\infty}^\infty g(x) \ \frac{\partial^2}{\partial x^2}(e^{-{\rm i} \xi x}) {\rm d}x \\
& = -\frac{1}{\sqrt{2\pi}}\int_{-\infty}^\infty \left(\frac{\rm d^2}{{\rm d} x^2}g(x)\right)\ e^{-{\rm i} \xi x} {\rm d}x \\
& = -{\mathcal F}\left[\frac{\rm d^2}{{\rm d}x^2} g\right]
\end{align*}

先ほどと同様にこれらの二式を組み合わせると次を得る。

\begin{align*}
\lambda\left(-\frac{\rm d^2}{{\rm d}\xi^2}+\xi^2\right)g = +{\mathcal F}[x^2 g] - {\mathcal F}\left[\frac{\rm d^2}{{\rm d}x^2}g\right] = {\mathcal F}\left[\left(-\frac{\rm d^2}{{\rm d}x^2}+x^2\right)g\right]
\end{align*}

よって $g$ が固有値 $\lambda$ を持つ固有関数であるとき、 $\left(-\frac{\rm d^2}{{\rm d}x^2}+x^2 \right) g$ も固有値 $\lambda$ を持つ固有関数であることが分かる。
ここで微分作用素 $H$ として

\begin{align*}
H = -\frac{\rm d^2}{{\rm d}x^2}+x^2
\end{align*}

とおくと、これは調和振動子のハミルトニアンになっている。
さらに、 $\lambda\left( -\frac{\rm d^2}{{\rm d}\xi^2}+\xi^2\right)g = {\mathcal F} \left[\left(-\frac{\rm d^2}{{\rm d}x^2}+x^2 \right)g\right]$ は

\begin{align*}
H{\mathcal F} = {\mathcal F}H
\end{align*}

と表すことができ、ハミルトニアン $H$ とフーリエ変換 $\mathcal F$ が可換で同時固有関数をもつことを表している。これはフーリエ変換 $\mathcal F$ が調和振動子のハミルトニアン $H$ を不変にする変換であることを示している。また、 $\mathcal F$ と $H$ が可換であるので、調和振動子の固有関数を使って両者を同時に対角化する形で整理できる。 $H$ に対する固有値方程式は量子力学を勉強した人にはよく見知った微分方程式になる。

\begin{align*}
\left(-\frac{\rm d^2}{{\rm d}x^2}+x^2 \right) g = \alpha g
\end{align*}

この式は調和振動子の固有値問題の微分方程式になっていて、固有関数はエルミート関数 $H_n(x)e^{-x^2/2}$ であることが知られている。ここで $H_n(x)$ はエルミート多項式である。
つまり、エルミート関数は、 $H$ と $\mathcal F$ の共通の固有関数を与える。

フーリエ変換の固有関数が調和振動子の固有関数でもあることが分かると

\begin{align*}
{\mathcal F}[xg + g'] & = {\rm i} \lambda (\xi g + g') \\
{\mathcal F}[xg - g'] & = {\rm i} \lambda g' - {\rm i}\lambda \xi g = -{\rm i} \lambda (\xi g - g')
\end{align*}

という式の見方が変わってくる。
量子力学と対応させるように $p = -{\rm i}\frac{\rm d}{{\rm d}x}$ とおくと、 $x - \frac{\rm d}{{\rm d}x} = x - {\rm i}p$ 、 $x + \frac{\rm d}{{\rm d}x} = x + {\rm i}p$ となる。これらはそれぞれ $1/\sqrt{2}$ 倍すると生成・消滅演算子になる。消滅演算子は $a = \frac{1}{\sqrt{2}}( x + {\rm i}p)$ であり、生成演算子は $a^\dagger = \frac{1}{\sqrt{2}}( x - {\rm i}p)$ であった。
ここで、 $\mathcal F$ の固有関数 $g_n$ として $g_n(x) = H_n(x)e^{-x^2/2}$ を考えるとする。このとき、 $a^\dagger$ は調和振動子の固有状態 $|n\rangle$ を $\sqrt{n+1}|n+1\rangle$ にして、逆に $a$ は調和振動子の固有状態 $|n\rangle$ を $\sqrt{n}|n-1\rangle$ にしたのを思い出すと、規格化定数を省いた $g_n(x) = H_n(x)e^{-x^2/2}$ を用いてるので、次を得る。

\begin{align*}
\frac{1}{\sqrt{2}}\left(x - \frac{\rm d}{{\rm d}x}\right)H_n(x)e^{-x^2/2} & \propto H_{n + 1}(x)e^{-x^2/2} \\
\frac{1}{\sqrt{2}}\left(x + \frac{\rm d}{{\rm d}x}\right)H_n(x)e^{-x^2/2} & \propto H_{n - 1}(x)e^{-x^2/2} \\
\end{align*}

さて、これまで $H_n(x)e^{-x^2/2}$ がフーリエ変換の固有関数になると考えてきたが、実際に $g_0(x) = e^{-x^2/2} = H_0(x)e^{-x^2/2}$ がフーリエ変換の固有関数になるかを確かめる。 $H_0(x)e^{-x^2/2}$ がフーリエ変換の固有関数であれば、これまでに見てきたように $\frac{1}{\sqrt{2}}\left(x - \frac{\rm d}{{\rm d}x}\right)$ を作用させていくことで、他の固有値をもつ固有関数を生成することができる。

\begin{align*}
{\mathcal F}[g_0](\xi) & = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^\infty g_0(x)e^{-{\rm i}\xi x} \ {\rm d}x \\
& = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^\infty e^{-\frac{1}{2}x^2-{\rm i}\xi x} \ {\rm d}x \\
& = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^\infty e^{-\frac{1}{2}(x+{\rm i}\xi)^2-\frac{1}{2}\xi^2} \ {\rm d}x \\
& = \frac{e^{-\frac{1}{2}\xi^2}}{\sqrt{2\pi}}\int_{-\infty}^\infty e^{-\frac{1}{2}(x+{\rm i}\xi)^2} \ {\rm d}x \\
& = \frac{e^{-\frac{1}{2}\xi^2}}{\sqrt{2\pi}}\sqrt{2\pi} = e^{-\frac{1}{2}\xi^2} = g_0(\xi)
\end{align*}

よって ${\mathcal F}[g_0](\xi) = g_0(\xi)$ となるので $g_0(x) = e^{-x^2/2}= H_0(x)e^{-x^2/2}$ は固有値 $1$ の固有関数になっている。
$H_0(x)e^{-x^2/2}$ に生成演算子 $a^\dagger$ を作用させると、固有関数は
$g_0 \rightarrow g_1 \rightarrow g_2 \rightarrow g_3 \rightarrow g_4$
と変化していく。また、それに伴って固有値は
$1 \rightarrow -{\rm i} \rightarrow -1 \rightarrow {\rm i} \rightarrow 1$
と変化していく。
すなわち、 $a^\dagger$ を $n$ 回作用させると定数倍をのぞいて $H_n(x)e^{-x^2/2}$ になり、この固有関数の固有値は $(-{\rm i})^n$ になることが分かる。

母関数による計算

エルミート多項式の母関数を使うとエルミート関数がフーリエ変換の固有関数であることがすっきりと分かる。エルミート多項式の母関数表示は

\begin{align*}
e^{-t^2+2xt} = \sum_{n=0}^\infty \frac{1}{n!}H_n(x)t^n
\end{align*}

であった。ここに $e^{-x^2/2}$ をかけると

\begin{align*}
e^{-t^2+2xt-\frac{1}{2}x^2} = \sum_{n=0}^\infty \frac{1}{n!}H_n(x)e^{-\frac{1}{2}x^2}t^n
\end{align*}

この左辺をフーリエ変換すると形式的に以下のように式変形できる。

\begin{align*}
{\mathcal F}[e^{-t^2+2xt-\frac{1}{2}x^2}] & = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^\infty e^{-t^2+2xt-\frac{1}{2}x^2} \ e^{-{\rm i} \xi x} \ {\rm d} x \\
& = \frac{e^{-t^2}}{\sqrt{2\pi}}\int_{-\infty}^\infty e^{-\frac{1}{2}x^2+(2t-{\rm i}\xi)x} \ {\rm d} x \\
& = \frac{e^{-t^2}}{\sqrt{2\pi}}\int_{-\infty}^\infty e^{-\frac{1}{2}(x-(2t-{\rm i}\xi))^2} \ e^{\frac{1}{2}(2t-{\rm i}\xi)^2} \ {\rm d} x \\
& = e^{-t^2}e^{\frac{1}{2}(2t-{\rm i}\xi)^2} = e^{t^2-2{\rm i} t\xi}e^{-\xi^2/2} \\
& = e^{(-{\rm i} t)^2+2(-{\rm i} t)\xi}e^{-\xi^2/2} \\
& = \sum_{n=0}^\infty \frac{1}{n!}H_n(\xi)e^{-\frac{1}{2}\xi^2}(-{\rm i} t)^n \\
& = \sum_{n=0}^\infty \frac{1}{n!}(-{\rm i})^nH_n(\xi)e^{-\frac{1}{2}\xi^2}t^n
\end{align*}

先ほどの式と比較すると

\begin{align*}
{\mathcal F}[H_n(x)e^{-x^2/2}](\xi)= (-{\rm i})^n H_n(\xi)e^{-\xi^2/2}
\end{align*}

となっていることが分かる。つまり、 $H_n(x)e^{-x^2/2}$ はフーリエ変換の固有関数で、固有値 $(-{\rm i})^n$ を持っていることが分かる。

まとめ

この記事ではフーリエ変換 $\mathcal F$ の固有値問題を調べて、 ${\mathcal F}^4 = {\rm id}$ という性質から、固有値の値が $\{\pm1, \pm {\rm i}\}$ をとるということをみた。
さらに、 $H = -\frac{\rm d^2}{{\rm d}x^2} + x^2$ はフーリエ変換 $\mathcal F$ と可換であり、 $H$ の固有関数であるエルミート関数 $g_n(x) = H_n(x)e^{-x^2/2}$ は ${\mathcal F}[g_n] = (-{\rm i})^ng_n$ を満たす。したがって、エルミート関数はフーリエ変換を対角化する自然な基底になっている。

参考文献

「新・フーリエ解析と関数解析学」新井仁之
「演習形式で学ぶ特殊関数・積分変換入門」蓬田清

2025-06-29

統計物理からみる所得分布

物理学

この記事の目的
年収の指数分布はどこから来たのか？
二人稼ぎ世帯のモデル
まとめ

この記事の目的

このブログ記事を書こうと思ったきっかけは 2025年5月27日、AERAのWeb記事がSNSを中心に話題になったことである。
そのタイトルは「平均年収460万円は『実感』とはほど遠い？実は6割が平均以下の『真実』」というものである。このWeb記事は国税庁の民間給与実態統計調査において平均年収は460万円になっていて、その平均年収よりも下の年収帯の人が6割近くいるものであるという内容である。

dot.asahi.com

多くの人が「自分の年収は平均より下だったのか」と驚いたこのニュースに対し、神戸大学で教鞭をとっている統計物理学の専門家である西野先生が興味深いコメントを寄せていた。「指数分布f(x)=exp(-ax) で計算してみると平均以下の割合は a によらず　0.63212056 自由経済が適正に保たれている証拠である」（なお、「自由経済が適正に保たれている証拠である」というのは西野先生個人の解釈であり、本記事ではその是非については立ち入らないことにする。）

日本の給与「６割が平均以下」

指数分布 f(x)=exp(-ax) で計算してみると
平均以下の割合は a によらず 0.63212056
自由経済が適正に保たれている証拠である
— DMRG.JP_T.Nishino_西野友年 (@IRFDMRG) 2025年5月27日

では、この「63%」という数字は一体どこから来たのか？
まずは、この文の意味するところを見ていきたい（変数を $x$ から $m$ に変更している）。規格化された指数分布

$f(m)=a\exp(-am)$

に対して、平均値は

$E[X] = \int_0^\infty mf(m) = \frac{1}{a}$

となる。このとき、 $m$ 以下の値が占める割合を表している累積分布関数は

$F(m) = \int_0^m f(t) {\rm d}t = 1-\exp(-am)$

で与えられる。さて、平均以下の割合はこの累積分布を用いると

$F\left(E[X] \right) = 1-\exp(-1)=0.63212056$

であることがわかる。つまり、指数分布ではパラメータ $a$ の値に関わらず、平均以下の人は63%ほどになる。これが上記のツイートの意味するところである。
また、こうした話題では「平均値」だけでなく、「中央値」もよくよりあげられるので、ここで話をしておきたい。指数分布

$f(m)=a\exp(-am)$

の平均値は

$E[X]=\frac{1}{a}$

なので、記事のように平均年収が460万円であるなら $\frac{1}{a} = 460$ となる（単位は万円である）。

中央値を $m_{\rm med}$ とすると、累積分布関数 $F(m)$ は $1/2$ となる。

$F(m_{\rm med}) = \frac{1}{2}$

指数分布の累積分布関数は

$F(m) = 1-\exp(-am)$

であったので、これを使うと
$\frac{1}{2}= 1-\exp(-am_{\rm med})$
となり、これを整理すると

$\exp(-am_{\rm med}) = \frac{1}{2}$
$am_{\rm med} = \log(2)$

$m_{\rm med} = \frac{\log(2)}{a} \approx 0.6931\times 460 = 318.826$

となる。つまり指数分布において中央値はおよそ 318.8 万円程度になり、平均よりも低い値であることが分かる。これは高所得者が平均値を引き上げるために中央値と平均値がずれるためである。
さて、指数分布での性質を見てきたがここで自然に次の問いが生まれる。
「なぜ指数分布が出てくるのか？」
「日本の年収分布は指数分布でどこまで説明できるのか？」
年収分布の背後には、統計物理的なモデルや社会構造の影響が隠れている。本記事では、その数理的背景を順を追って解き明かしていきたいと思う。

年収の指数分布はどこから来たのか？

結論からいうと指数分布は A. Dragulescu and V. M. Yakovenko による "Statistical mechanics of money" という論文がもとになる。

arxiv.org

この論文では、総所得が保存される経済（クローズドな系）を仮定して、統計力学のエネルギー分布と同じ考え方で分布を以下のように導出している（論文の数式展開を私のレベルでも分かるように詳しく書いくことにする）。所得 $m$ の分布が $P(m_i+m_j) = P(m_i)P(m_j)$

という独立性と加法性から

$P(m)=\exp\left(-m/T\right)$

という式が得られる。ここで、 $T$ は平均年収で $T = \sum_i m_i/N$ となっていて、統計力学における温度に対応している。 $T$ が大きいほど経済における所得の広がり（貧富の差）が大きいということになる。
このやり方の他にもエントロピー最大化原理を用いて、年収分布を求めることができる[1]。エントロピーを次のようにとり、
$S =- \int_0^\infty p(m) \log p(m) \ {\rm d}m-\lambda_1\left(\int_0^\infty p(m) \ {\rm d} m -1\right) - \lambda_2\left(\int_0^\infty mp(m) \ {\rm d} m -\langle m \rangle \right)$

$p \ , \ \lambda_1 \ , \ \lambda_2$ についてのオイラー・ラグランジュ方程式を解けば求める分布関数 $p(m)$ が求まる。オイラー・ラグランジュ方程式から
$-\log p(m) - 1 - \lambda_1 - \lambda_2 = 0$
これを整理すると
$p(m) = \exp(-1 - \lambda_1 - \lambda_2 m )$

となり、 $\lambda_1$ の拘束条件から
$\int_0^\infty p(m) \ {\rm d} m = e^{-(1+\lambda_1)}\int_0^\infty e^{-\lambda_2 m} \ {\rm d} m = \frac{e^{-(1+\lambda_1)}}{\lambda_2} = 1$ .

よって $e^{-(1+\lambda_1)} = \lambda_2$ を得る。また、 $\lambda_2$ の拘束条件から
$\int_0^\infty mp(m) \ {\rm d} m = \lambda_2\int_0^\infty m e^{-\lambda_2 m} \ {\rm d} m = \frac{1}{\lambda_2} = \langle m\rangle$
となり、分布が指数分布
$p(m) = \frac{1}{T}\exp(-m/T)$
となることがわかる。ここで $T = \langle m\rangle$ とおいた。

さて、この指数分布が実際の日本の年収データをどの程度一致するのかを比較するため、国税庁の「民間給与実態統計調査」（AERAの記事中で引用）と指数分布を比較したのが次の図である。

この図の横軸は所得階級をあらわし、縦軸はその所得帯をもつ人の割合(％)である。この図では所得区間 $[m_1, m_2)$ に含まれる人の割合と、その区間での指数分布の割合

$\int_{m_1}^{m_2} \frac{1}{T}\exp(-m/T) {\rm d} m = \exp(-m_1/T)-\exp(-m_2/T)$

を比較している。 $T$ はこの統計調査の結果から460万とわかっているので、その値を用いている。この図からは実データと指数分布の間に大きなズレがあり、指数分布によるモデル化が現実を記述できていないことがわかる。国税庁の民間給与実態統計調査は源泉徴収のある民間事業所を対象にしており、自営業・フリーランス・公務員は対象になっていない。そのため、標本にある種の偏りがあるのでこのように分布のズレが出ていると考えられる。

二人稼ぎ世帯のモデル

では、指数分布は完全に間違いなのか？必ずしもそうではない。
A. Dragulescu, V. M. Yakovenko による "Evidence for the exponential distribution of income in the USA" という別の論文で、一人あたりの年収が指数分布にしたがい、家族に稼ぎ手が二人いる場合、アメリカの世帯収入がガンマ分布でよく記述できることが示されている[2]。

arxiv.org

この二人稼ぎモデルの分布の導出を見ていく。確率変数 $X_1$ 、 $X_2$ が指数分布に従いかつ独立であるとき $X_1, X_2 \sim Exp(1/T)$ 、その和 $X_1+X_2$ が従う確率分布 $p_2(m)$ は指数分布の畳み込みになる。

$p_2(m) = \int_0^m p_1(m-m')p_1(m') {\rm d}m'$

ここで $p_1(m)$ は1人の所得分布をあらわしていて、この場合は指数分布である。
$p_1(m)=\frac{1}{T}\exp(-m/T)$

二人稼ぎの分布 $p_2(m)$ を計算すると

$p_2(m) = \int_0^m p_1(m-m')p_1(m') {\rm d}m'=\frac{1}{T^2}\int_0^m e^{-(m-m')/T}e^{-m'/T} {\rm d}m' = \frac{m}{T^2} e^{-m/T}$

となる。これは形状パラメータ2のガンマ分布になっている。

では、日本でこの分布は良い近似になっているのだろうか？この疑問への答えが、次の図である。

この図は令和5年の国民生活基礎調査の結果と二人稼ぎ世帯モデルのガンマ分布の理論値を比較したものである。横軸は世帯所得を表していて、縦軸はその所得をもつ世帯が占める割合を表している。令和5年の国民生活調査では平均世帯収入は524.2万円となっている。

この図では所得区間 $[m_1, m_2)$ に含まれる人の割合と、その区間でのガンマ分布の割合

$\int_{m_1}^{m_2} \frac{m}{T^2}\exp(-m/T) {\rm d} m = \frac{(m_1 + T)\exp(-m_1/T)-(m_2+T)\exp(-m_2/T)}{T}$

を比較している。ガンマ分布の平均値は

$E[X] = \int_0^\infty \frac{m^2}{T^2} e^{-m/T} = 2T$
であるので $2T = 524.2$ （万円）とおいている。

実データとガンマ分布を比較すると、ともに一峰性があり、先ほどよりも記述できているように見える。低所得帯で理論値よりも実データが大きいのは1人稼ぎ世帯が存在すること、高所得帯(2000万以上）で理論値よりも実データが大きいのは高所得帯では分布がべき的になるパレート分布 ( $p(m)\approx m^{-\alpha}$ ) になっていることを示唆している。
令和5年の世帯所得と二人稼ぎモデルのガンマ分布がどの程度近いのかをより直感的に見るために、ローレンツ曲線とジニ係数を導入する。

ローレンツ曲線というのは一定のグループへの集中度合いを視覚化したものである。世帯所得の集中度合いを視覚化するために、横軸を累積相対度数、縦軸を累積相対所得としてプロットした曲線となっている。数式で表すと (x, y) 座標が
$x = \int_0^r f(s) {\rm d}s \quad , \quad y =\frac{1}{\mu} \int_0^r s f(s) {\rm d} s$

となる曲線である。ここで、 $f(s)$ は分布関数、 $\mu$ は平均 $E[X]$ をあらわす。 $x$ は累積分布関数になっていることがわかる。また、 $y = x$ となる傾き1 の直線が完全平等線と呼ばれ、この直線にローレンツ曲線が近いほど所得が平等に分配されていて、逆に離れているほど不平等ということになる。

この不平等の度合いを数値化したものがジニ係数である。ジニ係数はローレンツ曲線と完全平等線とで囲まれた領域を A, ローレンツ曲線、 $y= 0$ 、 $x = 1$ で囲まれた領域を B とした時に、ジニ係数は $G := A/ (A+ B)$ で定義される。 $A+ B = 1/2$ なので $G = 2A$ である(下図参照）。

ここで、ガンマ分布

$p_2(m) = \frac{m}{T^2}e^{-m/T}$

のローレンツ曲線について考える[2]。媒介変数 $\tilde{r} \quad (0\le \tilde{r} < \infty)$ を用いると、ローレンツ曲線のx座標とy座標はそれぞれ次のように書ける。
$x(\tilde{r}) = 1-(\tilde{r}+1)e^{-\tilde{r}} \quad , \quad y(\tilde{r}) = 1-\left(\frac{1}{2}\tilde{r}(\tilde{r}+2) + 1\right)e^{-\tilde{r}}$

ここで $\tilde{r} = m/T$ となっている。ジニ係数はローレンツ曲線と完全平等線とで囲まれた領域 A を2倍したものだったので、
$G = 2\int_0^1 (x-y) dx = 2\int_0^\infty (x(\tilde{r})-y(\tilde{r})) \frac{{\rm d} x}{{\rm d} \tilde{r}}{\rm d} \tilde{r}$

$G = 2\int_0^\infty \frac{1}{2}\tilde{r}^2e^{-\tilde{r}}\cdot \tilde{r}e^{-\tilde{r}} {\rm d} \tilde{r}=2\frac{1}{2^5}\Gamma(4)=\frac{3}{2^3}$

となり $G = 3/8 = 0.375$ となる。

ちなみに、完全平等線はどのような分布になるのかについても触れておきたい。完全平等線は $y=x$ なので
$\int_0^rf(s){\rm d} s = \frac{1}{\mu}\int_0^rsf(s){\rm d}s$

$\int_0^r\left(1-\frac{s}{\mu}\right)f(s){\rm d} s = 0$

となる。これを任意の $r$ で満たすのはディラックのデルタ関数 $f(s) = \delta(s-\mu)$ になる。分布関数がデルタ関数になる状況はすべての人が同じ所得を得ている状況をあらわしている。

さて、令和5年のローレンツ曲線とジニ係数を比較すると次のようになる。

この図は世帯所得のデータから計算したローレンツ曲線とガンマ分布のローレンツ曲線を比較したものである。横軸は累積相対度数をあらわし、縦軸は累積相対所得をあらわしている。オレンジ色の曲線が令和5年の世帯所得のローレンツ曲線、水色の曲線がガンマ分布のローレンツ曲線、そして破線が完全平等線である。この図からかなり近い形にはなっているものの完全には一致していないことがわかる。

ではどの程度近いのだろうか？

所得分布のデータを各階級で階級値で代表させる近似 (ただし、2000万円以上の階級の代表値は便宜的に 2100万円とした）に基づいて階級近似ジニ係数を計算すると 0.405 となる。一方、先ほど計算したように形状パラメータ 2 のガンマ分布のジニ係数は0.375 となっている。これは 7.41% 程度の相対誤差がある。この値はおおむね良い近似になっていると言えるものになっている。二人稼ぎモデルは個人の収入が指数分布に従うという仮定があったので、ガンマ分布が現実の世帯収入を近似しているということは個人の収入も指数分布に従うことの傍証になっている。一方で、低所得帯と高所得帯でのずれが、ジニ係数の相対誤差をうみだしている。これは以上のようにシンプルな仮定から導出したボルツマン分布よりも現実を記述する分布が存在するかもしれないことを示唆している。

まとめ

統計力学を用いると所得分布は指数分布になることが分かったが、国税庁の民間給与実態統計調査のデータとは一致はしていないことが分かった。これは、調査対象にある種の偏りがありフリーランスの人などが含まれてないためと考えられる。一方で、世帯所得分布をみたときには、指数分布を仮定した二人稼ぎモデルのような分布になることが分かった。このズレがどの程度であるかを見るためにローレンツ曲線とジニ係数を比較した。その結果、相対誤差は7.41%で良い近似になっていることが分かった。
この記事ではほぼ論文を書き下してたような感じだったので機会があれば貯蓄や課税などの現実に近い仮定するとどのような分布になるのかなど、他のモデルについてもみてみたい。

[参考文献]

[1]A. Dragulescu and V. M. Yakovenko, "Statistical mechanics of money", Eur. Phys. J. B 17, 723 (2000). ( arXiv: cond-mat/0001432 )

[2]A. Dragulescu, V. M. Yakovenko, "Evidence for the exponential distribution of income in the USA", Eur. Phys. J. B 20, 585 (2001). ( arXiv: cond-mat/0008305 )

2025-04-13

誰が量子ゆらぎを見たか ~ Casimir 効果 ~

物理学

カシミール効果というのは、真空中に平板を二枚置いたときに量子ゆらぎの存在によってその平板の間に引力が発生するという現象である。逆に言うとこの現象が存在すると量子揺らぎによるエネルギーが存在することが分かる。

ちなみに、なぜこれに触れるのかというと、岡田勘三著「カシミール物理への招待」 (日本評論社) という本をちらちらとみていたからである。この本自体はいろいろなことを解説していて勉強にはなったが、肝心のカシミール力の導出をしておらず、参考文献に挙げられていた P. W. Millonni "The Quantum Vacuum : An Introduction to Quantum Electrodynamics" という本が高価でなかなか見ることがなさそうなので自分で導出しようと思ったためである。そこでこの記事を書き始めたのは良いものの、カシミール効果の英語版の Wikipedia の記事を見てみると、ゼータ関数を用いた正則化によってカシミール力を導出しているのを発見してしまった。ただ Wikipedia と同じ方法で導出をしてもしょうがないので、この記事ではゼータ関数による正則化を使った導出ではなく、ガウス関数を用いた正則化をつかってカシミール力を導出することを目指す。自分だけが分かるように数式を書きなぐってもしょうがないのでそれなりに順序立てして計算を書いていきたい。ちなみにこの記事のタイトルは江沢洋著「だれが原子をみたか」（岩波書店）のパクｒ…オマージュである。

量子揺らぎのエネルギー

いきなりカシミール効果に行く前にまずは調和振動子について考える。調和振動子のハミルトニアン $H = \frac{p^2}{2m} + \frac{m\omega^2}{2}x^2$ を生成消滅演算子 $a = \sqrt{\frac{m \omega}{2\hbar}}x + \frac{{\rm i}}{\sqrt{\hbar m\omega}}p$ , $a^\dagger = \sqrt{\frac{m \omega}{2\hbar}}x - \frac{{\rm i}}{\sqrt{\hbar m\omega}}p$ を用いて書き直す。ここで交換関係 $[a, a^\dagger] = 1$ に気をつけると $H = \hbar \omega \left(a^\dagger a + \frac{1}{2}\right)$ となる。このハミルトニアンの中に現れる $a^\dagger a$ は粒子数演算子になっていて、粒子の数が0である状態 (真空) が基底状態になる。この基底状態を $|0\rangle$ と書くことにすると $a|0\rangle = 0$ になる。基底状態エネルギー $E_0$ は $H|0\rangle = \frac{\hbar \omega}{2}|0\rangle$ より $E_0=\frac{\hbar \omega}{2}$ であることがわかる。この $E_0$ が量子ゆらぎによるエネルギーであり、零点エネルギーとも呼ばれる。

つぎに場の理論の話に一歩踏み入れる。以下では計算の簡単さのために自然単位系 ( $c=1 \ , \ \hbar = 1$ ) を採用している。場の理論での質量ゼロのスカラー場のハミルトニアンは波数表示で次のようになる。

$\displaystyle H = \int \frac{{\rm d}^3 k}{(2\pi)^3} \omega_{\bf k} (a^\dagger_{\bf k} a_{\bf k} + \frac{1}{2})$

ここで $\omega_{k} = |{\bf k}|$ である。

このとき真空状態 $|0 \rangle$ が任意の ${\bf k}$ について $a_{\bf k}|0 \rangle = 0$ を満たすことに注意すると真空エネルギーは

$\displaystyle E_0 = \langle 0| H|0\rangle = \int \frac{{\rm d}^3 k}{(2\pi)^3} \frac{\omega_{\bf k}}{2} = \infty$

となって無限大に発散することになる。これは非物理的な観測できない量になっていて量子ゆらぎの真空エネルギーが実在するのか分からない。以上で量子ゆらぎについての考察は終わり！めでたしめでたし！！としたい所だがそうは問屋が卸さない。

カシミール効果

1948年に理論物理学者のカシミールという人が真空揺らぎによって金属板の間に力が生じるということを示した。これをカシミール効果という。なぜ引力ができるのかというと、金属板の間の距離より長い波長の振動のモードの揺らぎはこの金属板間の空間には存在できず、その差によって引力が生じる。このカシミール効果が実験で観測されると量子揺らぎの零点エネルギーが存在することを実証できることになる。

では以下でカシミール効果を見ていきたい。設定として、距離 $L$ だけ離れた二枚の金属板1, 3があるとする。この間にもう一枚の金属板2をいれる(図を参照）。金属板1, 3で囲むことでこの真ん中の平板2が動く時に平板の外側の真空エネルギーの影響を考えずに済む。最終的には $L \rightarrow \infty$ とすることで金属板1, 2 の間の力を計算する。簡単のためにこれらの平板は $x$ 軸方向に並んでいるとする。

先ほど見たように、波数の大きな項まで計算するために真空エネルギーというのは発散することになる。こういった発散を抑えるためにガウス正則化というものを用いる。こういった正則化の一番分かりやすいものはハードカットオフだろう。ハードカットオフは和や積分を有限の大きさで切って計算して最後に無限大に飛ばしたりするやり方だ。また最初の方でも触れたが英語版の wiki ではゼータ関数を用いた正則化によって導出がされている。ちなみに、ゼータ関数の解析接続で出てくる $1+2+3+\cdots = -\frac{1}{12}$ という式は一次元のカシミール効果のゼータ正則化で出てくるためカシミール効果と関係があるという説明がよくされている。

この記事では高波数成分を滑らかに抑え、物理的に自然なカットオフを提供するガウシアン正則化を考える。ガウシアン正則化を用いると真空エネルギーは次のようになる。

$\displaystyle \langle E_0(r) \rangle =A \sum_{n=1}^\infty \int \frac{{\rm d}^2 k}{(2\pi)^2} \frac{\omega_k}{2}e^{-(\frac{\omega_k}{\pi \Lambda})^2}$

ここで、距離 $r$ の関数であることを明示するため $E_0(r)$ とした。 $A$ は金属板の面積である。さて、 $x$ 軸方向の波数に関しては境界条件から離散化されているので $\omega_k = \sqrt{(\frac{\pi n}{r})^2 + k_2^2 + k_3^2}$ である。この積分を極座標表示 ( $k_2^2 + k_3^2 = k^2 \ , \ \int {\rm d}^2 k = \int_{0}^{2\pi}{\rm d}\theta\int_0^\infty k{\rm d}k$ ) にして $t^2 = \frac{k^2 + (\frac{\pi n}{r})^2}{\pi^2\Lambda^2}$ と変数変換すると

$\displaystyle \langle E_0(r) \rangle / A =\frac{\Lambda^2 \pi}{4}\sum_{n=1}^\infty \int_{\frac{n}{r \Lambda}}^\infty {\rm d} t \ t^2e^{-t^2}$

という式を得る。ちなみに、エネルギーの両辺を面積 $A$ で割って左辺を単位面積あたりのエネルギーにしている。ここで積分を実行すると

$\displaystyle \langle E_0(r) \rangle / A = \frac{1}{2}\sum_{n=1}^\infty \left\{ \frac{(\pi \Lambda)^2}{4\pi} \left(\frac{\pi n}{r}\right) e^{-\frac{n^2}{(r \Lambda)^2}}+ \frac{(\pi \Lambda)^2}{4\pi} \frac{\Lambda \pi^{3/2}}{2} \left(1 - {\rm erf}\left(\frac{ n}{\Lambda r}\right)\right) \right \}$

となる。ここで和の計算のために次のオイラー・マクローリンの公式を用いる。

$\displaystyle \sum_{j = 0}^{n-1} f(j) = \int_0^n f(x) {\rm d}x + \sum_{k = 1}^m\frac{b_k}{k!}\left(f^{(k-1)}(n)-f^{(k-1)}(0)\right) + R_m$

ここで $b_k$ はベルヌーイ数で、 $R_m$ は補正項でベルヌーイ多項式というものを用いて計算される。ベルヌーイ数の値は $b_1 = -\frac{1}{2}, b_2 = \frac{1}{6}, b_3 = 0, b_4 = -\frac{1}{30}, b_5 = 0, b_6 = \frac{1}{42}, b_7 = 0, \dots$ である。さて、この系では

$\displaystyle f(x) = \frac{(\pi \Lambda)^2}{4\pi} \left(\frac{\pi x}{r}\right) e^{-\frac{x^2}{(r \Lambda)^2}} + \frac{(\pi \Lambda)^2}{4\pi} \frac{\Lambda \pi^{3/2}}{2} \left(1 - {\rm erf}\left(\frac{ x}{\Lambda r}\right)\right)$

であり、その導関数は

$\displaystyle f^{(1)}(x) = -\frac{\pi^2 x^2}{2r^3}e^{-\left(\frac{x}{\Lambda r}\right)^2}$ ,

$\displaystyle f^{(2)}(x) = \frac{\pi^2 x}{r^3}\left(\frac{x^2}{\Lambda^2 r^2}-1\right)e^{-\left(\frac{x}{\Lambda r}\right)^2}$ ,

$\displaystyle f^{(3)}(x) = \frac{\pi^2 }{r^3}\left(-1 + 5\frac{x^2}{\Lambda^2 r^2} -2 \frac{x^4}{\Lambda^4 r^4}\right)e^{-\left(\frac{x}{\Lambda r}\right)^2}$

であるのでこれらをオイラー・マクローリン公式に適用すると、距離 $r$ の金属板に挟まれた領域の単位面積あたりの真空エネルギーは次のように計算される。

$\displaystyle E_0(r) / A \sim \frac{1}{2}\left(\frac{\pi^2 \Lambda^4 r}{4} \rho + \frac{\Lambda^3 \pi^{5/2}}{8}-\frac{\pi^2}{720r^3}\right)$

ここで $\rho = \int_0^\infty te^{-t^2} + \frac{\pi^{1/2}}{2}\left(1 - {\rm erf}\left(t\right)\right) {\rm d}t$ である。またオイラーマクローリン公式の補正項 $R_m$ は十分小さい。

板1 と板3 で囲まれる空間の単位面積あたりの全真空エネルギーは次のようになる。

$\displaystyle E_0^{\rm tot} / A = (E_0(a) + E_0(L-a)) / A \sim \frac{1}{2}\left(\frac{\pi^2 \Lambda^4 L}{4} \rho + \frac{\Lambda^3 \pi^{5/2}}{4}-\frac{\pi^2}{720a^3}-\frac{\pi^2}{720(L-a)^3}\right)$

さて $F = -\frac{{\rm d}}{{\rm d}a} E_0^{\rm tot}$ を計算して $L \rightarrow \infty$ の極限をとると板1 と板2 の単位面積あたりのカシミール力を計算でき、その値は

$\displaystyle F / A= -\frac{\pi^2}{480 a^4}$

となる。以上の計算はスカラー場についての計算であったが電磁場の場合は偏極の自由度が 2 なので値が倍になって

$\displaystyle F / A= -\frac{\pi^2}{240 a^4}$

を得る。ここで負符号がついているのは引力であることを表している。

カシミール効果の実験

このカシミール効果の実験について触れておきたい。まず、1958年にSparnaay によって金属平板間に働く力を測定する実験が行われた[1]。この実験の結果、クロム-クロム間, クロム-スチール間では引力が生じたものの、アルミ-アルミ間では斥力が生じることが観測された。この斥力は洗浄しても落ちなかった埃塵や静電ポテンシャルに差があったことによるものと考えられている（カシミール力は小さいので静電気力が働けば無視できる）。この実験によってカシミール力の測定には静電ポテンシャルの差・不純物の除去・平板間の距離の精密な測定など実験するにあたって何が必要かという解像度が高くなった。

1978年には van Blokland と Overeek が同じく金属平板間でのカシミール効果を測定した[2]。この実験では Sparnaay による実験の問題点を改良したものだった。しかし実験の精度は50％程度と見積もられていてあまり精度は高くなかったようだ。金属板同士をかなり高い精度で平行に保つというのが難しいというのも理由の一つにある。

その後、1997年に Lamoreaux によってより高い精度の実験が行われた[3]。この実験は高い精度で平行を保つ必要のある金属板同士の間に働くカシミール力ではなく、金でコーティング ( 蒸着 ) された球面と金でコーティングされた平板との間のカシミール力を測定した。

ここでこの実験のように球と平板の間のカシミール力を導出したい。球の半径を $R$ とし、球面と平板の間の距離を $a$ とすると(下図を参照)球面上の点 $(x, y, z)$ の満たす式は $x^2 + y^2 + (z-a-R)^2 = R^2$ となり、 $z = a + R -\sqrt{R^2-x^2-y^2}$ を得る。

ここで $x^2 + y^2 \ll R^2$ という近似をすると $z = a + \frac{x^2+y^2}{2R}$ となる。この近似は近接力近似とよばれ $a \ll R$ で有効である。単位面積あたりの電磁場の真空エネルギーで距離の関数の部分は $e_0 = -\frac{\pi^2}{720r^3}$ であったので ( スカラー場の2倍 )、 $r = a + \frac{x^2+y^2}{2R}$ として $x$ , $y$ について積分すると真空エネルギーが得られる。

$\displaystyle E_0(a) =-\frac{\pi^2}{720}\int \int {\rm d}x {\rm d} y\frac{1}{\left(a + \frac{x^2 + y^2}{2R}\right)^3} =-\frac{\pi^2}{720}\int_0^\infty 2\pi \rho {\rm d} \rho \frac{1}{\left(a + \frac{\rho^2}{2R}\right)^3} =-\frac{\pi^3R}{720a^2}$

この結果を $a$ について微分すると球と金属板の間のカシミール力を得られる。

$\displaystyle F =-\frac{{\rm d}}{{\rm d} a}E_0(a) = -\frac{\pi^3R}{360a^3}$

実際の実験に使う式では金属の誘電率の周波数依存性を考慮した補正項がかかっている。

この実験では球面はピエゾ素子に取り付けられ、金属板との距離が精密に制御されている。ピエゾ素子というのは圧力をかけて変形して歪みが生じると電気的な応答(分極) をして、逆に電気分極によって変形する物質のことである。ピエゾ素子は走査型トンネル顕微鏡などにも応用されている。平板の方はねじれ振り子に取り付けられている。カシミール力の引力によって平板は球面にひきつけられ何もなければねじれ振り子がねじれるが、実験装置のフィードバック制御によって振り子の角度は一定に保たれる。角度を定に保つ電圧の大きさからカシミール力を計算して測定している。この Lamoreaux の実験によって5％の精度でカシミール効果が実証された。すなわち零点エネルギーが存在することが確かめられた。これらの実験によってカシミール効果が実証され、量子揺らぎが存在することが証明された。この後もより高い精度での実験が数多く行われた。

以上をまとめると次の表のようになる。

実験者	年	精度	対象
Sparnaay	1958	不明	平板-平板
van Blokland	1978	~50%	平板-平板
Lamoreaux	1997	~5%	球面-平板

まとめ

カシミール効果というのは物理学者のカシミールによって理論的に予言された現象で、この記事ではスカラー場の真空エネルギーの式からカシミール力をガウス正則化を用いて $F / A= -\frac{\pi^2}{480 a^4}$ と導出した。このカシミール効果は実験でも検証されていて、1997年に球面と平板との間の引力を測定することで理論と高い精度で一致することが示されるを解説した。この実験から量子ゆらぎが実在することが示された。

誘電体で満たされている場合やフェルミオン場ではカシミール力が斥力になったり、有限温度の場合にどうなるかや、平板の振動による動的カシミール効果など面白い現象はまだまだあり、ナノマシンなどの微小な領域での物理を考えるとカシミール効果が効いてくるので応用面でもカシミール効果がますます重要になってくると予想される[4]ので機会があればそれらについても触れてみたい。

[参考文献]

[1] M. J. Sparnaay, Phyisica 24, 751.

[2] P.H.G.M.van Blokland, J.T.G.Overbeek, J.Chem.Soc.Faraday Trans. 74, 2637.

[3] S. K. Lamoreaux, Phys. Rev. Lett. 78, 5.

[4] 岡田勘三, カシミール物理への招待, 日本評論社, 2020.

2025-02-01

高校数学でわかるリボ払い 2

数学

前回のブログでは定額支払い方式でのリボ払いを扱った。

tori-head.hatenablog.com

今回はリボ払いの定率返済方式と残高スライド方式について考えてみたいと思う。

まず、定率返済方式に関して、三井住友カードによると「定率返済方式は利用残高に対して一定割合をかけた金額を支払う方式」と書いてある。

ここで、借入残高を $S_m$ 、月利を $i$ 、毎月の返済率を $q$ とする。毎月の支払額が、その時点の利用残高に返済率 $q$ を掛けた $qS_m$ であるとすると、翌月の残高は

$S_{m+1}=(1+i)S_m-qS_m=(1+i-q)S_m$

となる。したがって、初期の借入残高を $S_0$ とすれば、 $m$ ヵ月後の残高は

$S_m=(1+i-q)^mS_0$

となる。残高が減少するためには

$q>i$

である必要がある。この条件を満たしていれば $0<1+i-q<1$ となり、残高は指数的に減衰する。

なお、利息を毎月の返済額とは別に支払い、元金についてのみ利用残高の一定割合 $q$ を返済する方式であれば、残高の漸化式は

$S_{m+1}=(1-q)S_m$

となり、

$S_m=(1-q)^mS_0$

となる。この場合でも、指数関数的な減衰である以上、数学的には有限回の返済で残高が厳密に 0 になることはない。

実際の返済では、残高が十分に小さくなれば最終回に残額をまとめて支払う、あるいは最低支払額が設定されるなどの処理が行われると考えられる。

例えば最低支払額を $p$ とし、通常の定率返済額 $qS_k$ が $p$ を下回った時点から毎月 $p$ を支払うと仮定する。この場合、利息を考慮するなら、それ以降の残高は単純に $S_k-np$ とはならず、

$S_{k+n+1}=(1+i)S_{k+n}-p$

という元利定額返済と同じ形の漸化式に従う。

したがって、定率返済方式を数式化する際には、金利 $i$ と返済率 $q$ を区別し、さらに最低支払額や最終返済の処理を別途考える必要がある。

さて、次は残高スライド方式について考える。三井住友カードによると、残高スライド方式は「借入残高に応じて月々の返済額が変動する返済方式のこと」とあり、実用上は残高スライド元利定額返済方式があるようだ。残高スライド元利定額返済方式は図1のように借入金の残高に応じて月々の支払い金額が変わるシステムである。前回の記事で少し触れたが元利定額支払いの場合には $rS_0/p \ge 1$ を満たすと残高が減らず支払い期限が無限に発散することになるのを防ぐシステムになっている。一方で残高が減るにつれて月々の支払い額が減るので元利定額方式よりも支払い期間が長くなる。

さて、残高スライド元利定額払いを数式にのせるために図１のように設定する。借入金額 $S_0$ は $\sigma_i \lt S_0 \le \sigma_{i+1}$ になっており、 $m$ ヵ月後の残高 $S_m$ が $\sigma_{j-1} \lt S_0 \le \sigma_{j}$ を満たす時の月々の支払い額が $p_j$ であるとする。また、 $0 \lt \sigma_1 \lt \cdots \lt \sigma_j \lt \cdots \lt \sigma_i \lt S_0$ かつ $p_1 \lt p_2 \lt \cdots \lt p_i \lt p_{i + 1}$ となっている。

さて、このとき $m_{i + 1}$ ヵ月で残高 $S_{m_{i+1}}$ が $S_{m_{i+1}} \lt \sigma_i$ をみたし、そのひと月前の残高 $S_{m_{i+1}-1}$ は $S_{m_{i+1}-1} \ge \sigma_i$ となるとする。このとき、前記事の結果を用いると

$\displaystyle S_{m_{i+1}} = (1+r)^{m_{i+1}}S_0-\frac{p_{i+1}}{r}\left((1+r)^{m_{i+1}}-1\right)$

となる。そこからさらに $m_i$ ヵ月たち、残高 $S_{m_{i+1} + m_i}$ は $S_{m_{i+1} + m_i} \lt \sigma_{i-1}$ かつ $S_{m_{i+1} + m_i-1} \ge \sigma_{i-1}$ をみたすとする。この時の $S_{m_{i+1} + m_i}$ の表式は $S_{m_{i+1}}$ の式で $S_0 \rightarrow S_{m_{i+1}}$ 、 $m_{i+1} \rightarrow m_i$ 、 $p_{i+1} \rightarrow p_i$ 　としたものなので

$S_{m_{i+1}+m_i} = (1+r)^{m_i}S_{m_{i+1}}-\frac{p_i}{r}\left((1+r)^{m_i}-1\right)$

となり、 $S_{m_{i+1}}$ に先ほどの式を代入すると

$\displaystyle S_{m_{i+1}+m_i} = (1+r)^{m_{i+1}+m_i}S_0-\frac{p_{i+1}}{r}\left((1+r)^{m_{i+1}+m_i}-(1+r)^{m_i}\right)-\frac{p_i}{r}\left((1+r)^{m_i}-1\right)$

となる。これを繰り返していくと次の式を得る。

$\displaystyle S = (1+r)^{\sum_j m_j}S_0-\sum_j\frac{p_j}{r}\left\{(1+r)^{\sum_{k=0}^{j} m_k}-(1+r)^{\sum_{k=0}^{j-1} m_k}\right\}$

これが厳密な式なので $S = 0$ として方程式をとくと支払い期間などが分かるはずだが人間にはおそらく解くことは難しいと思われる。

厳密な結果は解けないので次のような近似を考えることにする。 $S_0$ から $m_{i + 1}$ ヵ月で残高が $\sigma_i$ になり（ $S_{m_{i+1}} = \sigma_i$ ）、次に $m_{i}$ ヵ月で残高が $\sigma_i$ から $\sigma_{i-1}$ に減り、それを繰り返すことで最終的に借入残高が 0 になる。

最初の $m_{i+1}$ ヵ月で $S_{m_{i+1}} = \sigma_i$ より $S_{m_{i+1}} = (1+r)^{m_{i+1}}S_0-\frac{p_{i+1}}{r}\left((1+r)^{m_{i+1}}-1\right) = \sigma_i$ となる。これを $m_{i+1}$ についてとくと、

$\displaystyle m_{i+1} = \frac{ \log\left(\frac{p_{i+1}}{r}-\sigma_i\right) - \log\left(\frac{p_{i+1}}{r}-S_0\right)}{\log(1+r)}$

を得る。ここで支払い期間が無限大に発散しないために $\frac{p_{i+1}}{r} \gt S_0$ となっているので $\log$ の中身が上記のようになっている。さて、次に、残高が $\sigma_j$ から $\sigma_{j-1}$ へと下がっていくにとき（ただし $j = 1, \cdots, i$ で $\sigma_0 := 0$ と定義する）、月々の支払い金額は $p_j$ なので先ほどと同様にして支払い期間は

$\displaystyle m_j = \frac{ \log\left(\frac{p_j}{r}-\sigma_{j-1}\right) - \log\left(\frac{p_j}{r}-\sigma_j\right)}{\log(1+r)}$

となる。ここでも支払い期間が無限大に発散しないために $\frac{p_j}{r} \gt \sigma_j$ となっているので $\log$ の中身が上記のようになっている。これらを足し合わせた $\sum_{k=1}^{i+1} m_k$ がトータルの支払い期間になる。また、トータルの支払い金額 $T$ は $T = \sum_{k=1}^{i+1} p_km_k$ で与えられる。すべての $j$ について $p_j = p$ が成り立つとき、足し合わせた項同士が打ち消しあって支払い期間と支払い総額が前回の記事でみた元利定額支払いと一致することが分かる。

これまでの結果を応用するために三井住友銀行の例で出てる残高スライド元利定額返済方式の例を使って計算してみたい。借入金を240 万円 $S_0 = 240$ , 年利 15％の利子率つまり月利は 1.25% $r = 15/1200 = 0.0125$ であるとする。

この時、 $m_4 = \frac{ \log\left(\frac{p_4}{r}-\sigma_3\right) - \log\left(\frac{p_4}{r}-S_0\right)}{\log(1+r)}$ より $m_4 \fallingdotseq 32.64$ となる。また、 $m_j = \frac{ \log\left(\frac{p_j}{r}-\sigma_{j-1}\right) - \log\left(\frac{p_j}{r}-\sigma_j\right)}{\log(1+r)}$ $(j = 1, 2, 3)$ より $m_3 \fallingdotseq 100.84$ 、 $m_2 \fallingdotseq 48.79$ 、 $m_1 \fallingdotseq 78.95$ となる。したがって全期間は $\sum_{k=1}^4 m_k \fallingdotseq 261.22$ ヵ月となる。これは 21年10ヵ月もの期間になる。また、このときの支払い総額 $T$ は $T = \sum_{k = 1}^4 p_k m_k \fallingdotseq 609.61$ 万円となる。これではもはや元々の借入金よりも支払った利息の方が大きくなってしまっている。やはりここら辺がリボ払いが悪名高い所以だろう。

数値計算をすると支払い期間は 259 ヵ月で、支払い総額は 606.26 万円程度となる。上の支払い期間と支払い総額を比較したときにそれぞれの誤差は 0.86% と 0.55% となっていて上の計算は良い近似になっていることがわかる。下記に実際に使ったC++のコードを載せておく（お恥ずかしながら綺麗な良いコードでないので申し訳ありません）。

ちなみに、上記の借入金 240万円と月利 1.25% で毎月の支払いが 4万円での元利定額支払いであれば支払い期間は 112ヵ月程度(9年4ヵ月) で支払い総額は 446.38 万円程度になる。つまり残高スライド式にすることで支払い期間が元利定額方式より2倍以上で150万近く余分に払っていることになる。

以上でリボ払いについての考察を終えたいと思う。前回と今回で扱った支払い方式の他にも支払い方法がいくつかあるのでそれを数式に乗せたり、残高スライド元利定額方式で支払い総額や支払い期間を最大化するような $p_j$ 、 $\sigma_j$ の組み合わせは何かなど興味は尽きないがキリが良いのでここまでとしたい（また何か思いついたら記事を書くかもしれない）。

``` C++
#include <iostream>

using namespace std ;

int main(){

double balance = 240 ;

double r = 0.0125 ;

double plist[4] = {1, 2, 3, 4} ;

double sigmalist[4] = {50, 100 ,200, 300};

double total_pay = 0 ;

int total_month = 0 ;

while(balance >= 0){

if(balance < sigmalist[3] && balance >= sigmalist[2]){

balance = (1+r) * balance - plist[3];

total_pay += plist[3] ;

total_month++ ;

} else if(balance < sigmalist[2] && balance >= sigmalist[1]){

balance = (1+r) * balance - plist[2];

total_pay += plist[2] ;

total_month++ ;

} else if(balance < sigmalist[1] && balance >= sigmalist[0]){

balance = (1+r) * balance - plist[1];

total_pay += plist[1] ;

total_month++ ;

} else {

if((1+r) * balance - plist[0] < 0 ){

total_pay += balance ;

total_month++ ;

break ;

}

balance = (1+r) * balance - plist[0];

total_pay += plist[0] ;

total_month++ ;

}

}

cout << "total month = " << total_month << endl ;

cout << "total payment = " << total_pay << endl ;

return 0 ;

}
```

参考文献：

www.smbc-card.com

2025-01-29

高校数学でわかるリボ払い

数学

世に悪名を轟かせているものには世襲政治家や官僚やワンパンパスタなどがある。とくに最後のものは不味いし時間はかかるのに持て囃されている意味不明なシロモノである。

さて、他にも世に悪名を轟かせているモノがある。それがリボ払いである。実際、「リボ払い」で検索しようとすると「リボ払いヤバイ」が候補に上がってくるぐらい悪名を轟かせている。あまりにも印象が悪いものだからカード会社はリボ払いという単語をごまかして一見するとリボ払いだと分からないようにしているようだ。しかしながら、あまりにもヤバいという情報が先走っていてリボ払いの中身がよくわからない。リボ払いのシミュレーターを使えばどれくらいヤバいかは分かるかもしれないが、どういう計算をしてるのか中身がよくわからないのでこの記事で見ていきたいと思う。具体的にはリボ払いでの支払い総額を計算して元の金額との差額を計算してどの程度違うのかを見てみたい。

三井住友カードによるとリボ払いには大まかには定額返済方式と残額スライド方式がある。この記事では定額返済方式のみを扱うことにする（計算が楽だから）。

さらに定額返済方式にも元金定額方式と元利定額方式との二種類がある。

まずは元金定額方式について見てみる。三井住友銀行の解説によると「元金定額方式は、元金を一定額に設定して、そこに手数料を計算して上乗せする方式」とのことである。

ここで、これを数式にするために以下のように設定する。元金の金額を $p$ , 利率を $r$ , 借入から m ヵ月後の残高を $S_m$ とする。このとき手数料は残高と利率の積 $r S_m$ で与えられる。 $m-1$ ヵ月後の残高から $p$ だけ引いた分が翌月の残高になるので $S_m$ に関する漸化式は

$\displaystyle S_m = S_{m-1}-p$

となってこれをとくと $S_m = S_0 - mp$ となる。ここで $S_0$ は最初に借りたお金あるいはカードで買った金額である。簡単のために $N = S_0/p$ とおいて、 $S_0$ は $p$ で割り切れるとする ( $N \in {\mathbb N}$ ) 。 $m$ ヵ月後の支払い金額は元金と手数料の和 $p+rS_m$ で与えられるのでトータルの支払い金額 $T$ は

$\displaystyle T = \sum_{m=0}^{N-1}(p+rS_m)=\sum_{m=0}^{N-1}(p + r\left(S_0-pm\right))$

これを計算すると $T = N(p + rS_0)-rp\frac{(N-1)N}{2}$ となる。 $N = S_0/p$ 　であったので整理すると　 $T = \left(\frac{S_0}{p}\right)\left(\left(1+\frac{r}{2}\right)p+\frac{rS_0}{2}\right)$ となる。支払い額と借りたお金との差は $T - S_0 = \frac{rS_0}{2} + \frac{rS_0^2}{2p}$ となる。

以上の計算を応用してみることを考える。仮に年率 15％のところで 20万円借りて月々の元金が 5千円の元金定額方式で返済をするとすると、 $r$ はひと月の利率なので $r = \frac{15}{12 \times 100}=0.0125$ 、 $p =0.5$ 、 $S_0 = 20$ から $T = 25.125$ となるので、リボ払いになると5万1250円ほど利息として支払っていることになる。この時、 $N = 20/0.5 = 40$ なので支払い期間は 3年4ヵ月程度になる。

次に元利定額方式について見る。三井住友カードの解説によると「元利定額方式は、元金に手数料を含めて計算して、毎月一定額を返済していく方式」とのことである。

先ほどと同様に元金の金額を $p$ , 利率を $r$ , 借入から mヵ月後の残高を $S_m$ とする。このとき、手数料は残高と利率の積 $r S_m$ で与えられる。元金と手数料の差 $\Delta_m = (p-rS_m)$ によって残高が減ることになるので翌月の残高は $S_m - \Delta_m$ になる。これから $S_m$ に関する漸化式を以下のようにつくることができる。

$\displaystyle S_m = S_{m-1}-\Delta_{m-1} = S_{m-1}-(p-r S_{m-1})$

この式を整理すると

$\displaystyle S_m = (1+r)S_{m-1}-p$

となる。この式に $m$ を1だけずらした漸化式 $S_{m-1}=(1+r)S_{m-2}-p$ を代入すると

$\displaystyle　S_m = (1+r)^2S_{m-2}-p(1+(1+r))$

を得る。これを繰り返すことで最終的に

$\displaystyle　S_m = (1+r)^mS_{0}-p(1+(1+r)+\cdots +(1+r)^{m-1})=(1+r)^mS_{0}-p\frac{(1+r)^m-1}{r}$

を得ることができる。（もっとエレガントには $S_m - p/r = (1+r)(S_{m-1}-p/r)$ と変形して一般項を求めることができる。）ここで $S_0$ は最初にリボ払いで借りた金額である。簡単のために支払いが $n$ ヵ月後に終わるつまり $S_n = 0$ となるとすると、

$\displaystyle　(1+r)^nS_{0}=p\frac{(1+r)^n-1}{r}$

が成り立つ。これを $n$ についてとくと

$\displaystyle　n=-\frac{\log\left(1-\frac{rS_0}{p}\right)}{\log(1+r)}$

となる。ここで $\frac{rS_0}{p} \ge 1$ の時 $n$ が発散することがわかる。すなわち、借金が払いきれず支払いが完了できなくなることが分かる。トータルの支払い金額 $T$ は $T = np=-p\frac{\log\left(1-\frac{rS_0}{p}\right)}{\log(1+r)}$ となる。

先ほどと同様に年率 15％のところで 20万円借りて月々 5千円の返済をするとすると、 $r$ はひと月の利率なので $r = \frac{15}{12 \times 100}=0.0125$ 、 $p =0.5$ 、 $S_0 = 20$ から $T = 27.899$ となるので、リボ払いになると7万8990円ほど利息として支払っていることになる。また、 $n = 55.7976\dots$ なので支払い期間は 56ヵ月 (4年 8ヵ月) 程度になることがわかる。

こちらの金融経済教育指導教材の5章にある問題では16％になっているのでその数値でも計算してみると良いかもしれない。

www.fsa.go.jp

この記事では元金定額方式と元利定額方式での支払い金額の総額を計算したが残金スライド方式についても余裕があれば考えたいと思う。また、元利定額方式の漸化式は複利のローンの残高の漸化式と同じなので利子率などを変えるとローンの理解にも役立つと思う。とくに今回は支払い総額に注目して $p$ , $S_0$ , $r$ に具体的な値を入れて $n$ を決めたが、例えば $n$ , $p$ , $r$ に値を入れて許容できる最大の借入額 $S_0$ などを求めてみて遊んでも良いと思う。(他の値から利子率 $r$ を決める計算は数値的にやらなければならなかったと思うが)

追記：

こちらの記事で残高スライド元利定額支払い方式をメインに扱った記事を書いたので興味があればこちらもぜひ読んでみてください。

tori-head.hatenablog.com

参考文献 :