8. ジョルダン標準形とスペクトル集合
Pythonで学ぶ線形代数学
ジョルダン標準形
code: jordan.py
from sympy import *
from numpy.random import seed, permutation
A = diag(1, 2, 2, 2, 2, 3, 3, 3, 3, 3)
A1, 2 = A3, 4 = A5, 6 = A7, 8 = A8, 9 = 1
seed(123)
for n in range(10):
P = permutation(10)
for i, j in [(P2 * k, P2 * k + 1) for k in range(5)]:
A:, j += A:, i
Ai, : -= Aj, :
B = Lambda(S('lmd'), A - S('lmd') * eye(10))
x = Matrix(var('x0, x1, x2, x3, x4, x5, x6, x7, x8, x9'))
y = Matrix(var('y0, y1, y2, y3, y4, y5, y6, y7, y8, y9'))
z = Matrix(var('z0, z1, z2, z3, z4, z5, z6, z7, z8, z9'))
a1 = x.subs(solve(B(1) * x))
a2 = x.subs(solve(B(2) * x))
b2 = y.subs(solve(B(2) * y - a2))
a3 = x.subs(solve(B(3) * x))
b3 = y.subs(solve(B(3) * y - a3))
c3 = z.subs(solve(B(3) * z - b3))
v0 = a1.subs({x9: 1})
v1 = b2.subs({y6: 1, y7: 0, y8: 0, y9: 0})
v2 = B(2) * v1
v3 = b2.subs({y6: 0, y7: 1, y8: 0, y9: 0})
v4 = B(2) * v3
v5 = c3.subs({z5: 1, z6: 0, z7: 0, z8: 0, z9: 0})
v6 = B(3) * v5
v7 = B(3) * v6
v8 = b3.subs({y6: 1, y7: 0, y8: 0, y9: 0})
v9 = B(3) * v8
if __name__ == '__main__':
V = v0
for v in v1, v2, v3, v4, v5, v6, v7, v8, v9:
V = V.row_join(v)
print(repr(V.inv() * A * V))
ソースコード: jordan.py
実行結果:
code: python
Matrix([
35, 28, 24, 6, 26, 16, -6, 14, 26, 42,
-5, -8, -9, -24, -21, -2, -5, -3, -5, -17,
11, 14, 10, 5, 7, 8, 2, 6, 13, 13,
5, 2, 4, 2, 5, 1, -3, 1, 2, 7,
1, 4, 4, 6, 8, 0, -1, 2, 0, 5,
19, 11, 14, 19, 25, 13, 6, 7, 16, 31,
8, 11, 5, 7, 6, 7, 7, 5, 11, 10,
6, 16, 14, 40, 34, 2, 8, 7, 6, 26,
-27, -19, -16, -6, -19, -16, -2, -11, -22, -33,
-20, -21, -19, -11, -22, -9, 5, -10, -15, -28])
Matrix([
-1/3, 0, -1/3, 2, 47/6, -15/11, -6/11, 24/11, -1/2, 1,
5/3, 1/6, -1/6, 1/12, -7/12, 3/11, -10/11, 6/11, 0, 1/6,
-1, 5/6, 1/2, 5/12, 19/4, -1/11, -5/11, 0, 7/6, -1/6,
0, -1/2, -1/3, 5/4, 1/3, -8/11, 1/11, 6/11, -2/3, 1/3,
-1/3, 0, 0, -2, 3/2, 1, 0, 0, 1/6, 0,
8/3, -1, -1, -5/2, -2, 1, -3, 30/11, -1/3, 1,
-1/3, 1, 1/2, 0, 13/4, 0, -5/11, 0, 1, -1/6,
-8/3, 0, 1/3, 1, 7/6, 0, 20/11, -12/11, 0, -1/3,
-2, 0, 1/2, 0, -11/4, 0, 27/11, -30/11, 0, -1,
1, 0, 1/6, 0, -83/12, 0, 3/11, -12/11, 0, -1/2])
Matrix([
1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 2, 0, 0, 0, 0, 0, 0, 0, 0,
0, 1, 2, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 2, 0, 0, 0, 0, 0, 0,
0, 0, 0, 1, 2, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 3, 0, 0, 0, 0,
0, 0, 0, 0, 0, 1, 3, 0, 0, 0,
0, 0, 0, 0, 0, 0, 1, 3, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 3, 0,
0, 0, 0, 0, 0, 0, 0, 0, 1, 3])
固有値1の固有空間と一般化固有空間
https://gyazo.com/2d4cc316f1f848ae6231eb22c9e6b5ab
固有値2の固有空間と一般化固有空間
https://gyazo.com/5ca613b71641256c28d6b5d6149ea006
固有値3の固有空間と一般化固有空間
https://gyazo.com/9b81fed8e9a4ab5f529c9568a2063e01
ジョルダン標準形及びジョルダン分解の問題作成プログラム
code: jordan2.py
from sympy import Matrix, diag
from numpy.random import permutation
X = Matrix(1, 1, 0], 0, 1, 0, [0, 0, 2)
Y = Matrix(2, 1, 0], 0, 2, 1, [0, 0, 2)
Z = Matrix(2, 1, 0], 0, 2, 0, [0, 0, 2)
while True:
A = X.copy()
while 0 in A:
i, j, _ = permutation(3)
A:, j += A:, i
Ai, : -= Aj, :
if max(abs(A)) >= 10:
break
if max(abs(A)) < 10:
break
U, J = A.jordan_form()
print(f'A = {A}')
print(f'U = {U}')
print(f'U**(-1)*A*U = {J}')
C = U * diag(J0, 0, J1, 1, J2, 2) * U**(-1)
B = A - C
print(f'B = {B}')
print(f'C = {C}')
ソースコード: jordan2.py
実行例:
code: python
A = Matrix(2, 1, 1], -2, -2, -4, [1, 2, 4)
U = Matrix(1, 1, 0], -2, 0, -1, [1, 0, 1)
U**(-1)*A*U = Matrix(1, 1, 0], 0, 1, 0, [0, 0, 2)
B = Matrix(1, 1, 1], -2, -2, -2, [1, 1, 1)
C = Matrix(1, 0, 0], 0, 0, -2, [0, 1, 3)
行列のスペクトル集合
code: spectrum.py
from numpy import matrix, pi, sin, cos, linspace
from numpy.random import normal
from numpy.linalg import eig, eigh
import matplotlib.pylab as plt
N = 100
B = normal(0, 1, (N, N, 2))
A = matrix(B:, :, 0 + 1j * B:, :, 1)
Real = A + A.conj()
Hermite = A + A.H
PositiveDefinite = A * A.H
PositiveComponents = abs(A)
Unitary = matrix(eigh(Hermite)1)
X = PositiveComponents
Lmd = eig(X)0
r = max(abs(Lmd))
T = linspace(0, 2 * pi, 100)
plt.plot(r * cos(T), r * sin(T))
plt.scatter(Lmd.real, Lmd.imag, s=20)
plt.axis('equal'), plt.show()
ソースコード: spectrum.py
複素行列
https://gyazo.com/c151fce2f3971edaaaa1f169a710da4f
実行列
https://gyazo.com/d984bda740b6307aa359ef55ed3e8de4
エルミート行列
https://gyazo.com/82e2bf899628825d957c0c2949858047
正定値行列
https://gyazo.com/a9db4e3e99578704b6fe1b36426e0157
ユニタリ行列
https://gyazo.com/840a32bc03320cfab58f9569abc915c7
正成分行列
https://gyazo.com/c7dccc0c52d047875cd9971c548e5b1f
ノルム公式
code: norm.py
from numpy import matrix
from numpy.linalg import eig, norm
from numpy.random import normal
import matplotlib.pyplot as plt
def power(ax, m):
A = matrix(normal(0, 1, (m, m)))
lmd = max(abs(eig(A)0)) * 1.1
N = range(50)
P = (A / lmd)**n for n in N
for i in range(m):
for j in range(m):
ax.plot(N, [Bi, j for B in P])
ax.plot(N, norm(B, 2) for B in P)
fig, axs = plt.subplots(2, 5, figsize=(20, 5))
[power(axsij, 3) for i in range(2) for j in range(5)]
plt.show()
ソースコード: norm.py
https://gyazo.com/085359f1be61afdcd25697ac04745247
ゲルファントの公式
code: gelfand.py
from numpy import matrix
from numpy.linalg import eig, norm
from numpy.random import normal
import matplotlib.pyplot as plt
def gelfand(ax, m):
A = matrix(normal(0, 1, (m, m)))
lmd = max(abs(eig(A)0))
N = range(50)
P = A**n for n in N
ax.plot(N1:, [norm(Pn, 2)**(1 / n) for n in N1:])
ax.plot([N1, N-1], lmd, lmd)
fig, axs = plt.subplots(2, 5, figsize=(20, 5))
[gelfand(axsij, 3) for i in range(2) for j in range(5)]
plt.show()
ソースコード: gelfand.py
https://gyazo.com/48fdb811abb6b2229f0d2c31e7dd2626