分数型漸化式の問題をPythonで解いてみた。

　知恵袋の分数型漸化式の一般項を求める問題をPythonで解いてみました。(^_^;
　特性方程式が重解をもつ場合(1)と異なる2解をもつ場合(2)について、公式から一般項を求めて、その式から求めた値と漸化式から漸次求めた値を比較して検算しました。

　 $(1)\quad a_1=-17,\quad a_{n+1}=\frac{17a_n+9}{-a_n+23}$
　 $(2)\quad a_1=8,\qquad a_{n+1}=\frac{3a_n+2}{a_n+2}$

　分数型漸化式の一般項の式は、次式を用いてsympyを使って求めました。(^_^;

　分数型漸化式
　 $a_{n+1}=\frac{p a_n+q}{r a_n+s}$
(ただし、r≠0)の特性方程式
　 $x=\frac{px+q}{rx+s} \qquad\therefore rx^2+(s-p)x-q=0$
の解をα、βとすると、漸化式の一般項は、次式で表される。
(Ⅰ) α≠βのとき、
　 $a_n=\frac{(a_1-\alpha)\beta(r\beta+s)^{n-1}-(a_1-\beta)\alpha(r\alpha+s)^{n-1}}{(a_1-\alpha)(r\beta+s)^{n-1}-(a_1-\beta)(r\alpha+s)^{n-1}}$
　 $\quad =\alpha+\frac{(\alpha-\beta)\left(\frac{a_1-\alpha}{a_1-\beta}\right)}{\left(\frac{r\alpha+s}{r\beta+s}\right)^{n-1}-\left(\frac{a_1-\alpha}{a_1-\beta}\right)}$
(Ⅱ) α＝βのとき、
　 $a_n=\frac{r\alpha(a_1-\alpha)(n-1)+a_1(r\alpha+s)}{r(a_1-\alpha)(n-1)+(r\alpha+s)}=\frac{r\alpha(a_1-\alpha)n+r\alpha^2+sa_1}{r(a_1-\alpha)n+r(2\alpha-a_1)+s}$
　 $\quad =\alpha+\frac{(a_1-\alpha)(r\alpha+s)}{r(a_1-\alpha)(n-1)+(r\alpha+s)}$

ただし、初項を $a_1=a$ としてもいいのですが、aとα(ｱﾙﾌｧ)を見分けやすいように数字を付けておきました。(^_^;
　ちなみに、自力で解くとき公式が憶えられない場合には、A,Bを定数として、一般項が次式の形になるのがわかっているので、n=1,2の場合についてAとBの連立1次方程式を立ててA,Bを決定すればよい。

(Ⅰ) α≠βのとき、
　 $a_n=\alpha+\frac{1}{A \left(\frac{r\alpha+s}{r\beta+s}\right)^{n-1}+B} \qquad \rightarrow \qquad A \left(\frac{r\alpha+s}{r\beta+s}\right)^{n-1}+B=\frac{1}{a_n-\alpha}$
(Ⅱ) α＝βのとき、
　 $a_n=\alpha+\frac{1}{A(n-1)+B} \qquad \qquad \rightarrow \qquad A(n-1)+B=\frac{1}{a_n-\alpha}$

　数学の解答には使えませんが、検算や結果だけが知りたい場合などには便利かと思います。(^_^;

● fracRecurrenceRel1.py

# coding: UTF-8
# fracRecurrenceRel1.py

from time import time
from fractions import Fraction
from sympy import *

# 漸化式を使って数列の初項から第n項までを列挙
# lc = [a1,(p,q,r,s)]
def la(n,lc):
    [a1,(p,q,r,s)] = lc
    a = Fraction(a1,1)
    l = [a]
    for i in range(2,n+1):
        a = Fraction(p*a+q,r*a+s)
        l.append(a)
    return l

# 約数を列挙
def divisor(p):
    return [n for n in range(1,p+1) if p%n==0]

# 因数定理で多項式の有理数の解を求める
def factorThm(la):
    result = []
    n = len(la)
    if n==0 or la[0]==0: return []  # Error処理
    if la[-1]==0:
        result.append(0)
    m = n-1
    while m> 0 and la[m]==0:
        m-=1
    nums = divisor(abs(la[m]))
    dens = divisor(abs(la[0]))

    for sgn in (-1,1):
        for num in nums:
            for den in dens:
                x = sgn*Fraction(num,den)
                if x in result: continue
                if x.denominator==1:
                    x = x.numerator
                f = 0    # Horner法
                for i in range(n):
                    f = la[i]+f*x
                if f==0:
                    result.append(x)
    return result

# 公式を使って数列の第n項の値を求める
def a(n,lc):
    [a1,(p,q,r,s)] = lc
    li = factorThm([r,s-p,-q])
    if len(li)==0: return float('nan')  # Error処理
    if len(li)==1: li=li*2
    Alp,Bet = sorted(li)
    if Alp==Bet:
        num = r*Alp*(a1-Alp)*n+r*Alp**2+s*a1
        den = r*(a1-Alp)*n+r*(2*Alp-a1)+s
        return Fraction(num,den)
    num = (a1-Alp)*Bet*(r*Bet+s)**(n-1)-(a1-Bet)*Alp*(r*Alp+s)**(n-1)
    den = (a1-Alp)*(r*Bet+s)**(n-1)-(a1-Bet)*(r*Alp+s)**(n-1)
    return Fraction(num,den)

# 公式を使って数列の第n項の式を求める
def fa(lc):
    [a1,(p,q,r,s)] = lc
    li = factorThm([r,s-p,-q])
    if len(li)==0: return float('nan')  # Error処理
    if len(li)==1: li=li*2
    Alp,Bet = sorted(li)
    n = Symbol('n')
    if Alp==Bet:
        num = r*Alp*(a1-Alp)*n+r*Alp**2+s*a1
        den = r*(a1-Alp)*n+r*(2*Alp-a1)+s
        return simplify(num/den)
    num = (a1-Alp)*Bet*(r*Bet+s)**(n-1)-(a1-Bet)*Alp*(r*Alp+s)**(n-1)
    den = (a1-Alp)*(r*Bet+s)**(n-1)-(a1-Bet)*(r*Alp+s)**(n-1)
    return simplify(num/den)

def test(n,lc):
    [a1,(p,q,r,s)] = lc
    print('x = %s'%factorThm([r,s-p,-q]))   # 特性方程式の有理数の解
    print('a_n = %s'%fa(lc))                # 一般項の式
    lA1 = la(n,lc)                          # 漸化式から求めた数列
    lA2 = [a(i,lc) for i in range(1,n+1)]   # 一般項から求めた数列
##    print(lA1[n-1])
##    print(lA2[n-1])
    print(lA1==lA2)                         # 2通りの方法で求めた数列を比較

def main():
    tm = time()  # Timer Start
    N = 30
    print('(1)')
    test(N,[-17,(17,9,-1,23)])
    print('(2)')
    test(N,[8,(3,2,1,2)])
    print("Runtime : %.3f [sec]"%(time()-tm))   # Timer Stop & Disp

if __name__ == '__main__':
    main()

●実行結果

(1)
x = [3]
a_n = 3 - 20/n
True
(2)
x = [-1, 2]
a_n = (6*4**n + 8)/(3*4**n - 8)
True
Runtime : 0.250 [sec]

※参考URL
●［PDF］分数型の漸化式
●完全数もPythonでやってみた - 牌語備忘録 -pygo
●3.2. Sympy : Python での代数計算 − Scipy lecture notes - turbare.net