カテゴリー: 数学

単振動の微分方程式の解

/ tau / コメントする

概要

以下の単振動に関する方程式の解法を整理した。

(1)   \begin{equation*} \frac{d^2x(t)}{dt^2} + \omega^2 x(t) = 0 \quad \left( t = 0 \Rightarrow \left\{ \begin{array}{l} x = X_0 \\ \dot{x} = 0 \end{array} \right) \end{equation*}

変数変換による方法

式(1)の両辺にdx/dtを掛けて変形する。

(2)   \begin{align*} &\frac{dx}{dt} \frac{d^2x}{dt^2} + \omega^2 x \frac{dx}{dt} = 0 \\ &\frac{1}{2} \frac{d}{dt} \left( \frac{dx}{dt} \right)^2 + \frac{1}{2} \omega^2 \frac{d}{dt} \left( x^2 \right) = 0 \end{align*}

これをtで積分して以下を得る。

(3)   \begin{equation*} \left( \frac{dx}{dt} \right)^2 + \omega^2 x^2 = C \end{equation*}

ここでC = ω2A2とおいて

(4)   \begin{equation*} \frac{dx}{dt} \right = \pm \omega \sqrt {A^2 - x^2} \quad (x \le A) \end{equation*}

さらにx = A cos θとおいて、

(5)   \begin{equation*} (- A \sin \theta) \frac{d \theta}{dt} = \pm \omega A \sin \theta \quad \Rightarrow \quad \frac{d \theta}{dt} = \mp \omega A \end{equation*}

θの解は以下の様になる。

(6)   \begin{equation*} \theta = \mp \omega t + \alpha \end{equation*}

x = A cos θだったので、以下を得る。

(7)   \begin{equation*} x = A \cos (\omega t + \alpha) \end{equation*}

初期条件から

(8)   \begin{align*} x|_{t=0} &= A \cos (\alpha) = X0 \\ \dot{x}|_{t=0} &= -A \sin (\alpha) = 0 \end{align*}

これより以下の解を得る。

(9)   \begin{equation*} x = X_0 \cos \omega t \end{equation*}

オイラーの公式を使う方法

式(1)の解を以下の様に置く。

(10)   \begin{equation*} x = e^{\lambda t} \end{equation*}

これを式(1)に代入して

(11)   \begin{equation*} \lambda^2 e^{\lambda t} + \omega^2 e^{\lambda t} = 0 \quad \Rightarrow \quad \lambda ^2 + \omega^2 = 0 \quad \Rightarrow \quad \lambda = \pm i \omega \end{equation*}

この解を式(10)に代入して2つの関数が得られるが、式(1)の一般解はこれらの線形結合として得られる。

(12)   \begin{equation*} x = A e^{i \omega t} + Be^{-i \omega t} \end{equation*}

ここで以下のオイラーの公式を適用する。

(13)   \begin{equation*} e^{\pm i \theta} = \cos \theta \pm i \sin \theta \end{equation*}

これにより式(12)は以下のように書ける。

(14)   \begin{align*} x &= C_1 (\cos \omega t + i \sin \omega t) + C_2 (\cos \omega t - i \sin \omega t) \\ &= (C_1 + C_2) \cos \omega t + (C_1 - C_2) i \sin \omega t \end{align*}

xが実数解を持つために、cos, sin双方の係数が実数となるよう、複素係数C1, C2を以下の様に置く。

(15)   \begin{align*} C_1 = \frac{A_1 + i A_2}{2} \\ C_2 = \frac{A_1 - i A_2}{2} \end{align*}

これにより、式(16)は以下の様に変形される。

(16)   \begin{equation*} x = A_1 \cos \omega t + A_2 \sin \omega t \end{equation*}

三角関数の和の関係から、上式は以下の様に表せる。

(17)   \begin{equation*} x = A \cos (\omega t + \alpha) \end{equation*}

初期条件(8)より式(9)と同じ解を得る。

 

 

Runge-Kutta法による1階常微分方程式の解

/ tau / コメントする

概要

1階常微分方程式について、4次のRunge-Kutta法(RK4)による数値解と解析解を比較する。Euler法による解についても重ねてみた。

  • RK4は、かなり粗い分割間隔でもよい結果になる
  • Euler法は、分割間隔を小さくとらないと精度が上がらない

Logistic関数

方程式と解析解

Logistic関数に関する微分方程式は以下で表される。

(1)   \begin{equation*} \frac{N(t)}{dt} = r \left( 1 - \dfrac{N(t)}{K} \right) N(t) \end{equation*}

この方程式は解析的に解けて、以下の解を得る。

(2)   \begin{equation*} N(t) = \dfrac{K}{1 + \left( \dfrac{K}{N_0} - 1 \right) e^{-rt}} \end{equation*}

数値解

式(1)をRK4、Euler法によって解き、解析解と重ねる（解き方は1階常微分方程式の数値解法を参照）。

r=10, K=1000として計算した結果は以下の通り。

  • RK4は上限に収束するまでの時間(t=1.0)の20等分くらい(dt=0.05)でもよく合っている
  • Euler法では凹な部分で計算値が過少となり、凸の部分で解析解に追いついているが、適合度はよろしくない。

dtを上記の10分の1(dt=0.05)として計算した結果が以下の通り。Euler法もかなり適合度合いがよくなったが、それでもまだ少し右方にずれている。

逆にdt=0.1と大きくしてみた。折れ線がはっきりわかるくらいの粗さだが、RK4はよくフィットしている。Euler法のずれは大きくなっている。

コード

上記の結果は、以下のコードを実行して得ている。

三角関数

方程式と解析解

簡単な正弦関数を考える。

(3)   \begin{equation*} y = \sin t \end{equation*}

この関数は、以下の微分方程式の解である。

(4)   \begin{equation*} y' = \cos t \quad ( t=0 \Rightarrow y=0 ) \end{equation*}

数値解

式(4)をEuler法、Runge-Kutta法によって解き、式(3)のグラフと重ねる。

以下は分割間隔を1/40周期で計算した結果。RK4の精度はいいが、Euler法では上に凸のときに過大となり、下に凸のときに過少となって初期値に戻っている。

分割幅を1/200周期にするとEuler法でも精度は上がってくるが、まだ少し誤差が残っている。

分割数を10とかなり粗くとっても、RK4の適合度合いはよい。

繰り返し数が増えた時に誤差が累積するかどうか確認してみた。

以下は10万周期分計算した最後の2周忌を取り出した結果で、RK4、Euler法とも傾向は変わらない。規則正しい周期関数の場合は、繰り返し数を増やしても誤差は拡散しないと思われる。

コード

以下はこれまでの計算をしたPythonのコード。

計算区間を長くとって一部だけ表示させるため、Logistic関数のときとは異なる処理をしている。

  • ループの各段階で計算結果をリストに保存せず、各ステップで計算された次のステップの結果を元の変数に入れ替える
  • 予め指定した計算区間の間の結果だけ、表示用のリストに保存する

 

2階常微分方程式の数値解法

/ tau / コメントする

概要

2階常微分方程式の数値解法、Euler法とRunge-Kutta法をPythonのモジュールとしてまとめた。

1階常微分方程式の解法についてはこちら。

基本形

以下のような方程式を考える。

(1)   \begin{equation*} \dfrac{d^2 y}{dx^2} = f\left ( x, y, \dfrac{dy}{dx}, ... \right ) \end{equation*}

(2)   \begin{equation*} y'' = f(x, y, y', ... ) \end{equation*}

y‘(x)をv(x)と置いて以下のように変形する。

(3)   \begin{align*} \frac{dy}{dx} &=y' = v \\ \frac{dv}{dx} &= v' = f(x, y, y', ...) \end{align*}

Euler法

Euler法では、iステップ目のΔv→Δyと計算し、それらを使ってi+1ステップ目のyとvを計算する。

(4)   \begin{align*} \Delta v_i &= \Delata x \cdot f(x_i, y_i, y'_i, ...) \\ \Delta y_i &= \Delta x \cdot v_i \\ y_{i+1} &= y_i + \Delta y_i \\ v_{i+1} &= v_i + \Delta v_i \end{align*}

Runge-Kutta法

Runge-Kutta法も同様に、iステップ目のΔv1→Δy1を計算し、それらを使ってΔvi、Δyi (i = 1, 2, 3, 4)を計算していく。

(5)   \begin{align*} \Delta v_{1, i} &= \Delta x \cdot f(x_i, y_i, v_i, ...) \\ \Delta y_{1, i} &= \Delta x \cdot v_i \\ \Delta v_{2, i} &= \Delta x \cdot f \left ( x_i + \dfrac{\Delta x}{2}, y_i + \dfrac{\Delta y_{1, i}}{2}, v_i + \dfrac{\Delta v_{1, i}}{2}, ... \right ) \\ \Delta y_{2, i} &= \Delta x \cdot \left (v_i + \dfrac{\Delta v_{1, i}}{2} \right ) \\ \Delta v_{3, i} &= \Delta x \cdot f \left ( x_i + \dfrac{\Delta x}{2}, y_i + \dfrac{\Delta y_{2, i}}{2}, v_i + \dfrac{\Delta v_{2, i}}{2}, ... \right ) \\ \Delta y_{3, i} &= \Delta x \cdot \left (v_i + \dfrac{\Delta v_{2, i}}{2} \right ) \\ \Delta v_{4, i} &= \Delta x \cdot (v_i + \Delta v_{3, i}) \\ \Delta y_{4, i} &= \Delta x \cdot f(x_i + \Delta x, y_i + \Delta y_{3, i}, v_i + \Delta v_{3, i}, ...) \\ x_{i+1} &= x_i + \Delta x \\ y_{i+1} &= y_i + \frac{\Delta y_{1, i} + 2 \Delta y_{2, i} + 2 \Delta y_{3, i} + \Delta y_{4, i}}{6} \\ v_{i+1} &= v_i + \frac{\Delta v_{1, i} + 2 \Delta v_{2, i} + 2 \Delta v_{3, i} + \Delta v_{4, i}}{6} \end{align*}

数値解法モジュール

コード

Euler法とRunge-Kutta法で2階常微分方程式を解くモジュールを、Pythonで以下の様に書いた。

関数euler2、rk4_2の使い方はいずれも同じで、以下の通り。

  • x, y, v：現在の計算ステップでの独立変数x、xにおけるyの値、xにおけるyの1次微分値
  • dx：計算の分割幅
  • func：微分方程式の右辺の関数
  • params：関数funcに与えるパラメーター

パラメーターは、呼び出し側の関数定義で設定したパラメーター名をキーとした辞書で与える。その前提で、モジュール内で**paramsとしてアンパックしている。

戻り値は実数型の数値。

使い方

モジュールの使い方の概要は以下のとおり。

  1. 右辺の関数を定義する
  2. 関数に渡すパラメーターを設定する
  3. 計算開始点と終了点のxの値、分割幅dx, yとvの初期値を設定する
  4. x, y, vの計算用リストを確保する
  5. 計算開始点のyとvのリストに初期値をセットする
  6. 計算実行

単純な放物線の微分方程式を例にしてモジュールの使い方を示す。まず、解の方程式として以下を考える。

(6)   \begin{equation*} y = a x^2 \quad (-1 \le x \le 1) \end{equation*}

この式を2階微分した、以下の常微分方程式を解いていく。

(7)   \begin{equation*} y'' = 2a \quad \left( x = -1 \Rightarrow \left\{ \begin{array}{l} y = a\\ y' =-2a \end{array} \right) \end{equation*}

まず、必要なパッケージをインポート。diffeqは上で示したモジュールを収めている。

次に、微分方程式右辺の関数を定義する(以下のderivative)。この関数はパラメーターとして係数aを1つ持っている。また、パラメーターの値を辞書で定義する。

数値計算に計算開始点x_startと終了点x_end、分割幅dx、yとvの計算初期値を設定。

次に、予めxの計算点のリストを用意して、これに対応するy, vの値を計算・保存していく。計算点数は植木算で分割数＋1とし、ループ計算に用いるリストは計算終了点を除いている。結果を保存しないなら、リストではなく逐次各ステップの計算値を次のステップの入力値にする。

また、計算開始点におけるy, vの初期値をセットしている。

計算自体はシンプルで、定義域リストの値を1つずつ進めながら、現在の値や関数への参照、パラメーターを与えてモジュールを呼び出す。

実行例

元の解析解と併せて表示した実行結果は以下の通りで、1回微分方程式の場合と同じく、Euler法の解は計算が進むごとに乖離が大きくなっている。Runge-Kutta法の解はかなり正確に解析解と一致している。

以下は、解析解も含めてグラフ表示するための、全体の実行コード。

 

1階常微分方程式の数値解法

/ tau / コメントする

概要

1階常微分方程式の数値解法、Euler法とRunge-Kutta法をPythonのモジュールとしてまとめた。

2階常微分方程式の解法についてはこちら。

基本形

以下のような方程式を考える。

(1)   \begin{equation*} \dfrac{dy}{dx} = y' = f(x, y, ...) \end{equation*}

Euler法

Euler法はシンプルな方法で、ある計算ステップiの時の関数上の点(xi, yi)とその点における微分係数y‘iを使って、計算幅dxだけ先にある点xi+1 = xi + dxのに対応するyiの値を近似的に計算する。

(2)   \begin{align*} \Delta y_i &= \Delta x \cdot f(x_i, y_i, ...) \\ x_{i+1} &= x_i + \Delta x \\ y_{i+1} &= y_i + \Delta y_i \end{align*}

概念図からもわかるように、たとえば下に凸の関数の場合、微分係数で線形的に予測した値は常に真値より小さくなるため、誤差が拡大していくことが予想される。

Runge-Kutta法

Runge-Kutta法には1次、2次など異なる次数の解法があるが、4次の方法がよく紹介されている。

(3)   \begin{align*} \Delta y_{1, i} &= \Delta x \cdot f(x_i, y_i, ...) \\ \Delta y_{2, i} &= \Delta x \cdot f \left ( x_i + \dfrac{\Delta x}{2}, y_i + \dfrac{\Delta y_{1, i}}{2}, ... \right ) \\ \Delta y_{3, i} &= \Delta x \cdot f \left ( x_i + \dfrac{\Delta x}{2}, y_i + \dfrac{\Delta y_{2, i}}{2}, ... \right ) \\ \Delta y_{4, i} &= \Delta x \cdot f(x_i + \Delta x, y_i + \Delta y_{3, i}, ...) \\ x_{i+1} &= x_i + \Delta x \\ y_{i+1} &= y_i + \frac{\Delta y_{1, i} + 2 \Delta y_{2, i} + 2 \Delta y_{3, i} + \Delta y_{4, i}}{6} \end{align*}

これを説明するのに、x–y平面上での関数を描き、その微分係数を複数描いているものが多い。

本来は、右辺のz = f(…)上のある曲線があって、それをx–y平面に投影した曲線のある点における微分係数と、その点におけるzの値が等しくなるような曲線を求めるという意味となる。実際にzの値と投影曲線の微分係数が等しいというのを直感的に示すのはやっかいだなと思っていた。

いろいろ探してみると、どうやらこの式形は元の関数をTaylor展開して、次数を上げて精度を高めるというものらしい。文献を見てみると、その導出過程でかなり長々とした式が出ていた。時間がある時に試してみたい。

数値解法モジュール

コード

Euler法とRunge-Kutta法で1階常微分方程式を解くモジュールを、Pythonで以下の様に書いた。

関数euler1、rk4_1の使い方はいずれも同じで、以下の通り。

  • x, y：現在の計算ステップでの独立変数xと関数値yの値
  • dx：計算の分割幅
  • func：微分方程式の右辺の関数
  • params：関数funcに与えるパラメーター

パラメーターは、呼び出し側の関数定義で設定したパラメーター名をキーとした辞書で与える。その前提で、モジュール内で**paramsとしてアンパックしている。

戻り値は実数型の数値。

使い方

モジュールの使い方の概要は以下のとおり。

  1. 右辺の関数を定義する
  2. 関数に渡すパラメーターを設定する
  3. 計算開始点と終了点のxの値、分割幅dx, yの初期値を設定する
  4. x, yの計算用リストを確保する
  5. 計算開始点のyのリストに初期値をセットする
  6. 計算実行

単純な放物線の微分方程式を例にしてモジュールの使い方を示す。まず、解の方程式として以下を考える。

(4)   \begin{equation*} y = a x^2 \quad (-1 \le x \le 1) \end{equation*}

この式を1階微分した、以下の常微分方程式を解いていく。

(5)   \begin{equation*} y'= 2ax \quad (x = -1 \Rightarrow y = 1) \end{equation*}

まず、必要なパッケージをインポート。diffeqは上で示したモジュールを収めている。

次に、微分方程式右辺の関数を定義する(以下のderivative)。この関数はパラメーターとして係数aを1つ持っている。また、パラメーターの値を辞書で定義する。

数値計算に計算開始点x_startと終了点x_end、分割幅dx、yの計算初期値を設定。

次に、予めxの計算点のリストを用意して、これに対応するyの値を計算・保存していく。計算点数は植木算で分割数＋1とし、ループ計算に用いるリストは計算終了点を除いている。結果を保存しないなら、リストではなく逐次各ステップの計算値を次のステップの入力値にする。

また、計算開始点におけるyの初期値をセットしている。

計算自体はシンプルで、定義域リストの値を1つずつ進めながら、現在の値や関数への参照、パラメーターを与えてモジュールを呼び出す。

実行例

元の解析解と併せて表示した実行結果は以下の通りで、予想通りEuler法の解は計算が進むごとに乖離が大きくなっている。Runge-Kutta法の解はかなり正確に解析解と一致している。

以下は、解析解も含めてグラフ表示するための、全体の実行コード。

 

放物線の焦点

/ tau / コメントする

概要

放物線が、その対称軸に平行な線に対して焦点を持つことを確認する。イメージ図は以下の通り。

証明

下図を考える。

放物線Cを以下の式で定義する。

(1)   \begin{equation*} C: \: y_0(x) = a x^2 \end{equation*}

任意の点(X, Y)における接線l1の式は以下の様になる。

(2)   \begin{align*} l_1: \: y_1(x) &= y_1' (X) \cdot (x - X) + Y &= 2aX (x - X) + X^2 \end{align*}

ここで、y1‘は、図中の角θを用いて以下のように表せる。

(3)   \begin{equation*} y_1'(X) = \tan \theta \quad \left( -\frac{\pi}{2} \le \theta < \frac{\pi}{2} \right) \end{equation*}

また、この点を通ってl1に垂直な直線l2の式は、この点におけるl2の微分係数をy2‘として以下のように表せる。

(4)   \begin{equation*} l_2: \: y_2(x) = y_2' (X) \cdot (x - X) + Y \end{equation*}

ここでy2‘は、図中の角θを用いて以下のように導かれる。

(5)   \begin{align*} y_2'(X) &= - \tan \left( \dfrac{\pi}{2} - 2 \theta \right) \\ &= - \frac{\cos 2 \theta}{\sin 2 \theta} \\ &= - \frac{\cos^2 \theta - \sin^2 \theta}{2 \sin \theta \cos \theta} \\ &= - \frac{-1 + \tan^2 \theta}{2 \tan \theta} \end{align*}

この式でθがπ/2よりも大きい場合でも、式形は変わらない。

ここで式(3)を考慮して、l2の式は以下のように変形される。

(6)   \begin{equation*} y_2(x) = \frac{-1 + 4 a^2 X^2}{4aX} (x - X) + a X^2 \end{equation*}

l2がy軸と交わる点のy座標値を求める。

(7)   \begin{align*} y_2(0) &= \frac{-1 + 4 a^2 X^2}{4aX} (-X) + a X^2 \\ &= \frac{X - 4 a^2 X^3 + 4a^2 X^3}{4aX} \\ &= \frac{1}{4a} \end{align*}

この値はXに寄らず一定なことから、y軸に平行な線が放物線で反射された線は、y軸上の1点を通る。

なお、入射線と反射線が直角になるときの入射線のx座標値は以下の通り。

(8)   \begin{align*} a x^2 = \frac{1}{4a} \quad \Leftrightarrow \quad x = \plusminus \frac{1}{2a} \end{align*}

斜め入射について

たとえば先の図で、入射角が右側に傾いた場合を考える。

y軸の右側（xが正の側）からの入射線の放物線に対する角度は、y軸に平行な場合に対してより浅くなることから、反射方向はより下側に傾く。したがって、反射線とy軸の交点は下側に移動する。

一方y軸の左側（xが負の側）からの入射線の放物線に対する角度は、y軸に平行な場合に対してより深くなることから、反射方向はより上側に傾く。したがって、反射線とy軸の交点は上側に移動する。

少なくともy軸の両側で反射線のy軸との交点が一致しないことから、焦点は存在しない。より正確には、入射線のy軸に対する傾きをαなどと置いて上記のような計算をするとよい。

 

ロジスティック関数

/ tau / コメントする

増加率一定

個体の増加速度が個体数に比例する場合、以下のような方程式になる。

(1)   \begin{equation*} \frac{dN}{dt} =m N \end{equation*}

この解は以下のような指数関数となり、ある程度以上から先は急激に個体が増加し、無限に増えていく。

(2)   \begin{equation*} N(t) = N_0 e^{mt} \end{equation*}

ロジスティック関数～増加に対する歯止め

増加率を一定ではなく個体数Nに応じた値となるよう、以下の関数とする。

(3)   \begin{equation*} m(N) = r \left( 1 - \frac{N}{K} \right) \end{equation*}

個体数が増えるにしたがって増加率は減少し、個体数がKのときに増加率はゼロ、それよりも多いときには個体数は減少する。

この増加率を適用した微分方程式は以下の通り。

(4)   \begin{equation*} \frac{dN}{dt} = r \left( 1 - \frac{N}{K} \right) N \end{equation*}

まず以下のように変形。

(5)   \begin{equation*} \frac{dN}{dt} = \frac{r}{K} (K - N) N \end{equation*}

変数分離形。

(6)   \begin{equation*} \frac{K}{(K - N) N} dN = r dt \end{equation*}

分数部分を以下のように変形。

(7)   \begin{equation*} \left( \frac{1}{N} - \frac{1}{K - N} \right) dN = r dt \end{equation*}

積分した一般解は以下の通り。

(8)   \begin{equation*} \ln N - \ln (K-N) = rt + const. \end{equation*}

指数関数の形に。

(9)   \begin{equation*} \frac{N}{K - N} = A e^{rt} \end{equation*}

Nをまとめる。

(10)   \begin{equation*} \frac{K - N}{N} = \frac{K}{N} - 1 = A e^{-rt} \end{equation*}

Nについて解いた一般解を得る。

(11)   \begin{equation*} N = \frac{K}{1 + A e^{-rt}} \end{equation*}

これにt=0のときの初期値N0を適用して、以下の解を得る。

(12)   \begin{equation*} N = \frac{K}{ 1 + \left( \dfrac{K}{N_0} - 1 \right) e^{-rt} } \end{equation*}

t=0のときは初期値N0、t→∞のときはN=Kに収束する。

(13)   \begin{equation*} \left\{ \begin{array}{l} N(0) = N_0 \\ \lim\limits_{t \to \infty} N(t) = K \end{array} \end{equation*}

N0=2, r=10, K=1000のグラフを描くと以下の通り。

また先の解を微分して増加速度の式を求めると以下の様になる。

パラメーターの影響

初期値N0

  1. 初期値がKより小さいときはよく見られるロジスティック曲線でKに漸近する
  2. Kと等しいときは変化しない
  3. Kより大きいときは指数関数的にKに漸近する

増加率r

増加率が大きいほど立ち上がりが急になり、小さいほど緩やかになる。

収束値K

収束値Kに向かって漸近する。

増加速度

Nをtで微分して増加速度の式を得る。

(14)   \begin{equation*} \frac{dN}{dt} = r \frac{\left( \dfrac{K}{N_0} - 1 \right) e^{-rt}}{\left( 1 + \left( \dfrac{K}{N_0} - 1 \right) e^{-rt} \right)^2} \end{equation*}

増加速度のグラフは以下の様になり、極大値を一つ持つ。

増加速度が最大となるtの値は、K/N0 – 1 = Cとおいて以下の様に得られる。

(15)   \begin{equation*} t = \frac{1}{r} \ln C = \frac{1}{r} \ln \left( \frac{K}{N_0} - 1 \right) \end{equation*}

これは以下の様にしても導ける。式(4)の値がゼロとなるのはN=K/2のときなので、これを式(12)に代入して、

(16)   \begin{equation*} 1 + \left( \dfrac{K}{N_0} - 1 \right) e^{-rt} = 2 \end{equation*}

これを解いて(15)と同じ解を得る。

 

 

折れ線と曲線の近似誤差

/ tau / コメントする

円弧の長さと弦の長さの差

半径r、内角θの円弧の長さLaは

(1)   \begin{equation*} L_a = r \theta \end{equation*}

これに対する弦の長さは

(2)   \begin{align*} L_c &= \sqrt {r^2 + r^2 - 2 r^2 \cos \theta} \\ &= r \sqrt {2 (1 - cos \theta)} \end{align*}

ここでLcとLaの比を計算すると、rを含まないθのみの関数となる。

(3)   \begin{equation*} \frac{L_c}{L_a} = \frac{\sqrt {2 ( 1 - \cos \theta ) }}{\theta} \end{equation*}

θを度単位としてこの比をグラフにすると以下のようになる。

角度変化が30度弱で誤差が2%程度、20度で0.5%程度、10度になるとかなり誤差は小さくなっている。

そこでこの付近の誤差を計算してみる。

角度 5度 10度 15度 20度 25度 30度
Lc/La 0.9997 0.9987 0.9971 0.9949 0.9921 0.9886
誤差 0.03% 0.13% 0.29% 0.51% 0.79% 1.14%

道路線形の例

道路ネットワークのリンクを折れ線で近似する場合の誤差を考える。

式((3)の誤差率の計算結果から、円弧の内角が20～30度程度以内であれば、弧と弦の長さの差はかなり小さくなりそうだと予想される。

そこで角度に対して曲線半径rを掛けて曲線部の道路の長さを計算し、さらに誤差から曲線の道路の長さと直線で近似した長さの差を計算してみる。

角度 rad 弧長(60) 弧長(710) 差分(60) 差分(710)
5 0.08727 5.236m 61.959m 1.7mm 1.9cm
10 0.17453 10.472m 123.918m 1.3cm 15.7cm
15 0.26180 15.708m 185.878m 4.5cm 53cm
20 0.34907 20.944m 247.837m 10.6cm 1.26m
25 0.43633 26.180m 309.796m 20.7cm 2.45m
30 0.52360 31.416m 371.755m 35.8cm 4.23m

ネットワークデータの用途により、許容される誤差が決まり、必要な補完点の密度が変わってくる。

たとえばメートル単位の経路案内であれば10cmオーダーの誤差は許容され、交差点での停止になると数cmなど。後者の様にミクロなコントロールや自動運転に直接データを活用する場合は、数cm程度以下の精度が必要となると想定される。

 

オイラー関数

/ tau / コメントする

定義

オイラーのトーシェント関数、オイラーのφ関数とも。totientとは数学用語で「与えられた数以下でこれと素な整数の個数」というこの関数の意味そのもの。

例

φ(6) = 2 (1, 5)

φ(12) = 4 (1, 5, 7, 11)

φ(20) = 8 (1, 3, 7, 9, 11, 13, 17, 19)

素数に対するオイラー関数値

素数はそれ自身と1以外に約数を持たないので、オイラー関数値はその値から自身の分1を減じた値となる。

φ(p) = p − 1 

2つの素数の積のオイラー関数値は、双方から1を減じた値の積になる。

φ(pq) = (p − 1)(q − 1)

これを納得するのに時間がかかった。まどろっこしいが以下のように考えた。

まず以下のことに気づく。

(1)   \begin{equation*} (p - 1)(q - 1) = pq - p - q + 1 = pq - (p - 1) - (q - 1) - 1 \end{equation*}

次に2つの素数5と3で以下のように考えてみる。

まずφ(15)に関して。15以下の数のうち15と互いに素であるのは以下の8個。

1, 2, 4, 7, 8, 11, 13, 14

逆に互いに素でない(15と公約数を持つ)数は以下の6個。

3, 5, 6, 9, 10, 12

ここで以下のように考える。

  • φ(5)に関してはカウントされる数は1, 2, 3, 4の4つで、φ(3)に関してはカウントされる数は1, 2の2つ
  • これらの数に任意の自然数を掛けて15以下の自然数を作り出すことができる
  • ところが3×5は素数の掛け算であり、上記の数に3や5以外の数をかけてもその結果は3×5と互いに素であることに変わりはない
  • φ(3)でカウントする1, 2に対して5を掛ければ、3×5以下でこれと互いに素でない数が得られる
  • またφ(5)でカウントする1, 2, 3, 4に対して3を掛ければ、3×5以下でこれと互いに素でない数が得られる
  • しかも3と5はいずれも素数で互いに素なので、上記の結果に重複は生じない

φ(3)に関しては、対象となる数に5を掛けて以下が得られる。

1, 2　→　5, 10

φ(5)に関しては、対象となる数に3を掛けて以下が得られる。

1, 2, 3, 4　→　3, 6, 9, 12

以上を併せて、15に対して1を除いて互いに素でない6個の数が得られる。

すなわち、φ(pq)は以下のようにして得られる。

  • 重複しないp − 1個、q − 1個の数がpqと互いに素でない
  • pq自身1個もpqと互いに素でない

したがって以下のようになる。

(2)   \begin{equation*} \phi (pq) = pq - (p - 1) - (q - 1) - 1 = (p - 1)(q - 1) \end{equation*}

 

主成分分析の定式化

/ tau / コメントする

概要

主成分分析では、複数の特徴量を持つデータセットから、そのデータセットの特徴を最もよく表す特徴量軸を発見していく。

ここで「特徴を最もよく表す」ことを数学的に「最も分散が大きくなる」と定義する。そして、分散が最も大きくなるような方向を探すことを目的とする。

ある軸に沿った分散が大きくなるということは、その軸に沿った性質のバリエーションが多いことになる。逆に分散が小さい場合は、その性質を表す数量によっては各データの特徴の違いが判別しにくい。

主成分分析では、分散が最大となるような軸の方向を発見することが目的となる。この軸は元の特徴量の線形和で表現されるもので、各特徴量の係数は、それぞれの特徴量の寄与を表す。

(1)   \begin{align*} \boldsymbol{v} &= a_1 \boldsymbol{x}_1 + \ldosts + a_m \boldsymbol{x_m} \\ &= a_1 \left( \begin{array}{c} \x_1 \\ 0 \\ \vdots \\ 0 \end{array} \right) + \cdots + a_m \left( \begin{array}{c} 0 \\ \vdots \\ 0 \\ x_m \end{array} \right) \\ v &= | \boldsymbol{v} | = a_1 x_1 + \cdots + a_m x_m \end{align*}

以後、複数の特徴量を持つデータを、特徴量を成分とするベクトルでx表し、多数のベクトルデータxiがデータセットを構成しているとする。

最大化すべき分散の導出

多数のデータの中のデータiが空間内の点に対応し、その位置ベクトルをxiであるとする。このxiの成分が特徴量に対応する。長さが1のあるベクトルdが与えられたとき、xiのdへの射影の長さは以下のように計算される。

(2)   \begin{align*} x_{i | \boldsymbol{d}} = {\boldsymbol{x}_i}^T \boldsymbol{d} = \boldsymbol{d}^T \boldsymbol{x}_i \quad (| \boldsymbol{d} | = 1) \end{align*}

たとえば特徴量が2つなら、2次元で以下のような計算になる。

(3)   \begin{align*} \boldsymbol{x}_i = \left( \begin{array}{C} x_{i1} \\ x_{i2} \end{array} \right) , \quad \boldsymbol{d} = \left( \begin{array}{C} d_1 \\ d_2 \end{array} \right) \end{align*}

(4)   \begin{align*} x_{i | \boldsymbol{d}} = ( d_1 , d_2 ) \left( \begin{array}{c} x_{i1} \\ x_{i2} \end{array} \right) = ( d_1 x_{i1} + d_2 x_{i2} ) \end{align*}

n個のデータ(i = 1～n)について、射影の平均は以下のように計算される。これは全データのベクトルdの方向に沿った値の平均となる。

(5)   \begin{align*} E( x_{i | \boldsymbol{d}} ) = E \left( \boldsymbol{d}^T \boldsymbol{x}_i \right) = \boldsymbol{d}^T E \left( \boldsymbol{x}_i \right) = \boldsymbol{d}^T \boldsymbol{\mu}_i \end{align*}

これも2次元の場合で確認すると以下の通り。

(6)   \begin{align*} E(x_{i | \boldsymbol{d}} ) &= E\left[ (d1, d2) \left( \begin{array}{c} x_{i1} \\ x_{i2} \end{array} \right) \right] = (d_1, d_2) \left( \begin{array}{c} E(x_{i1}) \\ E(x_{i2}) \end{array} \right) \\ &= (d_1, d_2) \left( \begin{array}{c} \mu_{i1} \\ \mu_{i2} \end{aray} \right) \end{align*}

式(5)を使ってベクトルdの方向に沿ったデータの分散を計算する。

(7)   \begin{align*} V( x_{i | \boldsymbol{d}} ) &= V \left( \boldsymbol{d}^T \boldsymbol{x}_i \right) \\ &= E \left[ \left( \boldsymbol{d}^T \boldsymbol{x}_i - E \left( \boldsymbol{d}^T \boldsymbol{x}_i \right) \right)^2 \right] \\ &= E \left[ \left( {\boldsymbol{d}}^T \left( \boldsymbol{x}_i - E(\boldsymbol{x}_i) \right) \right)^2 \right] \\ &= E \left[ {\boldsymbol{d}}^T (\boldsymbol{x}_i - \boldsymbol{\mu}_i ) (\boldsymbol{x}_i - \boldsymbol{\mu}_i )^T \boldsymbol{d} \right] \\ &= \boldsymbol{d}^T E\left[ (\boldsymbol{x}_i - \boldsymbol{\mu}_i ) (\boldsymbol{x}_i - \boldsymbol{\mu}_i )^T \right] \boldsymbol{d} \\ &= \boldsymbol{d}^T \boldsymbol{\Sigma} \boldsymbol{d} \end{align*}

中央の平均の項が共分散行列Σとなっていることに留意。これより、あるベクトルが与えられたとき、その方向に沿った全データの成分の分散が、そのベクトルと元のデータの共分散行列を使って求めることができる。

こちらを2次元で確認すると以下の通り。

(8)   \begin{align*} &E\left[ (\boldsymbol{x}_i - \boldsymbol{\mu}_i ) (\boldsymbol{x}_i - \boldsymbol{\mu}_i )^T \right] \\ &= E \left[ \left( \begin{array}{c} x_{i1} - \mu_1 \\ x_{i2} - \mu_2 \end{array} \right) (x_{i1} - \mu_1, x_{i2} - \mu_2) \right] \\ &= \left[ \begin{array}{cc} (x_{i1} - \mu_1)^2 & (x_{i1} - \mu_1)(x_{i2} - \mu_2) \\ (x_{i2} - \mu_2)(x_{i1} - \mu_1) & (x_{i2} - \mu_2)^2 \end{array} \right] \end{align*}

分散の最大化

式(8)で計算された分散が最大となるようにベクトルdの方向を決定する。このとき、dの大きさが1であるという制約条件があるため、問題は制約条件付きの最大化問題となる。

(9)   \begin{gather*} {\rm max} \quad \boldsymbol{d}^T \boldsymbol{\Sigma} \boldsymbol{d} \quad \rm{s.t.} \; | \boldsymbol{d} | = 1 \end{gather*}

これをLagrangeの未定乗数法で解いていく。。

(10)   \begin{gather*} L( \boldsymbol{d}, \lambda ) = \boldsymbol{d}^T \boldsymbol{\Sigma} \boldsymbol{d} - \lambda (|\boldsymbol{d}|^2 - 1) = 0 \\ \frac{\partial L}{\partial d_i} = 0 \quad ( {\rm for \; all} \; i ) \end{gather*}

Lagrange関数の第1項については、

(11)   \begin{align*} \boldsymbol{d}^T \boldsymbol{\Sigma} &= \begin{array}{ccc} ( & d_1 V_1 + \cdots + d_n C_{n1} & , \\ & \vdots & ,\\ & d_1 C_{1j} + \cdots + d_n C_{n, j} & , \\ & \vdots & ,\\ & d_1 C_{1n} + \cdots + d_n V_n & ) \end{array} \end{align*}

より、以下のような長い式になる。

(12)   \begin{align*} \boldsymbol{d}^T \boldsymbol{\Sigma d} &= \begin{array}{c} {d_1}^2 V_1 + \cdots + d_j d_1 C_{j1} + \cdots + d_n d_1 C_{n1} + \\ \vdots \\ d_1 d_j C_{1j} + \cdots + {d_j}^2 V_j + \cdots + d_n d_j C_{nj} + \\ \vdots \\ d_1 d_n C_{1n} + \cdots + d_j d_n C_{jn} + \cdots + {d_n}^2 V_n \end{array} \end{align*}

また第2項の括弧の中については以下のようになる。

(13)   \begin{align*} | \boldsymbol{d} |^2 - 1 = ( {d_1}^2 + \cdots + {d_j}^2 + \cdots + {d_n}^2 ) -1 \end{align*}

これらを前提に、Lをdjで微分すると以下のようになる。

(14)   \begin{align*} 2 d_1 C_{1j} + \cdots + 2 d_j V_{j} + \cdots 2 d_n C_{2j} - 2 \lambda d_j = 0 \end{align*}

全てのdjについて考慮した連立方程式を行列形式で表すと以下のようになる。

(15)   \begin{gather*} \boldsymbol{\Sigma d} = \lambda \boldsymbol{d} \\ | \boldsymbol{d} | = 1 \end{gather*}

1つ目の式は共分散行列に関する固有値問題の式で、di (i=1～n)とλのn+1個の変数に対してn個の式となる。これに先ほど脇に置いていたdの大きさに関する制約式を加えて式の数もn+1個となり、dとλが求められる。

特徴量が2つの場合

特徴量が2つの場合を考え、以下のように記号を定義する。

(16)   \begin{align*} \boldsymbol{\Sigma} = \left( \begin{array}{cc} \sigma_{11} & \sigma_{12} \\ \sigma_{21} & \sigma_{22} \end{array} \right) , \quad \boldsymbol{d} = \left( \begin{array}{c} d_1 & d_2 \end{array} \right) \end{align*}

このとき、分散を最大化する方向の単位ベクトルdを求める方程式は以下のようになる。

(17)   \begin{equation*} \left\{ \begin{array}{l} \left( \begin{array}{cc} \sigma_{11} & \sigma_{12} \\ \sigma_{21} & \sigma_{22} \end{array} \right) \left( \begin{array}{c} d_1 & d_2 \end{array} \right) = \lambda \left( \begin{array}{c} d_1 & d_2 \end{array} \right) \\ {d_1}^2 + {d_2}^2 = 1 \end{array} \right. \end{equation*}

1つ目の式を解くと、

(18)   \begin{equation*} \left\{ \begin{array}{l} \sigma_{11} d_1 + \sigma_{12} d_2 = \lambda d_1 \\ \sigma_{21} d_1 + \sigma_{22} d_2 = \lambda d_2 \end{array} \rhight. \end{equation*}

この方程式は不定なのでd1、d2それぞれは求められないが、μ = d2/d1は計算できる。これは固有ベクトルの方向が定まる。具体的には下記の通り。

(19)   \begin{gather*} \left\{ \begin{array}{l} \sigma_{11} + \sigma_{12} \mu = \lambda \\ \sigma_{21} + \sigma_{22} \mu = \lambda \mu \end{array} \right. \\ \lambda = \sigma_{11} + \sigma_{12} \mu = \frac{\sigma_{21}}{\mu} + \sigma_{22} \\ \sigma_{12} \mu^2 + ( \sigma_{11} - \sigma_{22} ) \mu - \sigma_{21} = 0 \end{gather*}

これを解いてベクトルdの方向が定まる。これに制約条件|d|2 ＝１を加味することで、大きさ1の単位ベクトルとしてdが決定される。

この解き方は最大化問題ではないので、連立方程式から2つの固有ベクトルと固有値が求まる。

第2主成分以降

一般的な固有値問題では、元の変数と同じ数の固有ベクトルと固有値のセットが求まるが、最大化問題として解いた場合には主成分が1つだけ求まる。

scikit-learnのPCAインスタンス生成時にn_componentsで主成分の数に制約をかけることができるが、このことから、PCA.fit()の実行時には連立方程式を解いているのではなく、最大化問題で1つずつ主成分を計算しているのではないかと思われる。

第2主成分以降の計算についての紹介はあまり見られないが、以下の手順と考えらえれる。

  1. 各データについて、第1主成分の方向への射影を計算
  2. その射影の符号を逆にしたベクトルを各データに加える
  3. これで第1主成分に沿ったばらつきが全てゼロになるので、残りの成分の中で最大となるベクトルの方向を計算し、第2主成分とする
  4. 以上を繰り返し、順次最大主成分沿いの情報を消しながら、各主成分を計算

主成分の意味

主成分の意味の一つに、元のデータを構成する成分という捉え方がある。

たとえば特徴量の数がnである元データXがあり、主成分の数をm(<= n)でモデルを構築するとする。scikit-learnでPCAのインスタンスを生成するのにn_components=mと指定し、fit(X)を実行すると、m個の主成分が生成される。この主成分は共分散行列に対する固有ベクトルであり、要素数n個（特徴量数に等しい）の1次元配列がm行（主成分の数に等しい）並んだ2次元配列として、PCAインスタンスのプロパティーcomponents_に保存される。

(20)   \begin{equation*} \tt{components\_} = \left[ \begin{array}{ccc} (p_{0, 0} & \cdots & p_{0, n-1} ) \\ & \vdots &\\ (p_{m-1, 0} & \cdots & p_{m-1, n-1}) \end{array} \right] = \left[ \begin{array}{c} \boldsymbol{p}_0 \\ \vdots \\ \boldsymbol{p}_m \end{array} \right] \end{equation*}

元のデータは、各主成分（固有ベクトル）の重み付き和として表現される。

(21)   \begin{equation*} \boldsymbol{x} = (x_0, ..., x_n) = a_0 \boldsymbol{p}_0 + a_1 \boldsymbol{p}_1 + a_2 \boldsymbol{p}_2 + \cdots \end{equation*}

この様子を2次元で示したのが以下の図で、直行する2つの主成分から元データの１つxが定まる。

xの主成分1、2の方向の大きさはxの各主成分に対する射影で、それらの長さはxと各主成分の内積で得られる。

(22)   \begin{equation*} \boldsymbol{x} = (x_0, ..., x_n) = ( \boldsymbol{x} \cdot \boldsymbol{p}_0 ) \boldsymbol{p}_0 + ( \boldsymbol{x} \cdot \boldsymbol{p}_1 ) \boldsymbol{p}_1 + ( \boldsymbol{x} \cdot \boldsymbol{p}_2 ) \boldsymbol{p}_2 + \cdots \end{equation*}