コラム 人と星とともにある数学 数学

SymPyで漸化式を解く 番外編|Pythonで数学を学ぼう! 第32回

分数型漸化式 分数を出力する方法

前々回は漸化式を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

  • この記事を書いた人
  • 最新記事

桜井進(さくらいすすむ)様

1968年山形県生まれ。 サイエンスナビゲーター®。株式会社sakurAi Science Factory 代表取締役CEO。 (略歴) 東京工業大学理学部数学科卒、同大学大学院院社会理工学研究科博士課程中退。 東京理科大学大学院非常勤講師。 理数教育研究所Rimse「算数・数学の自由研究」中央審査委員。 高校数学教科書「数学活用」(啓林館)著者。 公益財団法人 中央教育研究所 理事。 国土地理院研究評価委員会委員。 2000年にサイエンスナビゲーターを名乗り、数学の驚きと感動を伝える講演活動をスタート。東京工業大学世界文明センターフェローを経て現在に至る。 子どもから大人までを対象とした講演会は年間70回以上。 全国で反響を呼び、テレビ・新聞・雑誌など様々なメディアに出演。 著書に『感動する!数学』『わくわく数の世界の大冒険』『面白くて眠れなくなる数学』など50冊以上。 サイエンスナビゲーターは株式会社sakurAi Science Factoryの登録商標です。

あわせて読みたい

-コラム, 人と星とともにある数学, 数学
-,