行列の対角化関係の手計算を検算

20230218_matrix-diagonalization.py

import numpy as np

A = np.matrix([
	[3, 0, 2],
	[-3, 2, 5],
	[4, 0, 1]
])
P = np.matrix([
	[3, 0, 3],
	[13, 1, 2],
	[-6, 0, 3]
])
D = np.matrix([
	[-1, 0, 0],
	[0, 2, 0],
	[0, 0, 5]
])
P_inv = P.I

# 逆行列
print(P_inv*9)
"""
[[  1.   0.  -1.]
 [-17.   9.  11.]
 [  2.   0.   1.]]
"""

# 対角化
print(P_inv.dot(A).dot(P))
"""
[[-1.00000000e+00  0.00000000e+00  2.49800181e-16]
 [ 0.00000000e+00  2.00000000e+00  1.77635684e-15]
 [ 0.00000000e+00  0.00000000e+00  5.00000000e+00]]
"""

# A^k
def check_A_power_k(k):
	# A^k, D^k 計算：愚直に掛け算
	A_k = A
	D_k = D
	for _ in range(k-1):
		A_k = A_k.dot(A)
		D_k = D_k.dot(D)
	# A^k を計算：P D^k P^{-1}
	PDkPinv_1 = P.dot(D_k).dot(P_inv)
	# A^k を計算：手計算の最終結果にkを代入
	l1_k = (-1)**k
	l2_k = 2**k
	l3_k = 5**k
	PDkPinv_2 = np.matrix([
		[3*l1_k+6*l3_k, 0, -3*l1_k+3*l3_k],
		[13*l1_k-17*l2_k+4*l3_k, 9*l2_k, -13*l1_k+11*l2_k+2*l3_k],
		[-6*l1_k+6*l3_k, 0, 6*l1_k+3*l3_k]
	]) / 9.0
	print('愚直に計算した A^{}'.format(k))
	print(A_k)
	print('A^{k} が機械的に行列演算した P D^{k} P^-1 に一致するか？'.format(k=k))
	print(np.abs(A_k - PDkPinv_1) < 1e-10)
	print('A^{k} が要素を手計算した P D^{k} P^-1 の式に一致するか？'.format(k=k))
	print(np.abs(A_k - PDkPinv_2) < 1e-10)

check_A_power_k(1)
"""
愚直に計算した A^1
[[ 3  0  2]
 [-3  2  5]
 [ 4  0  1]]
A^1 が機械的に行列演算した P D^1 P^-1 に一致するか？
[[ True  True  True]
 [ True  True  True]
 [ True  True  True]]
A^1 が要素を手計算した P D^1 P^-1 の式に一致するか？
[[ True  True  True]
 [ True  True  True]
 [ True  True  True]]
"""
check_A_power_k(3)
"""
愚直に計算した A^3
[[83  0 42]
 [39  8 39]
 [84  0 41]]
A^3 が機械的に行列演算した P D^3 P^-1 に一致するか？
[[ True  True  True]
 [ True  True  True]
 [ True  True  True]]
A^3 が要素を手計算した P D^3 P^-1 の式に一致するか？
[[ True  True  True]
 [ True  True  True]
 [ True  True  True]]
"""
check_A_power_k(10)
"""
愚直に計算した A^10
[[6510417       0 3255208]
 [4338345    1024 2171389]
 [6510416       0 3255209]]
A^10 が機械的に行列演算した P D^10 P^-1 に一致するか？
[[ True  True  True]
 [ True  True  True]
 [ True  True  True]]
A^10 が要素を手計算した P D^10 P^-1 の式に一致するか？
[[ True  True  True]
 [ True  True  True]
 [ True  True  True]]
"""