【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()