漸化式

例えば、次のような斉次線形漸化式

a_n+2 - 5a_n+1 + 6a_n = 0

には、次の方程式が対応しています。

x² - 5x + 6 = 0

なぜなら、上の漸化式を次の形にしようとすると、

a_n+2 - αa_n+1 = β(a_n+1 - αa_n)

これを展開して、

a_n+2 - (α + β)a_n+1 + αβa_n = 0

だから、

α + β = 5
αβ = 6

となり、上の方程式の解になっています。これは4項以上の漸化式でも同じです。
ここで、次の漸化式を解いてみましょう。

a_n+3 - 9a_n+2 + 26a_n+1 - 24a_n = 0

対応する方程式

x³ - 9x² + 26x - 24 = 0

の解は、x = 2, 3, 4なので、

a_n+3 - 7a_n+2 + 12a_n+1 = 2(a_n+2 - 7a_n+1 + 12a_n)
a_n+3 - 6a_n+2 + 8a_n+1 = 3(a_n+2 - 6a_n+1 + 8a_n)
a_n+3 - 5a_n+2 + 6a_n+1 = 4(a_n+2 - 5a_n+1 + 6a_n)

とりあえずこれを再帰的に解いて、

a_n+2 - 7a_n+1 + 12a_n = -5 • 2ⁿ
a_n+2 - 6a_n+1 + 8a_n = -4 • 3ⁿ
a_n+2 - 5a_n+1 + 6a_n = -3 • 4ⁿ

あとは連立方程式を解くだけです。

a_n = -3/2 • 4ⁿ + 4 • 3ⁿ -5/2 • 2ⁿ

重解があるときは少しややこしくなりますが、そうでないときはこのように簡単に解けます。
上の手順をコードにしてみました。

from fractions import Fraction
def 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 m
def 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 h
def 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_sol
def 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 ])