3×3行列の固有値・固有ベクトルを計算する。
固有値・固有ベクトルがすべて実数となることが前もって分かっている行列用。
Compute the eigenvalues and eigenvectors of a 3×3 matrix,
for which all the eigenvalues and eigenvectors
are known to be real numbers.
行列
\[\begin{equation}
\myvector{A}=
\begin{pmatrix}
a & b & c \\
d & e & f \\
g & h & i
\end{pmatrix}
\label{eq.matrixDef}
\end{equation}\]
に対する固有値問題は3×3の単位行列\(\myvector{I}\)を用いて
\[\begin{eqnarray}
det(\lambda \myvector{I}-\myvector{A})
&=& \lambda^3
\nonumber \\
& & -(a+e+i)\lambda^2
\nonumber \\
& & +(ae+ei+ia-bd-cg-fh)\lambda
\nonumber \\
& & -(aei+bfg+cdh-afh-ecg-ibd)
\nonumber \\
&=& 0
\label{eq.eigenEquation}
\end{eqnarray}\]
で与えられる。
これを3次方程式の解の公式
(関数eqsol3のマニュアル参照)
を使って解く。
The eigenvalue problem for a matrix \(\myvector{A}\)
(Eq. \ref{eq.matrixDef}) is given by Eq. (\ref{eq.eigenEquation}),
where \(\myvector{I}\) is a 3×3 unit matrix.
This equation is solved using the solution formula for cubic euqations
(see the
documentation of function eqsol3).
2. 規格化
(Normalization)
関数eqsol3では3次の項の係数が他の項の係数に比べて桁で小さい場合、
実質的に2次方程式になっていると判断して
プログラムをエラー終了する仕様になっている。
(\ref{eq.eigenEquation})式をそのまま関数eqsol3に渡した場合、
\(a\)-\(i\)の値が大きいと2次方程式であると判定されてエラーになりやすい。
これを避けるため、実際には\(a\)-\(i\)の中の絶対値最大値
\[\begin{equation}
A_{amp}\equiv \max\{|a|,|b|,|c|,|d|,|e|,|f|,|g|,|h|,|i|\}
\label{eq.Aamp}
\end{equation}\]
を用いて
\[\begin{equation}
\myvector{A}=A_{amp}\myvector{A’}, \hspace{1em}
\myvector{A’}=
\begin{pmatrix}
a’ & b’ & c’ \\
d’ & e’ & f’ \\
g’ & h’ & i’
\end{pmatrix}
\label{eq.Adash}
\end{equation}\]
と変形した上で\(\myvector{A’}\)の固有値を求め、
最後にそれを\(A_{amp}\)倍するアプローチを採っている。
When the coefficient for the cubic term is by orders of magnitudes less than
those for the other terms,
function eqsol3 regards the equation as a quadratic equation
and finishes the program as an error.
This error readily occurs when Eq. (\ref{eq.eigenEquation}) is
directly passed to function eqsol3
and the values of \(a\)-\(i\) are large.
To avoid this situation, function eigen3_r modifies the matrix as
Eq. (\ref{eq.Adash}), where \(A_{amp}\) is the absolute maximum value
of \(a\)-\(i\) (Eq. \ref{eq.Aamp}),
solves for the eigenvalues of \(\myvector{A}’\),
and multiplies the results with \(A_{amp}\).
3. 縮退の無い固有値に対する固有ベクトル
(Eigenvectors corresponding to undegenerated eigenvalues)
固有値\(\lambda\)に対する固有ベクトル\(\myvector{V}^T=(V_x,V_y,V_z)\)を
\[\begin{equation}
V_x=\sin\theta\cos\phi
\label{eq.Vx.pole}
\end{equation}\]
\[\begin{equation}
V_y=\sin\theta\sin\phi
\label{eq.Vy.pole}
\end{equation}\]
\[\begin{equation}
V_z=\cos\theta
\label{eq.Vz.pole}
\end{equation}\]
\[\begin{equation}
0\leq \theta < \pi, \hspace{1em} 0 \leq \phi < \pi
\label{eq.theta_phi_range}
\end{equation}\]
と極座標表示する。
固有ベクトルは定数倍しても固有ベクトルであるので
一般に固有ベクトルは一意に決まらない。
そこで上の表記では長さを1に規格化し、
\(-1\)倍の不確定を除くために\(\theta,\phi\))の範囲にも制約を設けている。
Let \(\myvector{V}^T=(V_x,V_y,V_z)\) be
the eigenvector corresponding to an eigenvalue \(\lambda\).
Let us introduce a polar-coordinate expression of this vector
(Eqs. \ref{eq.Vx.pole}-\ref{eq.Vz.pole}),
where \(0\leq \theta < \pi\) and \(0 \leq \phi < \pi\).
Here, the length of the eigenvector is assumed to be unity
and the ranges of \(\theta\) and \(\phi\) are limited
to avoid the ambiguity of eigenvectors
such that an eigenvector multiplied by a constant value is an eigenvector.
固有ベクトルを求める方程式は行列形式では
\[\begin{equation}
\begin{pmatrix}
a & b & c \\
d & e & f \\
g & h & i
\end{pmatrix}
\begin{pmatrix}
V_x \\
V_y \\
V_z
\end{pmatrix}
=\lambda
\begin{pmatrix}
V_x \\
V_y \\
V_z
\end{pmatrix}
\label{eq.eigenvector}
\end{equation}\]
と書けるが、これを成分ごとに書き下すと
\[\begin{equation}
a V_x + b V_y + c V_z = \lambda V_x
\label{eq.eigenvector.x}
\end{equation}\]
\[\begin{equation}
d V_x + e V_y + f V_z = \lambda V_y
\label{eq.eigenvector.y}
\end{equation}\]
\[\begin{equation}
g V_x + h V_y + i V_z = \lambda V_z
\label{eq.eigenvector.z}
\end{equation}\]
となる。これらに(\ref{eq.Vx.pole})-(\ref{eq.Vz.pole})を代入すると
\[\begin{equation}
a \sin\theta\cos\phi + b \sin\theta\sin\phi + c \cos\theta
= \lambda \sin\theta\cos\phi
\label{eq.eigenvector.x.polar}
\end{equation}\]
\[\begin{equation}
d \sin\theta\cos\phi + e \sin\theta\sin\phi + f \cos\theta
= \lambda \sin\theta\sin\phi
\label{eq.eigenvector.y.polar}
\end{equation}\]
\[\begin{equation}
g \sin\theta\cos\phi + h \sin\theta\sin\phi + i \cos\theta
= \lambda \cos\theta
\label{eq.eigenvector.z.polar}
\end{equation}\]
となり、両辺を\(\cos\theta\)で割って
\[\begin{equation}
a \tan\theta\cos\phi + b \tan\theta\sin\phi + c
= \lambda \tan\theta\cos\phi
\label{eq.eigenvector.x.polar.modified}
\end{equation}\]
\[\begin{equation}
d \tan\theta\cos\phi + e \tan\theta\sin\phi + f
= \lambda \tan\theta\sin\phi
\label{eq.eigenvector.y.polar.modified}
\end{equation}\]
\[\begin{equation}
g \tan\theta\cos\phi + h \tan\theta\sin\phi + i
= \lambda
\label{eq.eigenvector.z.polar.modified}
\end{equation}\]
を得る。
(\label{eq.eigenvector.x.polar.modified})
-
(\label{eq.eigenvector.z.polar.modified})
式をそれぞれ\(\tan\theta\)について解けば
\[\begin{equation}
\tan\theta
=\frac{c}{\lambda\cos\phi-a\cos\phi-b\sin\phi}
=\frac{f}{\lambda\sin\phi-d\cos\phi-e\sin\phi}
=\frac{\lambda-i}{g\cos\phi+h\sin\phi}
\label{eq.tantheta}
\end{equation}\]
が得られる。
(\ref{eq.tantheta})式の右辺には\(\theta\)が含まれないので、
(\ref{eq.tantheta})式は\(\phi\)が得られたときに\(\theta\)を陽に求める式
になっている。
The eigenvector satisfies Eq. (\ref{eq.eigenvector}),
which can be rewritten in scalar format as
Eqs. (\ref{eq.eigenvector.x})-(\ref{eq.eigenvector.z}).
Inserting Eqs. (\ref{eq.Vx.pole})-(\ref{eq.Vz.pole}) into these equations
results in (\ref{eq.eigenvector.x.polar})-(\ref{eq.eigenvector.z.polar}),
and dividing them by \(\cos\theta\) results in Eqs.
(\label{eq.eigenvector.x.polar.modified})
-
(\label{eq.eigenvector.z.polar.modified}).
Solving each of them for \(\tan\theta\) gives Eq. (\ref{eq.tantheta}).
Note that \(\theta\) is absent in this equation,
meaning that \(\theta\) can directly be calculated
from Eq. \ref{eq.tantheta} once \(\phi\) is obtained.
(\ref{eq.tantheta})式の右辺の1つ目と3つ目が等しくなる条件から
\[\begin{equation}
c(g\cos\phi+h\sin\phi)=(\lambda-i)(\lambda\cos\phi-a\cos\phi-b\sin\phi)
\label{eq.phi.derive1a}
\end{equation}\]
が得られ、(\ref{eq.tantheta})式の右辺の2つ目と3つ目が等しくなる条件から
\[\begin{equation}
f(g\cos\phi+h\sin\phi)=(\lambda-i)(\lambda\sin\phi-d\cos\phi-e\sin\phi)
\label{eq.phi.derive1b}
\end{equation}\]
が得られる。これらの式の両辺を\(\cos\phi\)で割ると
\[\begin{equation}
cg+ch\tan\phi=(\lambda-i)(\lambda-a)-(\lambda-i)b\tan\phi
\label{eq.phi.derive2a}
\end{equation}\]
\[\begin{equation}
fg+fh\tan\phi=-(\lambda-i)d+(\lambda-i)(\lambda-e)\tan\phi
\label{eq.phi.derive2b}
\end{equation}\]
となり、これらをそれぞれ\(\tan\phi\)について解けば
\[\begin{equation}
\tan\phi
=\frac{(\lambda-i)(\lambda-a)-cg}{(\lambda-i)b+ch}
=\frac{(\lambda-i)d+fg}{(\lambda-i)(\lambda-e)-fh}
\label{eq.tanphi}
\end{equation}\]
が得られる。
(\ref{eq.tanphi})式の右辺には\(\phi\)が含まれていない。
したがって固有値が求まれば(\ref{eq.tanphi})式によって\(\phi\)が求まり、
更にその結果を用いて(\ref{eq.tantheta})式により\(\theta\)が求まり、
それらを用いて(\ref{eq.Vx.pole})-(\ref{eq.Vz.pole})式により
固有ベクトルが求まるという流れである。
Requiring the 1st and 3rd right hand sides of Eq. (\ref{eq.tantheta})
being equal, we obtain Eq. (\ref{eq.phi.derive1a}).
Requiring the 2nd and 3rd right hand sides of Eq. (\ref{eq.tantheta})
being equal, we obtain Eq. (\ref{eq.phi.derive1b}).
Dividing the both sides of
Eqs. (\ref{eq.phi.derive1a}) and (\ref{eq.phi.derive1b}) by \(\cos\phi\)
results in Eqs. (\ref{eq.phi.derive2a}) and (\ref{eq.phi.derive2b}),
which can be solved for \(\tan\phi\) as Eq. (\ref{eq.tanphi}).
Note that \(\phi\) is absent in the right hand side of this equation.
Therefore, once an eigenvalue is obtained,
\(\phi\) is obtained from Eq. (\ref{eq.tanphi})
and then \(\theta\) is obtained from Eq. \(\theta\).
Finally, the eigenvector is obtained by inserting these results into
Eqs. (\ref{eq.Vx.pole})-(\ref{eq.Vz.pole}).
のように場合分けを行って求める。
(\ref{eq.tantheta})式を用いて\(\tan\theta\)を求める際も
同様の考え方で計算式が選択される。
Eq. (\ref{eq.tanphi}) consists of two expressions;
they would give the same result,
as requiring the two expressions being equal
results in an eigenvalue equation for \(\lambda\).
Function eigen3_r evaluates
the absolute values of the denominators and numerators of
the two expressions,
and adopts the expression that gives the large absolute values.
This procedure is formally explained as follows.
Let \(N_1\) and \(D_1\) be
the numerator and denominator of the 1st expression
given by Eqs. (\ref{eq.N1}) and (\ref{eq.D1}), respectively.
In the same way, let \(N_2\) and \(D_2\) be
the numerator and denominator of the 2nd expression
given by Eqs. (\ref{eq.N2}) and (\ref{eq.D2}), respectively.
If \(\max\{|N_1|,|D_1|\} \geq \max{|N_2|,|D_2|}\),
the program computes as \(\tan\phi=N_1/D_1\);
otherwise the program computes as \(\tan\phi=N_2/D_2\).
Selection of the expression for \(\tan\theta\)
from the three candidates of Eq. (\ref{eq.tantheta})
is based on the same approach.
Finally, \(\theta\) and \(\phi\) are adjusted to satisfy
Eq. (\ref{eq.theta_phi_range})
using the following two procedures:
add or subtract \(\pi\) to/from \(\theta\), or
add or subtract \(\pi\) to/from \(\phi\)
and at the same time replace \(\theta\) with \(\pi-\theta\).
These procedures cause the eigenvector to be multiplied with \(-1\)
(see Eqs. \ref{eq.Vx.pole}-\ref{eq.Vz.pole}),
and thus the result is still an eigenvector.
Actual steps in the function are as follows.
When \(\phi<0\),
\(\pi\) is added to \(\phi\)
and \(\theta\) is replaced by \(\pi-\theta\).
When \(\phi\geq\pi\),
\(\pi\) is subtracted from \(\phi\)
and \(\theta\) is replaced by \(\pi-\theta\).
When \(\theta<0\),
\(\pi\) is added to \(\theta\).
When \(\theta\geq\pi\),
\(\pi\) is subtracted from \(\theta\).
4. 縮退した固有値の扱い
(Treatment of degenerated eigenvalues)
固有値が縮退している場合には固有ベクトルは一意に求まらない。
関数eigen3_rでは現在のところ、
縮退した固有値に対する固有ベクトルを求めることはせず、
固有ベクトルの値として(0,0,0)を返す仕様にしている。
なお縮退しているか否かの判定には
固有値を\(A_{amp}\)で割った値に対して関数doublecmpを用いている。
When eigenvalues are degenerated,
eigenvectors are not obtained uniquely.
Currently, function eigen3_r does not solve for eigenvectors
corresponding to degenerated eigenvalues,
and instead returns (0,0,0) as the eigenvalues.
The degenerated eigenvalue are identified based on
eigenvalues normalized by \(A_{amp}\) with function doublecmp.