diagonalize2d.py ダウンロード/コピー
diagonalize2d.py
diagonalize2d.py
1#!/usr/bin/env python3
2"""
3diagonalize2d.py
4
5概要:
6 2x2実数行列の対角化を可視化するスクリプト。
7
8詳細説明:
9 このスクリプトは以下の機能を提供します:
10 1. 元の基底ベクトル (a1, a2) とそれらの行列Sによる像 (S a1, S a2) を描画します。
11 2. 2x2実数行列Sの固有値と固有ベクトルを計算します。
12 3. 行列が実数上で対角化可能でない場合、メッセージを表示して終了します。
13 4. 行列が実数上で対角化可能である場合、実固有ベクトル (v1, v2) とそれらの像 (S v1, S v2) を描画します。
14 5. 対角化前後の基底ベクトル間の角度を出力します。
15
16解釈:
17 - Sが対称行列の場合、対角化基底は直交するので、角度は90度を保ちます。
18 - Sが非対称であるが実数上で対角化可能な場合、対角化基底は一般に直交しません。
19 - 複素固有値が現れる場合、実数の主軸基底は存在しません。
20
21Usage:
22 python diagonalize2d.py S11 S12 S21 S22
23
24関連リンク:
25 :doc:`diagonalize2d_usage`
26"""
27
28from __future__ import annotations
29
30import argparse
31import math
32
33import matplotlib.pyplot as plt
34import numpy as np
35
36
37EPS = 1e-10
38
39
40def unit(v: np.ndarray) -> np.ndarray:
41 """
42 概要:
43 与えられたベクトルを単位ベクトルに正規化する。
44
45 詳細説明:
46 ベクトルのL2ノルムが非常に小さい（EPS以下）場合、ゼロベクトルとみなしValueErrorを発生させる。
47
48 :param v: 正規化する2Dベクトル（numpy.ndarray形式）。
49 :returns: 正規化された単位ベクトル（numpy.ndarray形式）。
50 :raises ValueError: ゼロベクトルを正規化しようとした場合。
51 """
52 n = np.linalg.norm(v)
53 if n < EPS:
54 raise ValueError("Zero vector cannot be normalized.")
55 return v / n
56
57
58def angle_deg(v: np.ndarray) -> float:
59 """
60 概要:
61 2Dベクトルが+X軸となす角度を度数で計算する。
62
63 詳細説明:
64 角度の範囲は (-180, 180] 度となる。
65
66 :param v: 角度を計算する2Dベクトル（numpy.ndarray形式）。
67 :returns: +X軸となす角度（度数）。
68 """
69 return math.degrees(math.atan2(float(v[1]), float(v[0])))
70
71
72def angle_between_deg(v1: np.ndarray, v2: np.ndarray) -> float:
73 """
74 概要:
75 2つの2Dベクトル間の最小角度を度数で計算する。
76
77 詳細説明:
78 計算される角度の範囲は [0, 180] 度となる。内積とアークコサインを用いて計算される。
79
80 :param v1: 1つ目の2Dベクトル（numpy.ndarray形式）。
81 :param v2: 2つ目の2Dベクトル（numpy.ndarray形式）。
82 :returns: 2つのベクトル間の最小角度（度数）。
83 """
84 u1 = unit(v1)
85 u2 = unit(v2)
86 c = float(np.clip(np.dot(u1, u2), -1.0, 1.0))
87 return math.degrees(math.acos(c))
88
89
90def is_parallel(v1: np.ndarray, v2: np.ndarray, tol: float = 1e-9) -> bool:
91 """
92 概要:
93 2つの2Dベクトルが平行であるかをチェックする。
94
95 詳細説明:
96 2Dベクトル `(x1, y1)` と `(x2, y2)` が平行である場合、`x1*y2 - y1*x2` (2D外積のZ成分) は0になる。
97 この値が指定された許容誤差 `tol` よりも小さい絶対値であれば、平行とみなす。
98
99 :param v1: 1つ目の2Dベクトル（numpy.ndarray形式）。
100 :param v2: 2つ目の2Dベクトル（numpy.ndarray形式）。
101 :param tol: 平行とみなすための許容誤差。
102 :returns: 2つのベクトルが平行であればTrue、そうでなければFalse。
103 """
104 return abs(float(v1[0] * v2[1] - v1[1] * v2[0])) < tol
105
106
107def classify_real_diagonalizable(S: np.ndarray):
108 """
109 概要:
110 2x2行列が実数上で対角化可能であるかを判定する。
111
112 詳細説明:
113 行列Sの固有値と固有ベクトルを計算し、以下の条件で対角化可能性を判定する。
114 1. 固有値が複素数を含む場合、実数上で対角化不可能。
115 2. 実数固有値のみの場合でも、固有ベクトルが線形独立でない場合、実数上で対角化不可能。
116 それ以外の場合は、実数上で対角化可能と判定する。
117
118 :param S: 判定対象の2x2実数行列（numpy.ndarray形式）。
119 :returns: (tuple)
120 - `diagonalizable_over_R` (bool): 実数上で対角化可能であればTrue、そうでなければFalse。
121 - `message` (str): 対角化可能性に関する説明メッセージ。
122 - `eigenvalues` (numpy.ndarray): 計算された固有値。複素数を含む場合がある。
123 - `eigenvectors` (numpy.ndarray): 計算された固有ベクトル。複素数を含む場合がある。
124 """
125 vals, vecs = np.linalg.eig(S)
126
127 if np.max(np.abs(vals.imag)) > 1e-9:
128 return False, "Complex eigenvalues: no real principal axes.", vals, vecs
129
130 vals = vals.real
131 vecs = vecs.real
132
133 rank = np.linalg.matrix_rank(vecs)
134 if rank < 2:
135 return False, "Real eigenvalues found, but eigenvectors are not independent: not diagonalizable over R.", vals, vecs
136
137 return True, "Diagonalizable over the real numbers.", vals, vecs
138
139
140def sort_eigensystem(vals: np.ndarray, vecs: np.ndarray):
141 """
142 概要:
143 固有値と固有ベクトルをソートし、表示用に正規化する。
144
145 詳細説明:
146 固有値を昇順にソートし、それに対応するように固有ベクトルの順序も変更する。
147 さらに、各固有ベクトルの最初の非ゼロ成分が正になるように符号を調整し、表示を統一する。
148
149 :param vals: ソート対象の固有値（numpy.ndarray形式）。
150 :param vecs: ソート対象の固有ベクトル（numpy.ndarray形式、各列が1つの固有ベクトル）。
151 :returns: (tuple)
152 - `vals` (numpy.ndarray): ソートされ、符号調整後の固有値。
153 - `vecs` (numpy.ndarray): ソートされ、符号調整後の固有ベクトル。
154 """
155 idx = np.argsort(vals)
156 vals = vals[idx]
157 vecs = vecs[:, idx]
158
159 for j in range(2):
160 v = unit(vecs[:, j])
161 if abs(float(v[0])) > EPS:
162 if float(v[0]) < 0:
163 v = -v
164 else:
165 if float(v[1]) < 0:
166 v = -v
167 vecs[:, j] = v
168
169 return vals, vecs
170
171
172def plot_vectors(
173 ax,
174 basis_vectors,
175 image_vectors,
176 basis_labels,
177 image_labels,
178 title,
179):
180 """
181 概要:
182 指定された軸上に基底ベクトルとそれらの像を描画する。
183
184 詳細説明:
185 プロットの軸範囲は自動的に調整され、全てのベクトルが見えるように設定される。
186 基底ベクトルは黒い実線矢印、像ベクトルは赤い実線矢印で表示される。
187 軸の目盛り、ラベル、タイトルは非表示になり、グリッドと枠線も描画されない。
188
189 :param ax: 描画対象のMatplotlib Axesオブジェクト。
190 :param basis_vectors: 基底ベクトルのリスト（各要素はnumpy.ndarray形式）。
191 :param image_vectors: 基底ベクトルの像のリスト（各要素はnumpy.ndarray形式）。
192 :param basis_labels: 基底ベクトルのラベルのリスト（文字列）。
193 :param image_labels: 像ベクトルのラベルのリスト（文字列）。
194 :param title: プロットのタイトル（未使用、コード内でコメントアウトされている）。
195 :returns: なし
196 """
197 all_vecs = basis_vectors + image_vectors
198 max_norm = max(np.linalg.norm(v) for v in all_vecs)
199 lim = max(1.2, 1.25 * max_norm)
200
201 for v, label in zip(basis_vectors, basis_labels):
202 ax.arrow(
203 0, 0, float(v[0]), float(v[1]),
204 length_includes_head=True,
205 head_width=0.06 * lim,
206 head_length=0.08 * lim,
207 linewidth=1.8,
208 color='black',
209 )
210 ax.text(float(v[0]) * 1.05, float(v[1]) * 1.05, label, fontsize=11)
211
212 for v, label in zip(image_vectors, image_labels):
213 ax.arrow(
214 0, 0, float(v[0]), float(v[1]),
215 length_includes_head=True,
216 head_width=0.03 * lim,
217 head_length=0.04 * lim,
218 linewidth=1.5,
219# linestyle="--",
220 color='red',
221 )
222 ax.text(float(v[0]) * 1.05, float(v[1]) * 1.05, label, fontsize=11)
223
224 ax.set_xlim(-lim, lim)
225 ax.set_ylim(-lim, lim)
226 ax.set_aspect("equal", adjustable="box")
227# ax.set_title(title)
228# ax.set_xlabel("x")
229# ax.set_ylabel("y")
230# ax.grid(True, alpha=0.3)
231
232 ax.set_xticks([])
233 ax.set_yticks([])
234 ax.set_xticklabels([])
235 ax.set_yticklabels([])
236
237 ax.grid(False)
238 ax.set_xlabel("")
239 ax.set_ylabel("")
240 ax.set_title("")
241
242 # 枠線（spines）を消す
243 for spine in ax.spines.values():
244 spine.set_visible(False)
245
246 # 軸線（axhline, axvline）を消したい場合はコメントアウト
247 # ax.axhline(0, color="0.75", lw=1)
248 # ax.axvline(0, color="0.75", lw=1)
249
250
251def main():
252 """
253 概要:
254 2x2実数行列の対角化を可視化し、結果を出力するメイン関数。
255
256 詳細説明:
257 コマンドライン引数から行列の要素を受け取り、以下の処理を実行する。
258 1. 元の基底ベクトルとその像を計算し、角度情報と共に表示する。
259 2. 行列が実数上で対角化可能かを判定し、その結果と理由を表示する。
260 3. 行列が対称であるかを判定する。
261 4. 実数上で対角化可能であれば、固有値と正規化された固有ベクトルを計算・表示し、
262 元の基底と固有ベクトルによる対角化基底の角度を比較する。
263 さらに、対角化された行列 `P^{-1} S P` を表示する。
264 5. 結果を2つのサブプロット（元の基底と対角化基底）または1つのサブプロット（対角化不可能時）に描画し、
265 指定されたファイル名で保存して表示する。
266
267 :returns: なし
268 """
269 parser = argparse.ArgumentParser(
270 description="Visualize diagonalization of a 2x2 real matrix."
271 )
272 parser.add_argument("S11", type=float)
273 parser.add_argument("S12", type=float)
274 parser.add_argument("S21", type=float)
275 parser.add_argument("S22", type=float)
276 parser.add_argument(
277 "--save",
278 type=str,
279 default="diagonalize2d.png",
280 help="Output figure filename (default: diagonalize2d.png)",
281 )
282 args = parser.parse_args()
283
284 S = np.array([[args.S11, args.S12], [args.S21, args.S22]], dtype=float)
285
286 a1 = np.array([1.0, 0.0])
287 a2 = np.array([0.0, 1.0])
288 Sa1 = S @ a1
289 Sa2 = S @ a2
290
291 print("=== Input matrix S ===")
292 print(S)
293 print()
294
295 print("=== Original basis ===")
296 print(f"a1 = {a1}, angle = {angle_deg(a1):.6f} deg")
297 print(f"a2 = {a2}, angle = {angle_deg(a2):.6f} deg")
298 print(f"angle(a1, a2) = {angle_between_deg(a1, a2):.6f} deg")
299 print()
300
301 print("=== Images of original basis ===")
302 print(f"S a1 = {Sa1}, parallel to a1? {is_parallel(a1, Sa1)}")
303 print(f"S a2 = {Sa2}, parallel to a2? {is_parallel(a2, Sa2)}")
304 print()
305
306 diag_ok, msg, vals, vecs = classify_real_diagonalizable(S)
307 print("=== Diagonalization test over R ===")
308 print(msg)
309 print()
310
311 symmetric = np.allclose(S, S.T, atol=1e-10)
312 print(f"Symmetric matrix? {symmetric}")
313 print()
314
315 if diag_ok:
316 vals, vecs = sort_eigensystem(vals.real, vecs.real)
317 v1 = vecs[:, 0]
318 v2 = vecs[:, 1]
319 Sv1 = S @ v1
320 Sv2 = S @ v2
321
322 print("=== Eigenvalues / eigenvectors ===")
323 print(f"lambda1 = {vals[0]:.12g}")
324 print(f"v1 = {v1}")
325 print(f"angle(v1) = {angle_deg(v1):.6f} deg")
326 print()
327 print(f"lambda2 = {vals[1]:.12g}")
328 print(f"v2 = {v2}")
329 print(f"angle(v2) = {angle_deg(v2):.6f} deg")
330 print()
331
332 ang_orig = angle_between_deg(a1, a2)
333 ang_diag = angle_between_deg(v1, v2)
334
335 print("=== Angle comparison ===")
336 print(f"angle(original basis) = {ang_orig:.6f} deg")
337 print(f"angle(diagonalizing basis) = {ang_diag:.6f} deg")
338 if symmetric:
339 print("For a symmetric matrix, this should stay 90 deg (orthogonal diagonalization).")
340 else:
341 print("For a non-symmetric matrix, the diagonalizing basis need not be orthogonal.")
342 print()
343
344 P = np.column_stack([v1, v2])
345 D = np.linalg.inv(P) @ S @ P
346 print("=== P^{-1} S P ===")
347 print(D)
348 print()
349
350 fig, axes = plt.subplots(1, 2, figsize=(11, 5))
351 ax_left, ax_right = axes
352
353 plot_vectors(
354 ax_left,
355 basis_vectors=[a1, a2],
356 image_vectors=[Sa1, Sa2],
357 basis_labels=["a1", "a2"],
358 image_labels=["S a1", "S a2"],
359 title="Original basis",
360 )
361 plot_vectors(
362 ax_right,
363 basis_vectors=[v1, v2],
364 image_vectors=[Sv1, Sv2],
365 basis_labels=["v1", "v2"],
366 image_labels=["S v1", "S v2"],
367 title="Diagonalizing basis (real eigenvectors)",
368 )
369 else:
370 print("=== Eigenvalues (possibly complex) ===")
371 print(vals)
372 print()
373
374 fig, ax = plt.subplots(1, 1, figsize=(5.5, 5))
375 plot_vectors(
376 ax,
377 basis_vectors=[a1, a2],
378 image_vectors=[Sa1, Sa2],
379 basis_labels=["a1", "a2"],
380 image_labels=["S a1", "S a2"],
381 title="Original basis (not diagonalizable over R)",
382 )
383
384 fig.suptitle("2D matrix visualization", fontsize=14)
385 fig.tight_layout()
386 fig.savefig(args.save, dpi=160)
387 print(f"Saved figure: {args.save}")
388 plt.show()
389
390
391if __name__ == "__main__":
392 main()