9 漸化式の計算

9 漸化式の計算v1.nb
1
9 漸化式の計算
9.1 ループを用いる方法
関数 For[ini, test, incr, expr] は、最初に ini を実行し、test の値が True である限りtest 
expr  incr の順で処理を繰り返す。ある変数 j の値を 1 増やす（減らす）演算子 ++ (--) などと共に使われる
ことが多い。例として関数 For を用いて，パラメーター v をもつ漸化式
yt
mt
st
y0
s0





150  0.75 yt1  mt ,
v yt1  st1 ,
yt1 ,
200,
200,
の第 n 項 yn , mn , sn  を求める関数 yms[v,n] を作成してみよう。ちなみに，この漸化式は「乗数効果」と「加
速度原理」をもとにして景気循環を説明するもっとも基本的な経済モデルである。（yt と mt はそれぞれ t 期の
所得と投資）
In[1]:=
ymsv_, n_ : Modulej, y, m, s,
y  200;
s  200;
Forj  1, j  n, j ,
m  v y  s;
s  y;
y  150  0.75 y  m;
y, m, s
上で用いた関数 Module[{x,y,...},f[x,y,...]] は， x, y, ... などの記号を局所化し，他で使われて
いるx, y, ... などの記号をに影響を与えないようにして，f[x, y,...] の計算をするための関数であ
る。v=0.8 のときの y1 , m1 , s1  ,y5 , m5 , s5  ,y30 , m30 , s30  などを求めてみる。
In[2]:=
Out[2]=
In[3]:=
Out[3]=
In[4]:=
Out[4]=
yms0.8, 1
300., 0., 200
yms0.8, 5
804.238, 99.51, 739.638
yms0.8, 30
614.09, 3.68248, 613.877
v  0.8 のケースで yn , mn  の様子をグラフで確かめてみよう。まず，n  1 から 50 までの
yn , mn , sn  のデータを作成する。出力が長いので，関数 Short で表示を短くしている。
In[5]:=
Shortdata0  Tableyms0.8, n, n, 50
300., 0., 200, 48, 599.257, 0.349036, 598.545
Out[5]//Short=
とくに y と m に注目して，この点列の様子をプロットしてみよう。まず，data0 から yn , mn 
のリストだけを取り出す。 data0 を行列としてみると 3 列目にあたるのが s のデータであるから，1,
2 列目だけを取り出せばよい。一般に，ある行列 A の行 i1 から i2 ，列 j1 から
2
9 漸化式の計算v1.nb
j2 で構成される部分行列を作成するには，関数 TakeA, i1 , i2 , j1 , j2 
を用いる。今回のケースでは次のようになる。（行や列をすべて選ぶときは，単に All と書けばよい。）
In[6]:=
Shortdata1  Takedata0, All, 1, 2
300., 0., 455., 80., 46, 598.545,  0.0366442, 599.257, 0.349036
Out[6]//Short=
関数 ListLinePlota1 , b1 , a2 , b2 , ... , an , bn  は， xy 座標上の点列 a1 , b1 ,
a2 , b2 , ... , an , bn  を順に結んだ折れ線をプロットする。オプション
Mesh  All を付けることで，点列が強調表示される。
In[7]:=
ListLinePlotdata1, Mesh  All
60
40
20
Out[7]=
550
600
650
700
-20
-40
-60
さらに y だけの動きを確認してみよう。data0 から 1 列目だけを取り出せばよい。先程紹介した関数
Take を用いても良いが，ここでは，別の方法を試してみよう。一般に，行列
A の i 行 j 列を取り出すには， Ai, j と書けばよい。
In[8]:=
Out[8]=
data03, 2
124.
さらに，行列の j 列だけを取り出す場合，A[[All, j]] とすればよい。（以上の方法はいろいろと応用が利く
ので，ヘルプで「行列の部分抽出と設定」を検索し，確認しておくこと。）
In[9]:=
Out[9]=
data2  data0All, 1
300., 455., 615.25, 739.638, 804.238, 804.859, 754.141, 675.031, 592.986,
529.103, 495.721, 495.085, 520.805, 561.18, 603.185, 635.992, 653.241, 653.729,
640.687, 620.082, 598.578, 581.73, 572.819, 572.485, 579.097, 589.613,
600.622, 609.273, 613.877, 614.09, 610.738, 605.372, 599.736, 595.294,
592.916, 592.785, 594.484, 597.222, 600.107, 602.388, 603.616, 603.695,
602.834, 601.436, 599.96, 598.788, 598.154, 598.108, 598.545, 599.257
関数 ListPlot の別用法として, 単に ListPlota1 , a2 , ... , an  とすると，点 a1 ,
a2 , ... , an を順に点にプロットするようになる。また, デフォルトでは,
極端に大きい点や小さい点が描画から除外される可能性があるが,
オプション PlotRange  Allを指定すれば，すべての点が描画される。
9 漸化式の計算v1.nb
In[10]:=
3
ListLinePlotdata2, Mesh  All, PlotRange  All
800
700
600
Out[10]=
500
400
10
20
30
40
50
ü 考察
ここで紹介した data0 を作る方法 Table[yms[0.8, n], {n, 50}] は極めて計算の効率が悪く実用的ではな
い。実際 n の最大値が大きくなると，下図のように直線を越えるペースで計算時間が長くなる。何故か？考えて
みて欲しい。
jikan  TableTimingTableyms0.8, n, n, t1, t, 500
ListPlotjikan, AxesLabel  "n の最大値", "計算時間秒"
計算時間秒
0.8
0.6
0.4
0.2
100
200
300
400
500
n の最大値
9.2 再帰的定義を用いる方法
フィボナッチ数列
xt  xt1  xt2 , x0  x1  1,
の第 n 項を求める方法を考えよう。9.1 で見たようなループを用いた方法でxn を求めることは，簡単である（練習
9.3）。ここでは上の数列を
x t  x t  1  x t  2, x 0  x 1  1,
として眺めてみる。そうすると，関数 x(t) はその定義に自分自身への参照を含んでいることが分かる。このよう
な定義を再帰的定義と呼ぶ。 Mathematica では再帰的定義を次のように簡単に扱うことが出来る。
4
9 漸化式の計算v1.nb
In[11]:=
In[14]:=
Out[14]=
Clearx;
xt_ : xt  1  xt  2;
x0  x1  1;
Tablext, t, 0, 20
1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144,
233, 377, 610, 987, 1597, 2584, 4181, 6765, 10 946
再帰的定義は直観的に扱いやすいという利点があるが，その反面，計算速度やメモリー効率の点で利用しづらい
場合もある。例えば，x200 の計算を上の定義のまま行うのは，教室のPCでは困難である
再帰的定義を用いて，n 段階の枝をもつ樹木図を描こう。最初に「根」の座標を x=(0,0), 枝の長さをlen=1, 枝の方
向（y軸方向を0とし，時計回りに測った角度）をang=0とする。また，枝は継ぎ足されるごとに長さが fac=0.7 倍
になるものとし，もとの枝に対して，角度が±turn (turn=0.5)開いた方向に2本伸びるものとする。具体的な手続き
を次の通りである。
1. xから角度 ang の方向に長さlen の線分を描く。
2. n=0 なら終了。
3. n≠n-1, len≠fac*len, ang≠ang±turnとし，1から実行する。
In[15]:=
fac  0.7;
turn  0.5;
Treen_, x_, len_, ang_ : Moduley,
y  x  len Sinang, Cosang;
Linex, y,
Ifn  0,
Treen  1, y, fac  len, ang  turn,
Treen  1, y, fac  len, ang  turn,

In[18]:=
tr0  Tree8, 0, 0, 1, 0;
In[19]:=
Graphicstr0, AspectRatio  Automatic
Out[19]=
9 漸化式の計算v1.nb
5
9.3 関数 Nest などを用いる方法
初項がz0 = cである漸化式
zn  zn1 
zn1 3  1
,
3 zn1 2
の第n項を求める関数z[n,c]を作成してみよう。ここでは上の数列を
f x  x 
x3  1
,
3 x2
zn  f f f  ... f c, c に関数 f を n 階適用,
として眺めてみる。つまり，数列の第 n 項は，値 c に関数 f を n 階ネスト（入れ子）して適用したものであると
考えることが出来るということである。 Mathematica には値 c に関数 f を n 階ネストして適用する関数
Nest[f, c, n] が用意されている。この関数 Nest を用いることにより、z[n, c] は，
In[20]:=
fx_ : Nx  x ^3  1  3 x^ 2
zn_, c_ : Nestf, c, n
と定義できる。（厳密解を計算すると非常に時間がかかるので，関数 N を用いて近似解の計算をしてい
る。）例として，z[5,2], z[4, -1+]、z[4, -1-] を求めてみる。
In[22]:=
z5, 2
Out[22]=
1.
In[23]:=
z4, 1  I
Out[23]=
 0.5  0.866025 
In[24]:=
z4, 1  I
Out[24]=
 0.5  0.866025 
ところで，今考えている数列 zn は，初項 z0 = c の違い（虚数を許す）によって 1,  2 
1

3
2
1
,2 

3
2
の何
れかに収束することが知られている。収束の様子を知るには，c, f(c), f(f(c)), ... を順に調べていくと便利であ
る。Mathematica の関数 FixedPointList [f, c] は，c, f(c), f(f(c)), ... を計算結果がほぼ変化しなくなるまで繰
り返し，その結果をリストで与える。例として、z0 = -1 + Âのケースについて調べてみよう。
In[25]:=
Out[25]=
FixedPointListf,  1  I
1  ,  0.666667  0.833333 ,  0.508692  0.8411 ,  0.49933  0.866269 ,
 0.5  0.866025 ,  0.5  0.866025 ,  0.5  0.866025 , 0.5  0.866025 
よって，この場合zn は， 2 
1

3
2
U  0.5  0.866025 に収束することが分かる。収束先のみを知りたい場
合は，関数 FixedPoint[f, c] を用いる。
In[26]:=
FixedPointf, 2
Out[26]=
1.
In[27]:=
FixedPointf,  1  I
Out[27]=
 0.5  0.866025 
In[28]:=
FixedPointf,  1  I
Out[28]=
 0.5  0.866025 
6
9 漸化式の計算v1.nb
初項 z0 を a  b とし，a, b   2, 2 の範囲で、数列 zn  が 1, 

1

2

3
1
2


3
,
2
のうち，どれに収束するか調べてみよう。３つの収束先は，それぞれ虚部が 0, 
2
3
,
2
と異なることに注目する。したがって、ab平面 2, 2   2, 2 のそれぞれの点 a, b を，初項
2
z0  a  b の数列の収束先の虚部の大きさに応じて色で塗り分けたグラフィックスが描ければ便利である。
Mathematica の関数 DensityPlotfx, y, x, s, t, y, u, v は，関数 fx, y
の値に応じて黒の濃淡で xy 平面 x  s, t, yu, v を塗り分ける.これを利用してみよう。
3
In[29]:=
DensityPlotImFixedPointf, a  b I, a,  2, 2, b, 2, 2
Out[29]=
（注） Mathematica のバージョンの違いにより，グラフィックスがかなり異なる可能性がある。
上のままでは，様子が分かりづらい。オプション PlotPointØ200 を追加することで，縦横それぞれを200分割した
図を描くことが出来る（ただし調べる点が増える分，結果が表示されるまでにかなりの時間がかかるようにな
る。）。
9 漸化式の計算v1.nb
In[30]:=
DensityPlotImFixedPointf, a  b I, a,  2, 2, b, 2, 2, PlotPoints  200
Out[30]=
練習
ü 練習 9.1
上の関数 yms の定義において，
m  v y  s;
s  y;
y  150  0.75 y  m
の部分を入れ替えて，例えば，
s  y;
m  v y  s;
y  150  0.75 y  m
などとしてはならない。何故か。理由を述べなさい。
ü 練習 9.2　v  0.2, 1.2, 2.3 のそれぞれケ  スについて yt の様子をグラフで確かめなさい.
ü 練習 9.3
フィボナッチ数列 xt の定義を,
In[31]:=
Clearx;
xt_ : xt  xt  1  xt  2;
x0  x1  1;
7
8
9 漸化式の計算v1.nb
とすると，計算が高速化される。（x200 を計算してみよ。）何故か。理由を考えなさい。
ü 練習 9.4
微分方程式
y '  f x, y,
y a0   b0 ,
の解を y   x とする。このとき x, x  x0  n h ｈは任意の実数 の値は，漸化式
yt  yt1  h f xt1 , yt1 ,
xt  xt1  h,
x0  a0 , y0  b0 ,
で定まる数列 {yt }の第 n 項 yn で近似されることが知られている。このような近似解法を前進オイラー法とい
う。微分方程式
y '  3 x2 y,
y 0  1,
について次のことを調べてみよ。
(3.1) h =0.1 として，前進オイラー法を用いたy(x), x=0, 0.1, 0.2, ...,1の値を赤色の折れ線図で図示せよ。（オプショ
ンPlotStyleRedを用いる。）
(3.2) h =0.01 として，前進オイラー法を用いたy(x), x=0, 0.01, 0.02, ...,1の値を青色の折れ線で図示せよ。（オプシ
ョンPlotStyleBlueを用いる。）
3
(3.3) この微分方程式の厳密な解は y(x) =‰x である。厳密な解と上で求めた２種類の近似解の様子を図に描いて比
較せよ。この結果，h の大きさと近似の精度はどのように関係するといえるか。
ü 練習 9.5
初項がz0 = cである漸化式
zn  zn1 
zn1 5 1
,
5 zn1 4
の収束先と初項の関係を 9.3 同様に図を描いて調べよ。（色の違いをより鮮明にするには，DensityPlot にオ
プション ColorFunction(Hue[0.8#+0.1]&) をつけて実行すると良い。）
参考文献
[1] 奥村晴彦，1991，C言語による最新アルゴリズム辞典，技術評論社．

Download Report