目次
分数型漸化式 分数を出力する方法
前々回は漸化式をdef文で数列を計算する方法、前回は漸化式をsympy.rsolveを用いて一般項を求める方法を紹介しました。
sympy.rsolveで解くことができる漸化式は線形に限定されています。次のような分数型の漸化式をsympy.rsolveで解くことはできません。そこで今回は分数型の漸化式の一般項解法プログラムをつくってみます。
上のコードと実行結果を見て気づくことは、数列が小数で出力されていることです。初項が整数ですから、つづく第2項以降は有理数のはずです。
その原因は3行目のコードです。漸化式右辺の有理式〇/〇が有理数としてではなく割り算として計算されるためです。
数列を小数として欲しいのならこのままでもいいでしょうが、正体は有理数ですから有理数として出力したいところです。有理数として計算・出力する方法は2つあります。
分数出力方法 その1 fractions.Fraction
標準ライブラリであるfractionsモジュールをインポートして関数fractions.Fraction ()を使います。
>>> import fractions
fractions.Fraction (分子, 分母)またはfractions.Fraction (分子/分母)として使います。
最初のコードの分数式の部分
>>> else:return (a(n-1)-8)/(a(n-1)-5)
を
>>> else:return fractions.Fraction (a(n-1)-8, a(n-1)-5)
または
>>> else:return fractions.Fraction ((a(n-1)-8)/(a(n-1)-5))
に置き換えればOK！
次がその出力結果です。
fractions.Fraction ()の詳細はPython documentationを参照。
分数出力方法 その2 sympy.Rational
SymPyライブラリにある関数sympy.Rational()を使う方法です。
>>> import sympy
sympy.Rational (分子, 分母)またはFraction(分子/分母)として使います。
最初のコードの分数式の部分
>>> else:return (a(n-1)-8)/(a(n-1)-5)
を
>>> else:return sympy.Rational (a(n-1)-8, a(n-1)-5)
または
>>> else:return sympy.Rational ((a(n-1)-8)/(a(n-1)-5))
に置き換えればOK！
次がその出力結果です。
sympy.Rational()の詳細はSymPy1.11 documentationを参照。
分数型漸化式の解法アルゴリズム
分数型漸化式の一般項は次の手順で求めます。
- STEP1. 漸化式の項をtと置いて特性方程式（tの2次方程式）をつくる
- STEP2. 特性解を求める t=t1,t2
- STEP3. 特性解が異なる2解（ケースa）か重解（ケースb）であるかで場合分け
- STEP4a. 特性解が異なる2解の場合 b(n)=(a(n)-t1)/(a(n)-t2)とおく
- STEP4b. 特性解が重解の場合 b(n)=1/(a(n)-t1)とおく
- STEP5a. b(n)が等比数列となる（漸化式b(n+1)=br×b(n)）。公比brを計算
- STEP5b. b(n)が等差数列となる（漸化式b(n+1)=b(n)+d）。公差dを計算
- STEP6a. 一般項b(n)をsympy.rsolveに初項b(1)と漸化式を渡して算出
- STEP6b. 一般項b(n)を初項b(1)と公差dを用いて算出
- STEP7a. STEP4の関係式にb(n)を代入した式をa(n)についての方程式としてsympy.solveにより解いてa(n)を算出
- STEP7b. STEP4の関係式からa(n)=1/b(n)+t1を算出
分数型漸化式の解法プログラム
これまでの経緯をすべて投入して分数型漸化式の一般項を算出するプログラムをつくることできます。
>>> # 分数型漸化式の解法プログラム
>>> # 分数型漸化式の解法プログラム
>>> import sympy
>>> sympy.init_printing()
>>> from IPython.display import Math
>>>
>>> # 漸化式定義
>>> print('''
>>> 例題1【特性解が異なる2整数解の場合】
>>> a1 = 3
>>> p q r s = 1 -8 1 -5
>>>
>>> 例題2【特性解が重解（整数）の場合】
>>> a1 = 1
>>> p q r s = 1 -1 1 3
>>>
>>> 例題3【特性解が重解（有理数）の場合】
>>> a1 = 1
>>> p q r s = 8 -1 4 4
>>>
>>> 例題4【特性解が2虚数解の場合】
>>> a1 = 1
>>> p q r s = 1 -2 1 1
>>>
>>> 例題5【初期値に有理数1/2を渡す場合】
>>> a1 = 1/2
>>> p r r s = 2 1 3 -1
>>> ''')
>>>
>>> display(Math(r'a_1'))
>>> a1 = input('分数1/2もOK >>>')
>>> display(Math(r'a_{n+1}=\frac{pa_n+q}{ra_n+s}'))
>>> p,q,r,s=(int(x) for x in input('4つの整数を半角スペース切りで入力 p q r s >>>').split())
>>>
>>> #漸化式表示
>>> a1=str(a1)
>>> display(Math(r'a_1 = %s ,~~~a_{n+1} = \frac{%d a_n + %d}{%d a_n + %d}' %(a1,p,q,r,s)))
>>> a1=sympy.Rational(a1)
>>>
>>> # def文で数列を計算
>>> def A(n):
>>> if n == 1:return a1
>>> else:return sympy.Rational(p*A(n-1)+q, r*A(n-1)+s)
>>>
>>> # a(1)からa(10)まで出力
>>> print('def文による計算')
>>> for n in range(1, 11):
>>> print("{:>7}".format(f'a({n})= ')+f'{A(n)}')
>>>
>>> # 一般項計算
>>>
>>> # 特性方程式を解く
>>> sympy.init_printing()
>>> t = sympy.Symbol('t', complex=True)
>>> toku = sympy.Eq(r*t**2+(s-p)*t-q,0)
>>> tokukai = sympy.solve(toku, t)
>>> # 特性解表示
>>>
>>> t1 = sympy.Symbol('t1', complex=True)
>>> t2 = sympy.Symbol('t2', complex=True)
>>> kosu = len(tokukai)
>>> if kosu == 2:
>>> print('\n特性解')
>>> display(sympy.Eq(t1,tokukai[0]))
>>> display(sympy.Eq(t2,tokukai[1]))
>>> else:
>>> print('\n特性解（重解）')
>>> display(sympy.Eq(t1,tokukai[0]))
>>>
>>> # 異なる2解の場合の解法
>>> def kotonaru():
>>> n = sympy.Symbol('n', integer = True)
>>> an = sympy.Symbol('an', complex=True)
>>> a = sympy.Function('a')
>>> b = sympy.Function('b')
>>> # b(n)=bnを解く
>>> br = (p-tokukai[0]*r)/(p-tokukai[1]*r)
>>> fb = b(n+1)-br*b(n)
>>> b1 = (a1-tokukai[0])/(a1-tokukai[1])
>>> bn = sympy.rsolve(fb, b(n), {b(1):b1})
>>> # b(n)=bnとa(n)=anの関係からa(n)=anについて解いたものをkaiとする
>>> eq = sympy.Eq(bn*(an-tokukai[1]), an-tokukai[0])
>>> kai = sympy.solve(eq, an)[0]
>>> # 一般項出力
>>> print('一般項')
>>> display(sympy.Eq(a(n),sympy.simplify(kai)))
>>> # 数列の値の検証
>>> print('一般項による計算')
>>> for i in range(1, 11):
>>> #print(f'a({i})={sympy.simplify(kai.subs(n, i))}')
>>> print("{:>7}".format(f'a({i})= ')+f'{sympy.simplify(kai.subs(n, i))}')
>>>
>>> # 重解の場合の解法
>>> def jukai():
>>> n = sympy.Symbol('n', integer = True)
>>> a = sympy.Function('a')
>>> b = sympy.Function('b')
>>>
>>> b1 = sympy.Rational(1, a1-tokukai[0])
>>> a2 = sympy.Rational(p*a1+q, r*a1+s)
>>> b2 = sympy.Rational(1, a2-tokukai[0])
>>>
>>> d = b2-b1
>>> b = b1+(n-1)*d
>>> kai = 1/b+tokukai[0]
>>> # 一般項出力
>>> print('一般項')
>>> display(sympy.Eq(a(n),sympy.simplify(kai)))
>>> # 数列の値の検証
>>> print('一般項による計算')
>>> for i in range(1, 11):
>>> #print(f'a({i})={sympy.simplify(kai.subs(n, i))}')
>>> print("{:>7}".format(f'a({i})= ')+f'{sympy.simplify(kai.subs(n, i))}')
>>>
>>> # 解（解の個数）で解法を場合分け
>>> #kosu = len(tokukai)
>>> if kosu == 2:
>>> kotonaru()
>>> else:
>>> jukai()
使い方
初項a1と漸化式の4つの係数p、q、r、s（整数）を代入します。上部に5つの例題が表示されるのでその数値を入力欄にコピー&ペーストできます。4つの係数はスペース切りで「1 -8 1 -5」のように入力します。
出力
def文による最初の10項の計算結果
特性解
一般項
一般項による最初の10項の計算結果
実行結果
この「分数型漸化式の解法プログラム」は特性解が整数、有理数、複素数でも計算できるようにしてあります。そこで以下の5つの例題を解いてみます。
例題1 特性解が異なる2解の場合
例題2 特性解が重解（整数）の場合
例題3 特性解が重解（有理数）の場合
例題4 特性解が2虚数解の場合
例題5 初期値に有理数を渡す場合
例題1から例題3までの問題は、過去に手計算で解いた経験があります。分数型漸化式の解法アルゴリズムは特性解が複素数でも解けるはずだとは分かっていましたが、実際に解いたことはありませんでした。
今回このプログラムをつくってみて、特性解が複素数になる場合の一般項をはじめて見ることができました。一般項に見える虚数は消えて、結果として整数になる様子が分かります。ぜひみなさんもこのプログラムを使って分数型漸化式を解いてみてください。
今回使用したプログラムファイル（Jupyter Notebook）は次からダウンロードできます。
https://drive.google.com/file/d/15fB5t5eBO7bnr5ODGiPbByRMelKd_Wcy/view?usp=share_link
