/linalg/8. ジョルダン標準形とスペクトル集合

generated at 2/17/2025, 12:24:49 PM
8. ジョルダン標準形とスペクトル集合
Pythonで学ぶ線形代数学

ジョルダン標準形
 jordan.pyfrom sympy import *
from numpy.random import seed, permutation

A = diag(1, 2, 2, 2, 2, 3, 3, 3, 3, 3)
A[1, 2] = A[3, 4] = A[5, 6] = A[7, 8] = A[8, 9] = 1

seed(123)
for n in range(10):
    P = permutation(10)
    for i, j in [(P[2 * k], P[2 * k + 1]) for k in range(5)]:
        A[:, j] += A[:, i]
        A[i, :] -= A[j, :]

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

実行結果:
 pythonMatrix([
[ 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の固有空間と一般化固有空間

固有値2の固有空間と一般化固有空間

固有値3の固有空間と一般化固有空間



ジョルダン標準形及びジョルダン分解の問題作成プログラム
 jordan2.pyfrom 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]
        A[i, :] -= A[j, :]
        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(J[0, 0], J[1, 1], J[2, 2]) * U**(-1)
B = A - C
print(f'B = {B}')
print(f'C = {C}')
ソースコード: jordan2.py

実行例:
 pythonA = 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]])


行列のスペクトル集合
 spectrum.pyfrom 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

複素行列

実行列

エルミート行列

正定値行列

ユニタリ行列

正成分行列


ノルム公式
 norm.pyfrom 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, [B[i, 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(axs[i][j], 3) for i in range(2) for j in range(5)]
plt.show()
ソースコード: norm.py


ゲルファントの公式
 gelfand.pyfrom 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(N[1:], [norm(P[n], 2)**(1 / n) for n in N[1:]])
    ax.plot([N[1], N[-1]], [lmd, lmd])

fig, axs = plt.subplots(2, 5, figsize=(20, 5))
[gelfand(axs[i][j], 3) for i in range(2) for j in range(5)]
plt.show()
ソースコード: gelfand.py