【1. 問題設定】
与えられた連立差分方程式：
x_{n+1} = 2x_n + y_n （初期値 x_0 = 7）
y_{n+1} = 2x_n + 3y_n （初期値 y_0 = -4）
行列表現：
[ x_{n+1} ] = [ 2 1 ] [ x_n ]
[ y_{n+1} ] [ 2 3 ] [ y_n ]
すなわち
x_{n+1} = A * x_n
（A は 2×2 の係数行列）
A = [[2, 1], [2, 3]]
x_0 = [7, -4]^T
───────────────────────────────
【2. z変換による式変形】
z変換の基本式：
zX(z) - zx_0 = Z{x_{n+1}} - zx_0
同様に Y(z) にも適用して、z領域での式を立てる：
[ (z - 2) (-1) ] [ X(z) ] = [ 7z ]
[ (-2) (z - 3) ] [ Y(z) ] [ 4z ]
すなわち
A_z * X(z) = B(z)
───────────────────────────────
【3. 行列方程式を解く】
A_z の行列式：
|A_z| = (z - 2)(z - 3) - 2
= z^2 - 5z + 4
= (z - 1)(z - 4)
逆行列：
A_z^{-1} = (1 / ((z - 1)(z - 4))) * [[z - 3, 1], [2, z - 2]]
したがって：
[ X(z) ] = (1 / ((z - 1)(z - 4))) * [[z - 3, 1], [2, z - 2]] * [7z, 4z]^T
───────────────────────────────
【4. 展開】
X(z) = (7z^2 - 25z) / ((z - 1)(z - 4))
Y(z) = (-4z^2 + 22z) / ((z - 1)(z - 4))
───────────────────────────────
【5. 逆z変換と部分分数分解】
分母 (z - 1)(z - 4) は、A の特性方程式
λ^2 - 5λ + 4 = 0 に対応。
固有値 λ1 = 1, λ2 = 4。
(1) X(z) の部分分数分解：
X(z) = Az/(z - 1) + Bz/(z - 4)
両辺に (z - 1)(z - 4) を掛けて比較：
7z^2 - 25z = Az(z - 4) + Bz(z - 1)
z = 1 のとき → 7 - 25 = A*(1)(-3) → A = 6
z = 4 のとき → 7(16) - 25*(4) = B*(4)*(3) → B = 7
よって
X(z) = 6z/(z - 1) + 7z/(z - 4)
(2) Y(z) の部分分数分解：
Y(z) = Cz/(z - 1) + Dz/(z - 4)
両辺に (z - 1)(z - 4) を掛けて比較：
-4z^2 + 22z = Cz(z - 4) + Dz(z - 1)
z = 1 のとき → -4 + 22 = C*(1)(-3) → C = -6
z = 4 のとき → -4(16) + 22*(4) = D*(4)*(3) → D = 2
よって
Y(z) = (-6z)/(z - 1) + (2z)/(z - 4)
───────────────────────────────
【6. 逆z変換】
基本公式：
z / (z - a) ⇔ a^n
したがって
x_n = 6*(1)^n + 7*(4)^n = 74^n + 6
y_n = (-6)(1)^n + 2*(4)^n = 2*4^n - 6
───────────────────────────────
【7. 結果】
x_n = 74^n + 6
y_n = 24^n - 6
───────────────────────────────
【8. 固有値との対応】
状態方程式：
[ x_{n+1} ] = [ 2 1 ] [ x_n ]
[ y_{n+1} ] [ 2 3 ] [ y_n ]
行列 A の固有値は：
det(A - λI) = |2 - λ, 1; 2, 3 - λ|
= (2 - λ)(3 - λ) - 2
= λ^2 - 5λ + 4
→ λ1 = 1, λ2 = 4
z領域の分母 (z - 1)(z - 4) と完全一致。
───────────────────────────────
【9. 図と数式の対応】
Step1: 差分方程式を立てる
Step2: (zI - A) * X(z) = z * x_0
Step3: 逆行列 A_z^{-1} で X(z) を解く
Step4: 部分分数分解して逆z変換、x_n, y_n を求める
───────────────────────────────
【10. まとめ】
z変換と固有値解析は等価構造を持つ。
x_{n+1} = A*x_n
⇔ (zI - A)X(z) = zx_0
固有値 λ_i は z変換の極 z = λ_i に対応し、
一般解は x_n = Σ C_i * λ_i^n の形で表される。
したがって、
z変換で解くことは「固有値展開による時間応答解析」に相当する。
───────────────────────────────
【最終結果】
x_n = 74^n + 6
y_n = 24^n - 6
import numpy as np
import matplotlib.pyplot as plt
# ==============================
# パラメータ定義 / Parameters
# ==============================
A = np.array([[2, 1],
[2, 3]]) # 係数行列 / Coefficient matrix
x0 = np.array([7, -4]) # 初期値 / Initial vector
n_max = 10 # 計算ステップ数 / Number of steps
# ==============================
# 差分方程式を逐次計算 / Iterative calculation of x_{n+1} = A * x_n
# ==============================
x_vals = [x0]
for n in range(n_max):
x_next = A @ x_vals[-1]
x_vals.append(x_next)
x_vals = np.array(x_vals)
# ==============================
# 固有値・固有ベクトル解析 / Eigenvalue and eigenvector analysis
# ==============================
eigvals, eigvecs = np.linalg.eig(A)
# ==============================
# 理論解（固有展開による解析解）/ Analytical solution via eigen decomposition
# ==============================
# 一般解: x_n = c1 * λ1^n * v1 + c2 * λ2^n * v2
v1, v2 = eigvecs[:,0], eigvecs[:,1]
λ1, λ2 = eigvals[0], eigvals[1]
# 初期条件から c1, c2 を決定
V = np.column_stack((v1, v2))
c = np.linalg.inv(V) @ x0
n = np.arange(n_max + 1)
x_theoretical = np.array([c[0]*λ1**n_i*v1 + c[1]*λ2**n_i*v2 for n_i in n])
# ==============================
# 結果表示 / Display results
# ==============================
print("Eigenvalues:", eigvals)
print("Eigenvectors:\n", eigvecs)
print("Analytical formula:")
print("x_n = 7*4^n + 6, y_n = 2*4^n - 6")
# ==============================
# グラフ描画 / Plot results
# ==============================
plt.figure(figsize=(8,5))
plt.plot(n, x_vals[:,0], 'o-', label='x_n (iterative)')
plt.plot(n, x_vals[:,1], 's-', label='y_n (iterative)')
plt.plot(n, x_theoretical[:,0], 'r--', label='x_n (theoretical)')
plt.plot(n, x_theoretical[:,1], 'g--', label='y_n (theoretical)')
plt.title("Difference Equation Solution (Eigenvalue Method)", fontsize=12)
plt.xlabel("n (step index)")
plt.ylabel("Value")
plt.legend()
plt.grid(True)
plt.show()