例えば、次のような斉次線形漸化式
an+2 - 5an+1 + 6an = 0
には、次の方程式が対応しています。
x2 - 5x + 6 = 0
なぜなら、上の漸化式を次の形にしようとすると、
an+2 - αan+1 = β(an+1 - αan)
これを展開して、
an+2 - (α + β)an+1 + αβan = 0
だから、
α + β = 5
αβ = 6
となり、上の方程式の解になっています。これは4項以上の漸化式でも同じです。
ここで、次の漸化式を解いてみましょう。
an+3 - 9an+2 + 26an+1 - 24an = 0
対応する方程式
x3 - 9x2 + 26x - 24 = 0
の解は、x = 2, 3, 4なので、
an+3 - 7an+2 + 12an+1 = 2(an+2 - 7an+1 + 12an)
an+3 - 6an+2 + 8an+1 = 3(an+2 - 6an+1 + 8an)
an+3 - 5an+2 + 6an+1 = 4(an+2 - 5an+1 + 6an)
とりあえずこれを再帰的に解いて、
an+2 - 7an+1 + 12an = -5 • 2n
an+2 - 6an+1 + 8an = -4 • 3n
an+2 - 5an+1 + 6an = -3 • 4n
あとは連立方程式を解くだけです。
an = -3/2 • 4n + 4 • 3n -5/2 • 2n
重解があるときは少しややこしくなりますが、そうでないときはこのように簡単に解けます。
上の手順をコードにしてみました。
from fractions import Fractiondef mat_inv(M):
n = len(M)
A = [ M[i] + [ int(i==k) for k in range(n) ] for i in range(n) ]
for i in range(n):
for k in range(i + 1, n * 2):
A[i][k] = A[i][k] / A[i][i]
for j in range(i + 1, n):
c = A[j][i]
for k in range(i + 1, n * 2):
A[j][k] -= A[i][k] * c
for i in range(n - 1, 0, -1):
for j in range(i):
for k in range(n, n * 2):
A[j][k] -= A[j][i] * A[i][k]
return [ A[i][n:] for i in range(n) ]def gen_divisors(n):
for m in filter(lambda m: n % m == 0, range(1, n + 1)):
yield mdef gen_solution_candidates(a):
for m in gen_divisors(abs(a[0])):
for n in gen_divisors(abs(a[-1])):
yield Fraction(n, m)
yield -Fraction(n, m)def mul(f, g):
h = [ 0 ] * (len(f) + len(g) - 1)
for k in range(len(f)):
for l in range(len(g)):
h[k+l] += f[k] * g[l]
return hdef inner_product(a, b):
return sum(map(lambda x: x[0] * x[1], zip(a, b)))def value(a, x):
return reduce(lambda p, q: p * x + q, a)def solve_eq(a):
return filter(lambda x: value(a, x) == 0, gen_solution_candidates(a))def print_solution(s, c, b):
str_sol = "a_n ="
first = True
for k in range(len(s)):
a = c[k] * b[k]
m = s[k]
if a == 0:
continue
if first:
str_sol += " " + str(a)
else:
if a > 0:
str_sol += " + " + str(a)
else:
str_sol += " - " + str(-a)
if m > 1 and m.denominator == 1:
str_sol += " * " + str(m) + "^n"
elif m != 1:
str_sol += " * (" + str(m) + ")^n"
first = False
print str_soldef solve(a, init):
M = [ ]
s = solve_eq(a)
for s2 in map(lambda k: s[:k] + s[k+1:], range(len(s))):
c = reduce(mul, map(lambda t: [ -t, 1 ], s2))
M.append(c)
b = [ inner_product(c, init) for c in M ]
Minv = mat_inv(M)
print_solution(s, Minv[0], b)solve([ 1, -9, 26, -24 ], [ 0, 1, 2 ])