線形代数で漸化式を解くやつ

2019-02-27

一部の漸化式は線形代数を用いて解くことができます。

教科書に載っているような正攻法と併せて見ていくことにしましょう。

備考

  • この記事の例題は引用元サイト様から拝借したので詳細が見たい方はそちらを参照ください。
  • 線形代数で解けるのは一次式で表せるような線形の漸化式です。 この記事で扱うのはこれに該当します。
目次

3項間漸化式

引用元:

隣接3項間型の漸化式 | 受験の月 (1)

まず普通に解きます。特性方程式と呼ばれるものを使います。

In [1]:
import sympy as sy
sy.init_printing()
$$ a_1 = 1 $$$$ a_2 = 3 $$$$ a_{n+2} - 2a_{n+1} - 8a_{n} = 0 $$

の一般項 $$ a_n $$ を求める

In [2]:
x, n, a_n, a_n1 = sy.symbols('x, n, a_n, a_n+1')
In [3]:
sy.solve(x ** 2 - 2 * x - 8)
Out[3]:
$$\left [ -2, \quad 4\right ]$$

特性方程式の解を用いて $$ \begin{eqnarray} \left\{ \begin{array}{l} a_{n+2} - 4a_{n+1} = -2(a_{n+1} - 4a_n) \\ a_{n+2} + 2a_{n+1} = 4(a_{n+1} + 2a_n) \end{array} \right. \end{eqnarray} $$

となる

$$ b_{n+1} = a_{n+1} - 4a_n $$$$ c_{n+1} = a_{n+1} + 2a_n $$

と置くと

$$ \begin{eqnarray} \left\{ \begin{array}{l} b_{n+1} = -2b_n \\ c_{n+1} = 4c_n \end{array} \right. \end{eqnarray} $$

となり、それぞれの公比は -2, -4 の等比数列と考えられる

公比を用いて n+1 のそれぞれの一般項を求める

$$ \begin{eqnarray} \left\{ \begin{array}{l} b_{n+1} = b_1 \cdot (-2)^n \\ c_{n+1} = c_1 \cdot 4^n \end{array} \right. \end{eqnarray} $$

b, c の第1項は求められないので、 第2項を求める

$$ \begin{eqnarray} \left\{ \begin{array}{l} b_2 = 3 - 4 = -1 \\ c_2 = 3 + 2 = 5 \end{array} \right. \end{eqnarray} $$

第2項を用いて n+1 の一般項を求める

$$ \begin{eqnarray} \left\{ \begin{array}{l} b_{n+1} = b_2 \cdot (-2)^{n-1} = -1 \cdot (-2)^{n-1} \\ c_{n+1} = c_2 \cdot 4^{n-1} = 5 \cdot 4^{n-1} \end{array} \right. \end{eqnarray} $$

先程の式と置き換えると

$$ \begin{eqnarray} \left\{ \begin{array}{l} a_{n+1} - 4a_n = -1 \cdot (-2)^{n-1} \\ a_{n+1} + 2a_n = 5 \cdot 4^{n-1} \end{array} \right. \end{eqnarray} $$

これらを連立して n の一般項を求める

In [4]:
ex1 = sy.Eq(a_n1 - 4 * a_n, -1 * (-2) ** (n-1))
ex2 = sy.Eq(a_n1 + 2 * a_n, 5 * 4 ** (n-1))
In [5]:
ex1
Out[5]:
$$- 4 a_{n} + a_{n+1} = - \left(-2\right)^{n - 1}$$
In [6]:
ex2
Out[6]:
$$2 a_{n} + a_{n+1} = 5 \cdot 4^{n - 1}$$
In [7]:
sy.solve([ex1, ex2], exclude=[n])
Out[7]:
$$\left \{ a_{n} : - \frac{\left(-2\right)^{n}}{12} + \frac{5 \cdot 4^{n}}{24}, \quad a_{n+1} : \frac{\left(-2\right)^{n}}{6} + \frac{5 \cdot 4^{n}}{6}\right \}$$
In [ ]:
 

なんかすっきりしませんね。 特性方程式って何なんでしょう。そもそもなぜクロスにかけてるんだろう?

という疑問を残しつつ、次は線形代数を使って解いてみます。

In [1]:
import numpy as np
import sympy as sy
sy.init_printing()
In [2]:
n = sy.Symbol('n')
$$ \left[ \begin{array}{rrr} a_{n+2} \\ a_{n+1} \end{array} \right] = \left[ \begin{array}{rrr} 2 & 8 \\ 1 & 0 \\ \end{array} \right] \left[ \begin{array}{rrr} a_{n+1} \\ a_{n} \end{array} \right] $$

を変形すると $$ \left[ \begin{array}{} a_{n+1} \\ a_{n} \end{array} \right] = \left[ \begin{array}{} 2 & 8 \\ 1 & 0 \\ \end{array} \right] ^ {n-1} \left[ \begin{array}{} a_{2} \\ a_{1} \end{array} \right] $$

In [3]:
A = sy.Matrix([[2, 8], [1, 0]])
A
Out[3]:
$$\left[\begin{matrix}2 & 8\\1 & 0\end{matrix}\right]$$
$$ \left[ \begin{array}{rrr} a_{n+1} \\ a_{n} \end{array} \right] = A^{n-1} \left[ \begin{array}{} a_{2} \\ a_{1} \end{array} \right] $$
In [4]:
A ** (n-1) * sy.Matrix([3, 1])
Out[4]:
$$\left[\begin{matrix}- \frac{\left(-2\right)^{n - 1}}{3} + \frac{10 \cdot 4^{n - 1}}{3}\\\frac{\left(-2\right)^{n - 1}}{6} + \frac{5 \cdot 4^{n - 1}}{6}\end{matrix}\right]$$

このように求めることができた(2行目が n の一般項)

でも、これでは面白くないし、人力では厳しいので対角化して解いてみることにする

In [5]:
# スペクトル分解して、固有値と固有ベクトルを求める
# 固有値, 重複度, 固有ベクトル(重複度分)
evs = A.eigenvects()
evs
Out[5]:
$$\left [ \left ( -2, \quad 1, \quad \left [ \left[\begin{matrix}-2\\1\end{matrix}\right]\right ]\right ), \quad \left ( 4, \quad 1, \quad \left [ \left[\begin{matrix}4\\1\end{matrix}\right]\right ]\right )\right ]$$
In [6]:
P = sy.Matrix([ev[-1][0].values() for ev in evs]).T
P
Out[6]:
$$\left[\begin{matrix}-2 & 4\\1 & 1\end{matrix}\right]$$
In [7]:
P ** -1
Out[7]:
$$\left[\begin{matrix}- \frac{1}{6} & \frac{2}{3}\\\frac{1}{6} & \frac{1}{3}\end{matrix}\right]$$
In [8]:
# 6倍したもの
P ** -1 * 6
Out[8]:
$$\left[\begin{matrix}-1 & 4\\1 & 2\end{matrix}\right]$$
$$ D = P^{-1} A P $$
In [9]:
D = P ** (-1) * A * P
D
Out[9]:
$$\left[\begin{matrix}-2 & 0\\0 & 4\end{matrix}\right]$$

D を用いて

$$ A^{n-1} = P D^{n-1} P^{-1} $$

のように表現できる

In [10]:
# これを計算すると同様の結果を得られる
(P) * (D ** (n-1)) * (P ** -1) * sy.Matrix([3, 1])
Out[10]:
$$\left[\begin{matrix}- \frac{\left(-2\right)^{n - 1}}{3} + \frac{10 \cdot 4^{n - 1}}{3}\\\frac{\left(-2\right)^{n - 1}}{6} + \frac{5 \cdot 4^{n - 1}}{6}\end{matrix}\right]$$

ちなみに $$ \left[ \begin{array}{rrr} a_{n+2} \\ a_{n+1} \end{array} \right]= A \left[ \begin{array}{rrr} a_{n+1} \\ a_{n} \end{array} \right] $$

を変形して

$$ \left[ \begin{array}{rrr} a_{n+2} \\ a_{n+1} \end{array} \right] = P D P^{-1} \left[ \begin{array}{rrr} a_{n+1} \\ a_{n} \end{array} \right] $$$$ P^{-1} \left[ \begin{array}{rrr} a_{n+2} \\ a_{n+1} \end{array} \right] = D P^{-1} \left[ \begin{array}{rrr} a_{n+1} \\ a_{n} \end{array} \right] $$
$$ \left[ \begin{array}{} -1 & 4 \\ 1 & 2 \\ \end{array} \right] \left[ \begin{array}{rrr} a_{n+2} \\ a_{n+1} \end{array} \right] = \left[ \begin{array}{} -2 & 0 \\ 0 & 4 \\ \end{array} \right] \left[ \begin{array}{} -1 & 4 \\ 1 & 2 \\ \end{array} \right] \left[ \begin{array}{rrr} a_{n+1} \\ a_{n} \end{array} \right] $$

(両辺を6倍してある)

漸化式の途中式と一致しているのがわかる $$ \begin{eqnarray} \left\{ \begin{array}{l} a_{n+2} - 4a_{n+1} = -2(a_{n+1} - 4a_n) \\ a_{n+2} + 2a_{n+1} = 4(a_{n+1} + 2a_n) \end{array} \right. \end{eqnarray} $$

In [ ]:
 

なるほど、特性方程式は固有方程式のことだったんですね。

そして固有値をクロスするように2式を組み立てたのは 固有ベクトルからなる行列の逆行列を掛けているのと同じことなんですね。 これが成り立つのは サイズが行列のサイズが 2x2 のときだけです。

(両辺にかけるので 1/6 は消えました)

重解を持つ3項間漸化式

引用元:

隣接3項間型の漸化式 | 受験の月 (3)

特性方程式により重解が得られる場合、それらの式を連立させることができませんが、 指数型の漸化式として解けます。( x^{n+1} で割るのがミソですね)

In [1]:
import numpy as np
import sympy as sy
sy.init_printing()
In [2]:
n, x = sy.symbols('n x')
$$ a_1 = 1 $$$$ a_2 = 5 $$$$ a_{n+2} = 6a_{n+1} - 9a_{n} $$
In [3]:
f = sy.Eq(x ** 2 - 6 * x + 9, 0)
f
Out[3]:
$$x^{2} - 6 x + 9 = 0$$
In [4]:
sy.solve(f)
Out[4]:
$$\left [ 3\right ]$$
$$ a_{n+2} - 3a_{n+1} = 3(a_{n+1} - 3a_{n}) $$
$$ b_{n+1} = a_{n+1} - 3a_{n} $$$$ b_{n+2} = 3b_{n+1} $$$$ b_{n+1} = b_1 \cdot 3^{n} $$
$$ b2 = 5 - 3 = 2 $$$$ b_{n+1} = b_2 \cdot 3^{n-1} $$$$ b_{n+1} = 2 \cdot 3^{n-1} $$
$$ a_{n+1} - 3a_{n} = 2 \cdot 3^{n-1} $$$$ a_{n+1} = 3a_{n} + 2 \cdot 3^{n-1} $$

ここで両辺を $ 3^{n+1} $ で割ると $$ \frac{a_{n+1}}{3^{n+1}} = \frac{a_{n}}{3^{n}} + 2 \cdot 3^{-2} $$

$$ b_n = \frac{a_n}{3^n} $$

とすると $$ b_1 = \frac{1}{3} $$ $$ b_{n+1} = b_{n} + \frac{2}{9} $$ $$ b_{n} = \frac{1}{3} + \frac{2}{9} \cdot (n-1) $$ $$ b_{n} = \frac{1}{3} + \frac{2}{9} n - \frac{2}{9} = \frac{2}{9} n + \frac{1}{9} $$ $$ \frac{a_n}{3^n} = \frac{2}{9} n + \frac{1}{9} = 2 \cdot 3^{-2}n + 3^{-2} $$ 両辺に $ 3^n $ をかけて $$ a_n = 2 \cdot 3^{n-2}n + 3^{n-2} $$ $$ a_n = 3^{n-2}(2n + 1) $$

In [ ]:
 

続いて線形代数を使って解きます。

独立した固有ベクトルが得られる場合、対角行列のn乗を使ってうまく計算できましたが、 今回の場合はそもそも対角行列を得られません。

ジョルダン標準形と呼ばれる形に整形された行列を使って計算したのが以下です。 対角行列ほどではありませんが、n乗の計算がしやすいです。

In [1]:
import numpy as np
import sympy as sy
sy.init_printing()
In [2]:
n = sy.Symbol('n')
$$ \left[ \begin{array}{rrr} a_{n+2} \\ a_{n+1} \end{array} \right] = \left[ \begin{array}{rrr} 6 & -9 \\ 1 & 0 \\ \end{array} \right] \left[ \begin{array}{rrr} a_{n+1} \\ a_{n} \end{array} \right] $$

を変形すると $$ \left[ \begin{array}{} a_{n+1} \\ a_{n} \end{array} \right] = \left[ \begin{array}{} 6 & -9 \\ 1 & 0 \\ \end{array} \right] ^ {n-1} \left[ \begin{array}{} a_{2} \\ a_{1} \end{array} \right] $$

In [3]:
A = sy.Matrix([[6, -9], [1, 0]])
A
Out[3]:
$$\left[\begin{matrix}6 & -9\\1 & 0\end{matrix}\right]$$
$$ \left[ \begin{array}{} a_{n+1} \\ a_{n} \end{array} \right] = A^{n-1} \left[ \begin{array}{} a_{2} \\ a_{1} \end{array} \right] $$
In [4]:
A ** (n-1) * sy.Matrix([5, 1])
Out[4]:
$$\left[\begin{matrix}6 \cdot 3^{n - 2} \left(n - 1\right) + 5 \cdot 3^{n - 1}\\2 \cdot 3^{n - 2} \left(n - 1\right) + 3^{n - 1}\end{matrix}\right]$$

先ほどと同様に求めることができた(2行目が n の一般項) やっぱり線形代数ってすごい

続いて固有値を求める

In [5]:
# 固有値と固有ベクトルを求めてみるが..
# 固有値, 重複度, 固有ベクトル(重複度分)
evs = A.eigenvects()
evs
Out[5]:
$$\left [ \left ( 3, \quad 2, \quad \left [ \left[\begin{matrix}3\\1\end{matrix}\right]\right ]\right )\right ]$$
In [6]:
sy.Matrix([[3, 3], [1, 1]]) ** -1

ValueErrorTraceback (most recent call last)
<ipython-input-6-1d1aea274b9b> in <module>
----> 1 sy.Matrix([[3, 3], [1, 1]]) ** -1

~/venv/lib/python3.7/site-packages/sympy/core/decorators.py in binary_op_wrapper(self, other)
    130                     else:
    131                         return f(self)
--> 132             return func(self, other)
    133         return binary_op_wrapper
    134     return priority_decorator

~/venv/lib/python3.7/site-packages/sympy/matrices/common.py in __pow__(self, num)
   2057                 if num < 0:
   2058                     num = -num
-> 2059                     a = a.inv()
   2060                 # When certain conditions are met,
   2061                 # Jordan block algorithm is faster than

~/venv/lib/python3.7/site-packages/sympy/matrices/matrices.py in inv(self, method, **kwargs)
   2916         if method is not None:
   2917             kwargs['method'] = method
-> 2918         return self._eval_inverse(**kwargs)
   2919 
   2920     def is_nilpotent(self):

~/venv/lib/python3.7/site-packages/sympy/matrices/dense.py in _eval_inverse(self, **kwargs)
    264         M = self.as_mutable()
    265         if method == "GE":
--> 266             rv = M.inverse_GE(iszerofunc=iszerofunc)
    267         elif method == "LU":
    268             rv = M.inverse_LU(iszerofunc=iszerofunc)

~/venv/lib/python3.7/site-packages/sympy/matrices/matrices.py in inverse_GE(self, iszerofunc)
   2831         red = big.rref(iszerofunc=iszerofunc, simplify=True)[0]
   2832         if any(iszerofunc(red[j, j]) for j in range(red.rows)):
-> 2833             raise ValueError("Matrix det == 0; not invertible.")
   2834 
   2835         return self._new(red[:, big.rows:])

ValueError: Matrix det == 0; not invertible.

固有ベクトルが重複していて対角化ができない(逆行列が求められない)ため

ジョルダン標準化して求めてみる

In [7]:
P, J = A.jordan_form()
In [8]:
P, J
Out[8]:
$$\left ( \left[\begin{matrix}3 & 1\\1 & 0\end{matrix}\right], \quad \left[\begin{matrix}3 & 1\\0 & 3\end{matrix}\right]\right )$$
$$ J = P^{-1} A P $$$$ A = P J P^{-1} $$
In [9]:
A * P == P * J
Out[9]:
True
$$ A^{n} = P J^{n} P^{-1} $$
$$ \left[ \begin{array}{} a_{n+1} \\ a_{n} \end{array} \right] = P J^{n-1} P^{-1} \left[ \begin{array}{} a_{2} \\ a_{1} \end{array} \right] $$
In [10]:
P ** -1
Out[10]:
$$\left[\begin{matrix}0 & 1\\1 & -3\end{matrix}\right]$$
In [11]:
J ** (n-1)
Out[11]:
$$\left[\begin{matrix}3^{n - 1} & 3^{n - 2} \left(n - 1\right)\\0 & 3^{n - 1}\end{matrix}\right]$$
$$ \left[ \begin{array}{} a_{n+1} \\ a_{n} \end{array} \right] = \left[ \begin{array}{} 3 & 1 \\ 1 & 0 \end{array} \right] \left[ \begin{array}{} 3^{n-1} & 3^{n-2}(n-1) \\ 0 & 3^{n-1} \end{array} \right] \left[ \begin{array}{} 0 & 1 \\ 1 & -3 \end{array} \right] \left[ \begin{array}{} 5 \\ 1 \end{array} \right] $$
In [12]:
P * (J ** (n-1)) * (P ** -1) * sy.Matrix([5, 1])
Out[12]:
$$\left[\begin{matrix}6 \cdot 3^{n - 2} \left(n - 1\right) + 5 \cdot 3^{n - 1}\\2 \cdot 3^{n - 2} \left(n - 1\right) + 3^{n - 1}\end{matrix}\right]$$
In [ ]:
 

ジョルダン標準形に分解した行列は指数型の漸化式との関係を見いだせませんでしたが、 これならなんとか人力で解けそうです。

2項間漸化式

引用元:

特殊解型の漸化式 | 受験の月

難しい方の3項間の漸化式からやりましたが、 2項間であってももちろん解けます。

まず普通の解き方から

In [1]:
import sympy as sy
sy.init_printing()
$$ a_{1} = 5 $$$$ a_{n+1} = 3a_n + 4 $$

の一般項 $$ a_n $$ を求める

In [2]:
x, n, a_n, a_n1 = sy.symbols('x, n, a_n, a_n+1')
$$ a_{n+1} = y , a_{n} = x $$

とおいて $$ y = x $$ との交点を求めるために y に x を代入して $$ x = 3x + 4 $$ を解く

In [3]:
sy.solve(sy.Eq(x, 3 * x + 4))
Out[3]:
$$\left [ -2\right ]$$

解を用いて $$ y + 2 = 3(x + 2) $$ $$ a_{n+1} + 2 = 3(a_{n} + 2) $$

$$ b_{n} = a_{n} + 2 $$

と置いて

$$ b_{n+1} = 3b_{n} $$$$ b_{n} = b_{1}3^{n-1} $$
$$ b_1 = a_{1} + 2 = 5 + 2 = 7 $$

なので $$ b_{n} = 7 \cdot 3^{n-1} $$ $$ a_{n} + 2 = 7 \cdot 3^{n-1} $$ $$ a_{n} = 7 \cdot 3^{n-1} - 2 $$

In [ ]:
 

n, n+1 を x と置きます。 これによって y = x の直線との交点を求められるので 両辺の項を整えて、公比数列型の漸化式にします。

(これも特性方程式というんだけど、3項間漸化式の特性方程式とは全然違うのでなんか気持ち悪い)

これも行列を使って解けます

In [1]:
import sympy as sy
sy.init_printing()
In [2]:
n = sy.Symbol('n')
$$ a_{1} = 5 $$$$ a_{n+1} = 3a_n + 4 $$

における $$ a_n $$ の一般項を求める

$$ \left[ \begin{array}{rrr} a_{n+1} \\ 1 \end{array} \right] = \left[ \begin{array}{rrr} 3 & 4 \\ 0 & 1 \\ \end{array} \right] \left[ \begin{array}{rrr} a_{n} \\ 1 \end{array} \right] $$

を変形すると $$ \left[ \begin{array}{} a_{n} \\ 1 \end{array} \right] = \left[ \begin{array}{} 3 & 4 \\ 0 & 1 \\ \end{array} \right] ^ {n-1} \left[ \begin{array}{} 5 \\ 1 \end{array} \right] $$

In [3]:
A = sy.Matrix([[3, 4], [0, 1]])
A
Out[3]:
$$\left[\begin{matrix}3 & 4\\0 & 1\end{matrix}\right]$$
In [4]:
A ** (n-1) * sy.Matrix([5, 1])
Out[4]:
$$\left[\begin{matrix}7 \cdot 3^{n - 1} - 2\\1\end{matrix}\right]$$

1行目に n の一般項を求めることができた

これもやはり計算が難しいので固有値分解していく

In [5]:
evs = A.eigenvects()
evs
Out[5]:
$$\left [ \left ( 1, \quad 1, \quad \left [ \left[\begin{matrix}-2\\1\end{matrix}\right]\right ]\right ), \quad \left ( 3, \quad 1, \quad \left [ \left[\begin{matrix}1\\0\end{matrix}\right]\right ]\right )\right ]$$
In [6]:
P = sy.Matrix([[-2, 1], [1, 0]])
P
Out[6]:
$$\left[\begin{matrix}-2 & 1\\1 & 0\end{matrix}\right]$$
In [7]:
P ** -1
Out[7]:
$$\left[\begin{matrix}0 & 1\\1 & 2\end{matrix}\right]$$
$$ D = P^{-1} A P $$
In [8]:
D = P ** (-1) * A * P
D
Out[8]:
$$\left[\begin{matrix}1 & 0\\0 & 3\end{matrix}\right]$$

D を用いて

$$ A^{n-1} = P D^{n-1} P^{-1} $$

なので

$$ \left[ \begin{array}{} a_{n} \\ 1 \end{array} \right] = A ^ {n-1} \left[ \begin{array}{} 5 \\ 1 \end{array} \right] $$

$$ \left[ \begin{array}{} a_{n} \\ 1 \end{array} \right] = P D^{n-1} P^{-1} \left[ \begin{array}{} 5 \\ 1 \end{array} \right] $$

となり

In [9]:
# これを計算すると同様の結果を得られる
(P) * (D ** (n-1)) * (P ** -1) * sy.Matrix([5, 1])
Out[9]:
$$\left[\begin{matrix}7 \cdot 3^{n - 1} - 2\\1\end{matrix}\right]$$

ちなみに $$ P^{-1} \left[ \begin{array}{} a_{n+1} \\ 1 \end{array} \right] = D P^{-1} \left[ \begin{array}{} a_{n} \\ 1 \end{array} \right] $$ なので

$$ \left[ \begin{array}{rrr} 0 & 1 \\ 1 & 2 \\ \end{array} \right] \left[ \begin{array}{} a_{n+1} \\ 1 \end{array} \right] = \left[ \begin{array}{rrr} 1 & 0 \\ 0 & 3 \\ \end{array} \right] \left[ \begin{array}{rrr} 0 & 1 \\ 1 & 2 \\ \end{array} \right] \left[ \begin{array}{} a_{n} \\ 1 \end{array} \right] $$

2行目が漸化式の途中式と一致しているのがわかる (1行目は捨てていい)

$$ a_{n+1} + 2 = 3(a_n + 2) $$

実際に計算すると $$ \left[ \begin{array}{rrr} 1 \\ a_n + 2 \\ \end{array} \right] = \left[ \begin{array}{} 1 & 0 \\ 0 & 3^{n-1} \\ \end{array} \right] \left[ \begin{array}{} 1 \\ 7 \\ \end{array} \right] $$ $$ a_{n} + 2 = 7 \cdot 3^{n-1} $$ $$ a_{n} = 7 \cdot 3^{n-1} -2 $$

In [ ]:
 

3項間漸化式と比べて項が足りないので、 左辺の下行は 1 固定で計算しました。

連立漸化式

引用元:

連立漸化式 | 受験の月

2つの異なる数列の一般項がお互いを参照しあっているような漸化式です。 説明しづらいので詳しくは問題を読んでください。

両辺を揃えて等比数列型の漸化式に帰着させて解くのが一般的です。

In [1]:
import sympy as sy
sy.init_printing()
In [2]:
an, bn, n = sy.symbols('a_n, b_n, n')
$$ a_1 = 1 $$$$ b_1 = 2 $$$$ \begin{eqnarray} \left\{ \begin{array}{l} a_{n+1} = 2a_n + b_n \\ b_{n+1} = 4a_n - b_n \end{array} \right. \end{eqnarray} $$

における ab の一般項 をもとめる

方針: 右辺をβ倍したら左辺になるような等比数列を作る $$ a_{n+1} + αb_{n+1} = β(a_n + αb_n) $$ $$ a_{n+1} + αb_{n+1} = βa_n + αβb_n $$

(左辺が成り立つように)下式を α倍 して上式に加える $$ a_{n+1} + αb_{n+1} = 2a_n + b_n + α(4a_n - b_n) $$ $$ a_{n+1} + αb_{n+1} = 2a_n + b_n + 4αa_n - αb_n = a_n(2 + 4α) + b_n(1 - α)$$ よって $$ \begin{eqnarray} \left\{ \begin{array}{l} 2 + 4α = β \\ 1 - α = αβ \end{array} \right. \end{eqnarray} $$ 連立して解くと $$ (α, β) = (-1, -2), (\frac{1}{4}, 3) $$ となる

これを用い以下のように式を組み立てる $$ \begin{eqnarray} \left\{ \begin{array}{l} a_{n+1} - b_{n+1} = -2(a_n - b_n) \\ a_{n+1} + \frac{1}{4} b_{n+1} = 3(a_n + \frac{1}{4} b_n) \end{array} \right. \end{eqnarray} $$

$$ c_n = a_n - b_n $$

とおいて $$ c_1 = 1 - 2 = -1 $$ $$ c_{n} = -(-2)^{n-1} $$

$$ d_n = a_n + \frac{1}{4} b_n $$

とおいて $$ d_1 = 1 + \frac{2}{4} = \frac{3}{2} $$ $$ d_n = \frac{3^n}{2} $$

$$ \begin{eqnarray} \left\{ \begin{array}{l} a_n - b_n = -(-2)^{n-1} \\ a_n + \frac{1}{4} b_n = \frac{3^n}{2} \end{array} \right. \end{eqnarray} $$

を連立して解くと

In [3]:
ex1 = sy.Eq(an - bn, -(-2) ** (n-1))
ex2 = sy.Eq(an + bn / 4, (3 ** n) / 2)
In [4]:
sy.solve([ex1, ex2], exclude=[n])
Out[4]:
$$\left \{ a_{n} : \frac{\left(-2\right)^{n}}{10} + \frac{2 \cdot 3^{n}}{5}, \quad b_{n} : - \frac{2 \left(-2\right)^{n}}{5} + \frac{2 \cdot 3^{n}}{5}\right \}$$
In [ ]:
 

続いて行列を使ってときます。とっても簡単。

In [1]:
import sympy as sy
sy.init_printing()
In [2]:
n = sy.Symbol('n')
$$ a_1 = 1 $$$$ b_1 = 2 $$$$ \begin{eqnarray} \left\{ \begin{array}{l} a_{n+1} = 2a_n + b_n \\ b_{n+1} = 4a_n - b_n \end{array} \right. \end{eqnarray} $$

における ab の一般項 をもとめる

今までのやり方に沿って変換行列を考えてみる $$ \left[ \begin{array}{rrr} a_{n+1} \\ b_{n+1} \\ \end{array} \right] = \left[ \begin{array}{rrr} 2 & 1 \\ 4 & -1 \\ \end{array} \right] \left[ \begin{array}{rrr} a_{n} \\ b_{n} \\ \end{array} \right] $$ を変形すると $$ \left[ \begin{array}{rrr} a_{n+1} \\ b_{n+1} \\ \end{array} \right] = \left[ \begin{array}{rrr} 2 & 1 \\ 4 & -1 \\ \end{array} \right]^{n} \left[ \begin{array}{rrr} a_{1} \\ b_{1} \\ \end{array} \right] $$

求めたいのは n なので $$ \left[ \begin{array}{rrr} a_{n} \\ b_{n} \\ \end{array} \right] = \left[ \begin{array}{rrr} 2 & 1 \\ 4 & -1 \\ \end{array} \right]^{n-1} \left[ \begin{array}{rrr} a_{1} \\ b_{1} \\ \end{array} \right] $$

In [3]:
A = sy.Matrix([
    [2, 1],
    [4, -1],
])
A
Out[3]:
$$\left[\begin{matrix}2 & 1\\4 & -1\end{matrix}\right]$$
In [4]:
# 実際にこれを計算してみると
A ** (n-1) * sy.Matrix([1, 2])
Out[4]:
$$\left[\begin{matrix}- \frac{\left(-2\right)^{n - 1}}{5} + \frac{6 \cdot 3^{n - 1}}{5}\\\frac{4 \left(-2\right)^{n - 1}}{5} + \frac{6 \cdot 3^{n - 1}}{5}\end{matrix}\right]$$
In [5]:
evs = A.eigenvects()
evs
Out[5]:
$$\left [ \left ( -2, \quad 1, \quad \left [ \left[\begin{matrix}- \frac{1}{4}\\1\end{matrix}\right]\right ]\right ), \quad \left ( 3, \quad 1, \quad \left [ \left[\begin{matrix}1\\1\end{matrix}\right]\right ]\right )\right ]$$
In [6]:
P = sy.Matrix([ev[-1][0].values() for ev in evs]).T
P
Out[6]:
$$\left[\begin{matrix}- \frac{1}{4} & 1\\1 & 1\end{matrix}\right]$$
In [7]:
P ** -1
Out[7]:
$$\left[\begin{matrix}- \frac{4}{5} & \frac{4}{5}\\\frac{4}{5} & \frac{1}{5}\end{matrix}\right]$$
In [8]:
D = P ** (-1) * A * P
D
Out[8]:
$$\left[\begin{matrix}-2 & 0\\0 & 3\end{matrix}\right]$$
$$ \left[ \begin{array}{rrr} a_{n+1} \\ b_{n+1} \\ \end{array} \right] = A \left[ \begin{array}{rrr} a_{n} \\ b_{n} \\ \end{array} \right] $$$$ P^{-1} \left[ \begin{array}{rrr} a_{n+1} \\ b_{n+1} \\ \end{array} \right] = D P^{-1} \left[ \begin{array}{rrr} a_{n} \\ b_{n} \\ \end{array} \right] $$
In [9]:
# 両辺にかけるので分数を取る
P ** -1 * 5
Out[9]:
$$\left[\begin{matrix}-4 & 4\\4 & 1\end{matrix}\right]$$
$$ \left[ \begin{array}{rrr} -4 & 4 \\ 4 & 1 \\ \end{array} \right] \left[ \begin{array}{rrr} a_{n+1} \\ b_{n+1} \\ \end{array} \right] = \left[ \begin{array}{rrr} -2 & 0 \\ 0 & 3 \\ \end{array} \right] \left[ \begin{array}{rrr} -4 & 4 \\ 4 & 1 \\ \end{array} \right] \left[ \begin{array}{rrr} a_{n} \\ b_{n} \\ \end{array} \right] $$

組み立てると先ほどの連立方程式の途中式と一致する $$ \begin{eqnarray} \left\{ \begin{array}{l} -4a_{n+1} + 4b_{n+1} = -2(-4a_n + 4b_n) \\ 4a_{n+1} + b_{n+1} = 3(4a_n + b_n) \end{array} \right. \end{eqnarray} $$

In [ ]:
 
参考