Chapter 14
力学系（補講１）「力学系」とは「時間発展する世界」のごく単純なモデルである．本章では以下の 3 項目について解説する．

• グラフ解析と呼ばれる方法で実 1 次元力学系を可視化．
• 方程式の近似解法であるニュートン法を可視化．
• フラクタル図形として有名な複素 1 次元力学系のジュリア集合・マンデルブロー集合の描画．

注意． 本章は「補講」，すなわちオマケです．13 章までの内容は既知とさせていただ
きます．説明も親切さに欠ける部分があるかもしれませんがご了承ください．

14.1

力学系とは？

まずは力学系 (dynamical system) とは何か，具体例をみながら説明してみよう．
次のような漸化式を考える：

an+1 = 2an ; a0 = 1.
この漸化式を解くのは簡単だが（an = 2n ），ちょっと遠回りをして次のように解釈してみる．いま f (x) = 2x とおくと，

an+1 = f (an )と表されるから，a1 = f (a0 ), a2 = f (f (a0 ))，a3 = f (f (f (a0 )))，という具合に関数f を繰り返し合成（反復合成）することで数列が生成されていることがわかる．初項
a0 はたまたま 1 としているが，この値を変化させても数列を生成する仕組み自体は変化しないのである．

214

力学系（補講１）

ところで，ニュートン力学のような素朴な世界観では，世界（宇宙）を構成する要素として「空間」と「時間」，そして物体の運動を記述する「運動法則」（物理法則）*1 空間内に質量をもった複数の物体が配置されて，時間とともに全体のの 3 つがある．

「系」が変化するのである．
同様に，上の数列 {an } が棲む「世界」は次のように記述される：• 空間：数直線 R• 時間：0, ±1, ±2, · · · （秒）• 運動法則：x の位置にある動点は 1 秒後に f (x) = 2x の位置に移動する
この「世界」を「R 上の f による力学系 (dynamical system)」と呼ぶ．*2 また，点x0 ∈ R にたいしfffx0 7−→ f (x0 ) 7−→ f (f (x0 )) 7−→ · · ·fff（すなわち x0 7−→ 2x0 7−→ 22 x0 7−→ · · · ）で定まる数列を初期値 x0 の軌道 (orbit)と呼ぶ．たとえば数列 {an } は初期値 a0 = 1 の軌道だと考えられる．このようにして得られる軌道の振る舞いを調べるのが，力学系理論の目標である．
また，f ◦ f ◦ f ◦ · · · ◦ f と f を n 回合成したものを f n と書くことにすれば，上の*3 ，力学系理論では本質的に「関数 f の反復合成に軌道は {f n (x0 )} と表されるから，よって得られる数列の振る舞いを記述すること」が目標となる．しかし，一般にはこれが至って難しい．たとえば f (x) が 2 次関数の場合，f 10 (x) は x の 1024 次関数である．手計算で力学系の時間発展を追いかけ続けるのは至難の業，実験もままならない—— というのはあくまで昔の話．いまは便利なパソコンと，Mathematica があるではないか．

14.2

力学系のグラフ解析

与えられた実数 p と関数 y = f (x) に対し，軌道

p 7→ f (p) 7→ f (f (p)) 7→ f (f (f (p))) 7→ . . .
を視覚化する非常によい方法がある．y = f (x) のグラフを用いて，xy 平面上の点(p, p) から (f (p), f (p)) を作図する方法である：*1 相対性理論によれば，これらの 3 要素を明確に区別することはできない．
*2 より正確には離散力学系 (discrete dynamical system) と呼ばれるものである．時間が「ぱらぱら漫画」のように離散的だからである．歴史的にはポアンカレが微分方程式系で記述される（連続時間をもつ）力学系の理論を創始した．そこでは離散力学系が重要な「道具」として用いられる．
*3 f n (x) は f ◦ f ◦ f ◦ · · · ◦ f (x) のことで f (x) の n 乗 {f (x)}n ではない．

14.2

力学系のグラフ解析

3

Step 0 y = f (x) と y = x のグラフを描く．
Step 1 点 (p, p) からタテ方向に，y = f (x) のグラフまで線分を引く．その交点は(p, f (p)) である．
Step 2 (p, f (p)) からヨコ方向に，y = x のグラフとまで線分を引く．その交点が(f (p), f (p)) である．
この Step 1 と 2 を続けていき，直線 y = x を数直線だと思えば点 p の軌道が視覚化される（図 14.1）．この手続きをグラフ解析 (graphical analysis) と呼ぶ．また，グラフ解析で得られる図をウェブ・ダイアグラム（web diagram）と呼ぶ．

yy=xf (p)y = f (x)pOpf (p)x

図 14.1 f (x) のグラフ解析．(p, p) から (f 2 (p).f 2 (p)) までを作図したところ．

たとえば f (x) = 2x で初期値 x0 = ±1/4 とした場合は図 14.2 のようになる．
これを Mathematica に作図させる方法を紹介しよう．基本的にはウェブ・ダイア
グラムの折れ線を ListLinePlot で描画する．これに Show で関数のグラフを重ね合わせるのである．
[1] Step 0 のように，関数 y = f (x) = 2x と y = x のグラフを描く：
In[1]:=

f[x_] := 2 x;
gr = Plot[{f[x], x}, {x, -2, 8}, AspectRatio -> Automatic]

4力学系（補講１）1426-5-4-3-21-124-22-4-21-12345-6-2

図 14.2 f (x) = 2x のグラフ解析．左は x0 = 1/4 のとき，右は x0 = −1/4 のとき．

1510Out[1]=5-22468[2] Step 1 と Step 2 に対応する部分．点 (p, p) からタテに移動して (p, f (p)) へ，さらにヨコに移動して (f (p), f (p)) へ，という「１操作分」の折れ線を追加するための関数を定義：
In[2]:=

tateyoko[p_] := {{p, f[p]}, {f[p], f[p]}};[3] 初期値 p に対して，上記の「１操作分」を n 回繰り返したリストを作る関数を定義：In[3]:=

weblist[p_, n_] := (w = {{p, p}}; x = p;
Do[(w = Join[w, tateyoko[x]]; x = f[x]), {i, 1, n}];w);[4] たとえば…：
In[4]:=

Out[4]=

weblist[1/2, 3]
1 1
1:: , >, : , 1>, 81, 1<, 81, 2<, 82, 2<, 82, 4<, 84, 4<>
2 22

14.2

力学系のグラフ解析

5[5] ウェブ・ダイアグラムを作る関数：
In[5]:=

webdiag[p_, n_] := ListLinePlot[weblist[p, n],PlotStyle -> Thick, AspectRatio -> Automatic,PlotRange -> All];

最後の PlotRange -> All はおまじないのようなもので，ないとうまく描画されないことがある．
[6] たとえば…：
In[6]:=

webdiag[1/2, 3]
43Out[6]=

21

123

4

[7] Show を用いて最初のグラフと重ねて表示：
In[7]:=

Show[gr, webdiag[1/2, 3]]
15

10

Out[7]=

5-224

6

8[8] Manipulate で p と n を変化させる：
In[8]:=

Manipulate[Show[gr, webdiag[p, n]],{{p, 1}, -1, 4}, {{n, 3}, 0, 10, 1}]614

力学系（補講１）

pn

15

Out[8]=

105-22468n は整数のみを動くように指定していることに注意．
問題 14.1 (ロジスティック写像) 2 次関数 ga (x) = ax(1 − x) は 0 ≤ a ≤ 4 のとき区
間 [0, 1] から 区間 [0, 1] への関数となる．そのウェブ・ダイアグラムを Manipulateを用いて可視化せよ．ただし，a もパラメーターとして変化できるようにすること．
【解答】 基本的に [1] – [8] と同じである．パラメーター a が新たに加わる部分だけ注意すればよい．
In[ ]:=

g[a_, x_] := a x (1 - x);
gr2[a_] := Plot[{g[a, x], x}, {x, 0, 1},
AspectRatio -> Automatic]
tateyoko2[a_, p_] := {{p, g[a, p]}, {g[a, p], g[a, p]}};weblist2[a_, p_, n_] := (w = {{p, p}}; x = p;
Do[(w = Join[w, tateyoko2[a, x]]; x = g[a, x]),{i, 1, n}]; w);
webgr2[a_, p_, n_] := ListLinePlot[weblist2[a, p, n],
PlotStyle -> Thick, PlotRange -> All]
Manipulate[Show[gr2[a], webgr2[a, p, n]],{{a, 3}, 0, 4}, {{p, 0.2}, 0, 1}, {{n, 5}, 0, 10, 1}]

以上のように入力すると，次のような結果が得られる：

14.2

力学系のグラフ解析

7apn1.00.8

0.6

0.4

0.20.20.4

0.60.81.0■問題 14.2 (より見やすく) 問題 1 のウェブ・ダイアグラムを次のように変えて「見やすさ」を追求せよ：(1) ダイアグラムの線分は黒に．
(2) ダイアグラムで，y = x 上の点は赤く，グラフ上の点は青く．
【解答】 ListLinePlot のオプションではうまくいかないので，最初から折れ線を Graphicsで作ってしまおう．Graphics の使い方の詳細はドキュメントセンターを参照されたい．
まずはタテヨコの黒い線分と赤の点，青の点を作成するための関数を定義：
In[ ]:=

tateyoko3[p_, q_] := (pp = {p, p}; pq = {p, q}; qq = {q, q};{Line[{pp, pq, qq}],{Red, PointSize[Medium], Point[{pp, qq}]},{Blue, PointSize[Medium], Point[pq]}});

基本的には q に g[a, p] を入れるのだが，汎用性と式の見やすさを考慮してこのようにした．
次に関数 tateyoko3 を n 回施す関数を定義：
In[ ]:=

weblist3[a_, p_, n_] := (w = {{Red, PointSize[Large], Point[{p, p}]}}; x = p;
Do[(w = Join[w, tateyoko3[x, g[a, x] ]]; x = g[a, x]),{i, 1, n}]; w);

814

力学系（補講１）

Manipulate で描画：
In[ ]:=

Manipulate[ Show[gr2[a], Graphics[weblist3[a, p, n]]],{{a, 3}, 0, 4}, {{p, 0.25}, 0, 1}, {{n, 5}, 1, 10, 1}]

以上のように入力すると，次のような結果が得られる：

apn1.00.80.60.40.20.20.40.60.81.0■14.3

ニュートン法与えられた関数 f (x) に対し方程式 f (x) = 0 の解を数値的に求める方法として，ニュートン法 (Newton’s method) というものがある．幾何学的（直感的）には，次のような手続きによって解の近似値を得るのである：(1) 関数 y = f (x) のだいたいのグラフを描き，解の場所に見当をつける．
(2) 解の近くにあると思われる値 x0 を選び，グラフ上の点 (x0 , f (x0 )) における接線を引く．

(3) 接線と x 軸の交わる点を (x1 , 0) とする．

14.3

ニュートン法9

もし x0 が方程式の解 f (x) = 0 の解 α に十分に近ければ，x1 は解 α の「さらに良い近似値」を与えることが証明できる．*4 さらに (x1 , f (x1 )) からグラフの接線を引いて同様の操作を繰り返せば，近似はもっと良くなるであろう．簡単な計算により

x1 = x0 − f (x0 )/f ′ (x0 ) であることがわかるから，
Nf (x) := x −f (x)f ′ (x)とおいて得られる（関数 Nf によって得られる力学系の）軌道
Nf

Nf

fx0 7−→ x1 = Nf (x0 ) 7−→ x2 = Nf2 (x0 ) 7−→ · · ·は解 α への収束列を与えると期待される．この関数 Nf を f のニュートン写像とい図 14.3 ニュートン法う．この節では，ニュートン法で解が近似される様子を数値的に観察し，それをグラフ上で可視化する方法を紹介しよう．
[9] f (x) = x2 − 2 とすれば，ニュートン法で得られる列は√2 の近似列を与えるはずである．その様子を見てみよう．まずは関数 f を上書きし，さらにニュートン写像を定義：
In[9]:=f[x_]= x^2 - 2;

df[x_] = D[f[x], x];
newton[x_] = x - f[x]/df[x]*4 ここで，関数 f は最低でも接線を引ける程度に滑らかでなくてはならない．たとえば C 2 程度を仮定すると十分である．

1014

Out[9]=x-力学系（補講１）- 2 + x2
2xx2 + 2すなわち Nf (x) =であることがわかる．
2x[10] NestList を用いて，たとえば x0 = 1 から始めたときの x0 , x1 , . . . , x5 をリストにする．：
In[10]:=

Out[10]=app = NestList[newton, 1, 5]
3 17 577 665 857 886 731 088 897
,,,>
:1, ,
2 12 408 470 832 627 013 566 048[11] これでは近似の様子が見えてこないので，実用上はあまり意味がない．有効数字

20 桁で数値化する：
In[11]:=

N[app, 20] //TableForm

Out[11]//TableForm=

1.0000000000000000000
1.5000000000000000000
1.4166666666666666667
1.4142156862745098039
1.4142135623746899106
1.4142135623730950488

TableForm を用いて縦方向に並べた．
[12] 真の値 Sqrt[2] との誤差を見てみよう*5 ：In[12]:=

N[{app, app - Sqrt[2]}, 20] //Transpose//TableForm

Out[12]//TableForm=

1.0000000000000000000
1.5000000000000000000
1.4166666666666666667- 0.41421356237309504880
0.085786437626904951198
0.0024531042935716178650

1.4142156862745098039

2.1239014147551198799 ´ 10-6

1.4142135623746899106

1.5948618246068546804 ´ 10-12

1.4142135623730950488

8.9929283216504531005 ´ 10-25

右側に真の値との誤差が並んでいる．誤差は「急速に」小さくなっていることが分かるだろう．大雑把に言って，１ステップごとに真の値と合致する小数点以下の桁数が倍近くになる．方程式の近似解法の中でも，ニュートン法はたいへん効率がよい．

*5 この式はちょっと技巧的かもしれない．これはTranspose[N[{app, app - Sqrt[2]}, 20]] //TableForm と同じ意味であり，一旦転置を取ってから表にすることで誤差が右側に並ぶようにしたのである．次のように書いても同様の表示が得られる：

TableForm[N[{app, app - Sqrt[2]}, 20], TableDirections -> Row]

14.3

ニュートン法11[13] 次に図示してみよう．基本的にはウェブ・ダイアグラムを描く要領でやればよい：
In[13]:=

seg[p_] := {{p, f[p]}, {newton[p], 0}};seglist[p_, n_] := (w = {{p, 0}}; x = p;
Do[ (w = Join[w, seg[x]]; x = newton[x]), {i, 1, n}]; w);
seggr[p_, n_] := ListLinePlot[seglist[p, n],PlotRange -> All, PlotStyle -> Thick]gr = Plot[f[x], {x, -5, 5}];
Manipulate[ Show[gr, seggr[p, n]],{{p, 5}, -5, 5}, {{n, 3}, 1, 10, 1}]

pn20

Out[13]=15105-4-224√√
p が負のときは − 2 に収束し，p が正のときは 2 に収束することがわかる．*6
解答は省略するが，問題 1 を参考にぜひ次の問題 3 も考えてみていただきたい．
問題 14.3 (解に収束しない例) 3 次関数 ga (x) = x3 − 3x + a (a > 0) に対して，ニュートン法の軌道を Manipulate を用いて可視化せよ．ただし，a もパラメーターとして変化できるようにすること．また a の値によっては下の図のような軌道が現れて，必ずしも方程式 ga (x) = 0 の解に収束するとは限らない（と思われる）ことを確
認せよ．

*6 ちなみに，p がちょうど 0 の場合は f の微分がゼロになるので数列が定義できなくなる．実際は軌道が「無限遠点」という固定点に乗ってしまう．

1214

力学系（補講１）apn10

52-24-514.4

複素力学系とジュリア集合

今度は「複素平面 C 上の，関数 f (z) による力学系」を考える．いわゆる（1 次元）複素力学系 (complex dynamics) と呼ばれるものである．ここでは関数 f (z) をfc (z) = z 2 + c (c ∈ C) の形の 2 次多項式に制限して，その力学系における軌道のふるまい理解したいとしよう．ここで fc の定数項 c ∈ C は複素数のパラメーターであり，力学系における物理定数のようなものだと考えたらよいかもしれない．このパラメーターを変化させることで，「世界」の様相は意外なほど変化するのである．
無限の鉢と充填ジュリア集合． fc の力学系に共通の性質として，『|z| がある程度大きいとき，|fcn (z)| → ∞ (n → ∞) が成り立つ』ということが証明できる．*7 したがって「軌道（の絶対値）が無限大に発散する初期値」の集合

Bc := {z ∈ C | |fcn (z)| → ∞ (n → ∞)}の形をみれば，パラメーター c に依存する力学系の変化を捉えることができるかもしれない．この集合を fc の無限の鉢 (basin at infinity) と呼ぶ．また，Bc の補集合

Kc := C − Bc
*7

たと えば，|z| ≥ 2 + |c| のとき |fc (z)| ≥ 2|z| が成り立つ．実際，|f (z)| = |z 2 + c| ≥|z| · |z| − |c| ≥ (2 + |c|)|z| − |c| = 2|z| + |c|(|z| − 1) ≥ 2|z| + |c|(1 + |c|) ≥ 2|z|. よって |f n (z)| ≥ 2n |z| → ∞ (n → ∞)．もう少し精密に評価すると，|z| ≥ max{2, |c|} のとき|f n (z)| → ∞ (n → ∞) であることが証明できる．

14.4

複素力学系とジュリア集合

13を充填ジュリア集合 (filled Julia set) と呼ぶ．この集合は「軌道が有限の範囲にとどまるような初期値の集合」であることも知られている．
ちなみに Bc と Kc の境界部分は「発散する軌道」と「有界にとどまる軌道」がせめぎあう部分であり，わずかな誤差で前者にも後者にも変化しうる．すなわち，この力学系における「カオス部分」である．これを fc のジュリア集合 (Julia set) と呼び，Jc で表す．面白いことに，一般に Jc は自己相似性をもつフラクタル集合となるのである．
では，Mathematica でこれらの集合を描く方法を考えてみよう．いろんな描き方があるが，ここでは計算効率は度外視して数学的な単純さ・分かりやすさを優先した次の方法を紹介する．*8
まず話を簡単にするために，|c| ≤ 2 と仮定しておこう．このとき，与えられたz ∈ C に対し |fck (z)| ≥ 2 となる k が存在すれば，|fck+n (z)| は n に関して単調増加かつ発散することが証明できる．したがって，{}Bc (k) := z ∈ C | |fck (z)| ≥ 2とすれば，Bc (1) ⊂ Bc (2) ⊂ · · · かつ Bc =∪k≥1 Bc (k) である．よって k が十分大きいとき Bc (k) は Bc を近似していると考えてよいだろう．以下では k = 50 としている．
[14] まずは c = −0.122 + 0.745i に対して Bc と Kc を描いてみよう．f の定義を上書きする．汎用性を考えて c の定義と別にしておく：
In[14]:=

c = -0.122 + 0.745 I; f[z_] := z^2 + c;

[15] さて無限の鉢 Bc に色をつける関数を定義：
In[15]:=

col[z_] := (p = z; k = 0;
While[(Abs[p] < 2.0) && (k < 50),(p = f[p]; k = k + 1)];k)

While の中の意味は「|fck (z)| < 2 かつ k < 50 である限り，カッコ (...)

*8

内のふた

実際のところ，描画速度を優先するならば Mathematica を使うべきではない．C や Java のようなプログラミング言語を用いるべきである．複素力学系の計算に Mathematica を用いる利点は，• 数値がデフォルトで複素数であり，計算が極めて容易• グラフィクスやユーザー・インターフェイスが組込み関数で準備されているという部分であって，C 言語などに比べて「プログラムを書きあげるまでの時間」が大幅に短縮できるのである．

1414

力学系（補講１）

つの命令 p = f[p]; k = k + 1 を繰り返せ，という意味である．最終的には k が値として定まるが，• |fck (z)| ≥ 2 となる k ≤ 49 が見つかった．
• k = 50 になるまで |fck (z)| < 2 であり続けた．のいずれかである．すなわち，col は 0 以上 50 以下の整数値をとる関数であって，前者の場合（49 以下）は z ∈ Bc とみなし，後者の場合（ちょうど 50）は z ∈ Kc とみなすことにする．
[16] 複素平面の領域{x + yi | −2 ≤ x ≤ 2, − 2 ≤ y ≤ 2}
上を刻み幅 d = 0.01 で分割し，Table を用いて関数 col の値をリスト（のリスト，すなわち行列と同じ形のデータ）にする：
In[16]:=

d = 0.01;
tab = Table[col[x + y I], {x, -2, 2, d}, {y, -2, 2, d}];ここの計算が一番時間がかかる．数十秒は我慢しよう．*9[17] tab の値を ArrayPlot で表現するために，少し修正した complexAP (complexArray Plot の略) を導入する（この意味は図 14.4 参照．）：
In[17]:=*9complexAP[t_] := ArrayPlot[Reverse[Transpose[t]]];

400 × 400 = 160000 個の複素数についてそれぞれ col を計算させることになる．パソコンに長時間向かって肩もこっていることだろうから，ストレッチでもしながらのんびり待つとよい．ちなみに
col 関数の定義内の 50 を減らすか，刻み幅 d を増やせば計算は速くなる．そのかわり，得られる画像は粗くなる．14.4

複素力学系とジュリア集合15

図 14.4 complexAP の意味．左上の複素平面を刻み幅 d でグリッドに分割してcol 関数を適用するが，それを tab のように行列型のデータにすると右上のような配置になる．これに Transpose（行列の転置），Reverse（行の前後入れ替え）を施すことで，もとの複素平面のグリッドとデータの配置が一致する．これをArrayPlot で表示するのである．

[18] tab を complexAP で図示する：
In[18]:=complexAP[tab]Out[18]=col の値が小さいほど白くなる．したがって白からグレーの部分が Bc であり，黒い部分が Kc である．このように意外と複雑なフラクタル集合となる．

1614

力学系（補講１）

問題 14.4 (より美しく？) 上の Bc を白黒でなく，適当に色づけせよ．
【解答】 ArrayPlot のオプションには ColorFunction というのがあるので，それで指定できる．（詳しくはドキュメントセンターを参照せよ．
）たとえば，complexAP を修正してIn[ ]:=complexAP2[t_] := ArrayPlot[Reverse[Transpose[t]],ColorFunction -> "LightTemperatureMap"]; としてみよう．これに [16] で作ったデータ tab を放り込んで complexAP2[tab] とすればよい．結果は図 14.5 の左上のようになる．

■

図 14.5 ColorFunction のオプションはドキュメントセンターの『カラースキーム』の項に詳しく述べられている．左上："LightTemperatureMap"，右上："MintColors", 左下："WatermelonColors"，右下："RedBlueTones".

c の値をいろいろと変えても面白い．たとえば図 14.6 のようになる．
マンデルブロー集合.

fc の力学系において，z = 0 は fc′ (z) = 0 となる唯一の点である．そのような特殊性から，z = 0 の軌道によって力学系の性質がある程度分類できることが知られている．たとえば「fc の力学系において z = 0 の軌道（の絶対値）が発散する（すなわち 0 ∈ Bc ）」という性質を持ったパラメーター c の集合

H∞ := {c ∈ C | |fcn (0)| → ∞ (n → ∞)}14.4

複素力学系とジュリア集合17

図 14.6 左上：c = −1，右上：c = i, 左下：0.22 + 0.65i，右下：c = 0.25.

を考え，その補集合

M := C − H∞
*10をマンデルブロー集合 (the Mandelbrot set) と呼ぶ．

問題 14.5 (マンデルブロー集合) complexAP2 を用いて集合 M を描け．
【解答】 無限の鉢を描く方法を修正して，たとえばcolM[c_] := (p = 0.0; k = 0;
While[(Abs[p] < 2.0) && (k < 100),(p = p^2 + c; k = k + 1)];k)としてみよう．d = 0.01 はそのまま流用して，In[ ]:=In[ ]:=

tabM = Table[colM[a + b I],{a, -2, 0.6, d}, {b, -1.3, 1.3, d}];とおく．（これも実行に数十秒かかるだろう．）あとは得られたデータを上の問題で作ったcomplexAP2[tabM] のようにして表示すればよい．結果として次のような図が得られる：*10 この集合は「Kc が連結（ひとつながり）になっているような c」と一致することが知られている．

1814

力学系（補講１）■

参考：ニュートン法再訪. 方程式 f (z) = z − 1 = 0 について，複素数解もふくめて3考えよう．じつはニュートン法は複素数でもそのまま応用できて，初期値 z0 が方程式の解に十分近いとき，そのニュートン写像 Nf (z) による軌道は解に近づくことが知られている．
上で 2 次多項式の「無限の鉢」を描かせたときと同様の考え方で，「軌道が 1 の 3乗根に近づくかどうか」を判定するプログラムを描けば，有理関数 Nf (z) の「ジュリア集合」（カオス部分）を描くことができる．
たとえばこの f の場合（f, df, newton など使用済みの文字は適宜クリアしておこう），d = 0.01; f[x_] = x^3 - 1;
df[x_] = D[f[x], x];
newton[x_] := x - f[x]/df[x];colN[z_] := (p = z; k = 0;
While[(Abs[ f[p] ]> 0.1) && (k < 50),(p = newton[p]; k = k + 1)];k);
tabN = Table[colN[x + y I], {x, -2, 2, d}, {y, -2, 2, d}];complexAP2[tabN]としてみよう．次のような，意外と複雑なフラクタル図形が得られるはずである．

14.5

14.5

参考文献

参考文献

離散力学系理論と複素力学系理論の入門書として次を挙げておく．

(1) R. Devaney．「カオス力学系入門（第 2 版）」．共立出版．
(2) 川平友規．「マンデルブロー集合 — 2 次多項式の複素力学系入門」．
http://www.math.titech.ac.jp/~kawahira/courses/mandel.pdf

19