カオス力学系の数値計算：カオスを目で見てみる(1次元離散時間力学系)

はじめに

先日、書店の理工系図書のコーナーに行ったら劉慈欣による小説『三体』のフェアでカオスに関する書籍が並べられていました。筆者がカオスというものを知り、詳しく調べ始めて以来、このように注目を集めることはなかったように思います。もしかしたら小説『三体』の流行によってカオスに興味を持ってくれる人たちが増えるのではなかろうかと密かに期待しています。

カオスに興味を持ってくれた人たち、その中でもとくに学生やエンジニアのために筆者がQiitaでできることのひとつはカオスに関する数値計算のコードを紹介し解説することだと考えています。そのために、これから一連の記事として「カオス力学系の数値計算」シリーズを執筆し公開したいと思っています。記事の中には大学レベルの数学が出てくることもありますがなるべく平易な言葉で書くようにするつもりです。読者の方にとって得るものがあれば嬉しいです。

カオスという言葉はすでに広く知れ渡っていて、多くの人々に日常的に使われる言葉になっています。普段使われている意味はおおよそ「めちゃくちゃでよくわからないさま」だと思います。もともとカオスは力学系理論における学術用語ですが、誤解を恐れずに言えば学術用語としてのカオスも大体似たような意味です。ただし、もう少し限定された対象に使われますし、カオスかどうかを判定するための定量的な指標もいくつかあります。本記事を読めば、本来の意味でのカオスがわかるようにしたいと思います。

一連の記事では通してPython 3.x系を使う想定です。コードを簡単にし計算を高速にするため計算やプロットにいくつかモジュールを使用します。モジュールはコードの中に示します。コードは必ずしも最適とは言えないかもしれませんが、最低限動くものを載せるつもりです。本記事の内容・コードが読者の方々のお役に立てば幸いです。

導入

力学系理論は時間とともに変化する系を調べる数学の一分野です。系のふるまいは微分方程式や差分方程式で表現されることが多いです。（微分方程式・差分方程式をよく知らなくてもこの先の内容はわかると思いますので構えずに読み進めていただければと思います。）力学や時間という物理の言葉が出てくるのは力学系理論が力学における質点系の運動の解析に端を発するためだと思われます。物理学における力学系は19世紀・20世紀初頭にポアンカレらによって研究されました。余談ですがポアンカレは天体の運動を調べる天体力学における三体問題（小説『三体』の名前の由来となったそうです）に大きな寄与を残したエライ人です。数学の位相幾何（トポロジー）における有名な問題「ポアンカレ予想」で聞いたことがある人もいるかもしれません。

時間とともに変化する系は科学のさまざまな分野で現れます。その意味で力学系は物理学のみならずあらゆる科学との関連が非常に深い分野と言えます。20世紀も中盤を過ぎて電子計算機による数値計算が発展すると力学系の研究はさらに飛躍しカオスの発見につながっていきます。カオスの発見につながる結果は気象学や数理生物学で現れました。このあたりの歴史はたとえば山口昌哉著『カオスとフラクタル』に詳しく書かれています。

その発見以来、カオスの研究には数値計算がよく使われています。カオス発生下では系のふるまいを簡単な関数で記述できない場合が多く、そのような場合にふるまいを見るには数値計算を使わざるを得ない事情があります。カオスは非常に単純な力学系でも生じるので簡単な数値計算を行えば目で見ることができます。プログラミングの初心者でも簡単に計算できるので楽しんで学ぶことができるでしょう。ただし、カオスが発生する力学系の数値計算を甘く見てはなりません。典型的なカオスでは誤差の指数関数的増大が見られるためです。粗っぽく言えば、指数関数的に誤差が増大するということは実装者が本来計算して求めたかった結果ととんでもなくずれた結果が数値計算によって得られてしまうということです。しかもコンピュータでは数値はかならず有限のビット数で表さなければなりませんから、無理数のような小数点以下に無限の桁を持つ数はかならず誤差含みの表現をしなければなりません。無理数の例のみならず数値計算をする以上、誤差から逃れることはできないと言えるでしょう。指数関数的な誤差増大という特徴を持つカオスを数値計算で調べるのは危ないのですが、数値計算を使わないとカオスを調べるのは難しいというジレンマのような状態になっています。（ここにカオス研究の大変さがあります。）

筆者は学生でもエンジニアでも、数値計算をする人はみんなカオスのことを知っておいて損はないと考えています。前の段落で述べたように、カオスは非常に単純な力学系でも生じるのにもかかわらず、しばしば指数関数的な誤差増大という特徴を示します。実装者は単純な微分方程式や差分方程式を計算しているつもりなのに知らず知らずのうちにカオスが発生する力学系の計算をしていて誤った数値計算の結果を得てしまうことがあるかもしれません。とくに微分方程式を離散化したり単純化したりする場合は要注意です。微分方程式の解が規則的なふるまいを示していても、離散化・単純化して差分方程式にした途端にカオスが生じることはしばしばあります。

本記事の構成は以下のとおりです。最初に力学系そのもの、および一連の記事の中で共通して現れる用語を説明します。過不足なく説明することができればよいのですが本記事を執筆している段階では一連の記事で書くことをすべて決めていないので過不足があるかもしれません。それらは今後公開される記事の中で適宜説明したいと思います。次にもっともシンプルなカオス力学系としてロジスティック写像を例にとり、まずはカオスを目で見るということをしたいと思います。本記事ではただ見るということにとどめ、数学的な説明は今後の記事で書くこととします。

力学系の用語の説明

物体の運動などの興味のある対象の、ある時刻における状態を適当な空間の1点で表すことを考えます。解析力学を知っている人は質点系の運動状態が位置と運動量の空間上の1点で表せたことを思い出してください。力学系の分野では状態を表す空間のことを相空間と呼びます。ここでは相空間を$X$で表すこととします。

状態の移り変わりは相空間$X$上のある点から同じ空間$X$上の別の点への遷移として表現されます。対象の状態変化が何らかの一定の法則$f$にしたがう状況を考えましょう。法則$f$の性質によって力学系は離散時間力学系と連続時間力学系に分かれます。

本記事では離散時間力学系のみを定義します。一言で表すと離散時間力学系とは定義域と終域とが同じ空間であるような写像のことです。写像とは定義域の点を定めると終域の点が一点に定まるような変換のことです。誤解の恐れがない場合には離散時間力学系を指して写像と呼ぶこともあります。一般に写像$f$の定義域を$X$とし終域を$Y$とするとき、$f:X \to Y$のように書き表します。離散時間力学系は相空間$X$から同じ空間$X$への写像なので$f:X \to X$と表されます。

写像$f$によって点$x$が点$y$に移されることを$y=f(x)$や$f:x \mapsto y$のように書き表します。写像$f:X \to X$は定義域と終域が一致しているので写像$f$で移した点をさらに写像$f$で移すことができます。点$x$に対して写像$f$を$n$回繰り返し適用すると$f(f(f(\dots f(x)\dots )))$となってしまって長くなるので簡略化して$f^n(x)$で表します。$f^n$もまた空間$X$から$X$への写像となります。

以降では初期状態を表す相空間上の点には添え字をつけて$x_0$と表します。点$x_0$に写像$f$を$n$回繰り返し作用させて移した点を$x_n$で表します。初期点$x_0$とそれに写像$f$を繰り返して得られる点を並べた点列${ x_0, x_1, \dots, }$を点$x_0$を通る軌道と呼びます。（正確には正の半軌道と呼びます。）力学系の研究ではある初期点が与えられたときに軌道がどのようなふるまいを示すかにしばしば注目します。

カオスを目で見る

相空間として実数の区間$[0,1]$を考えます。写像$f:[0,1] \to [0,1]$のように1次元空間上で定義される写像は離散時間力学系としてもっとも簡単なものであり、よく調べられています。区間上で定義されているため区間力学系と呼ばれることもあります。

カオスが発生する区間力学系として非常に有名なロジスティック写像（もしくは2次写像と呼ばれる）は以下で定義されます¹。
$$
f_a(x) = ax(1-x)
$$
ここで$a$は実数のパラメータであり、本記事では$0<a \leq 4$とします。

ロジスティック写像における軌道の様子を可視化してみましょう。1次元の区間上を点が動くので単純に描いたら線分上にただ点が並ぶだけで、どの点からどの点に移ったのかがまったくわからなくなってしまいます。そこで縦軸に$x_{n+1}$、横軸に$x_{n}$をとり、すべての$n$について同じ平面上に点をプロットしてみます。$a=4$の場合にプロットした図を図1として載せました。$x_{n+1}=f_a(x_n)$ですから、点は$f_a(x_n)$上にプロットされます。どの点がどの点に移されたかはわかるようになりましたが、点の長時間の変遷はパッと見ではわかりません。

図1. ロジスティック写像の軌道を縦軸$x_{n+1}$・横軸$x_{n}$としてプロットしたもの($a=4$)

前の段落で実施したプロットでは写像を$N$回繰り返したときの軌道を表したとき点が打たれる位置は$(x_0,x_1),(x_1,x_2), \dots , (x_{N-2}, x_{N-1}), (x_{N-1}, x_{N})$です。縦軸の方の座標の値が、次の横軸の方の座標の値になっていることが見て取れます。この性質を利用して点の変遷を表すための補助線を引きましょう。そのためには次のようにします。

1. 関数$f_a(x_n)$と直線$x_{n+1}=x_n$を描きます
2. 直線$x_{n+1}=x_n$上の点$(x_0,x_0)$の位置から関数$f_a(x_n)$のグラフ上の点$(x_0,x_1=f_a(x_0))$に向かって補助線を引きます
3. 点$(x_0,x_1)$から直線$x_{n+1}=x_n$に向かって補助線を引きます
4. 3で引いた補助線と直線$x_{n+1}=x_n$の交点から関数$f_a(x_n)$のグラフに向かって補助線を引きます
5. 3,4を$N$回分、繰り返します

この方法はグラフ反復法と呼ばれます²。グラフ反復法は1次元離散時間力学系で軌道を可視化する際によく使われます。

ロジスティック写像における軌道をグラフ反復法で作図するコードを以下に記載します。今後使用することも考えて写像の定義はモジュール化しています。

one_dim_maps.py

# ロジスティック写像: f(x) = ax(1-x)
def logistic_map(a:float,x:float) -> float:
	return a * x * (1. - x)

graph_repetition.py

import one_dim_maps as maps
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
mpl.rcParams['text.usetex'] = True

# 初期点xを通るロジスティック写像の軌道を得るメソッド
def get_orbit(f, a, x, N):
    orbit = np.ones(N) * np.nan
    for i in range(N):
        orbit[i] = x
        x = f(a, x)
    return orbit

# グラフ反復法で軌道を描画するメソッド
def plot_graph_repetition(f, orbit):
    # 図の設定
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.set_aspect(1.0 / ax.get_data_ratio())
    ax.set_xlim(0,1)
    ax.set_ylim(0,1)
    ax.tick_params(labelsize=15)
    ax.set_xlabel("$x_n$", fontsize = 30)
    ax.set_ylabel("$x_{n+1}$", fontsize = 30)
    
    # 手順1：fのグラフと直線x_{n+1}=x_nをプロット
    X = np.linspace(0., 1., 10**3)
    Y = np.array([f(a, x) for x in X])
    ax.plot(X, Y, '-', color='k')
    ax.plot(X, X, '-', color='k')

    # 手順2および手順3,4の反復
    for i in range(orbit.size - 1):
        ax.plot([orbit[i], orbit[i]], [orbit[i], orbit[i + 1]], '-', color = 'k')
        ax.plot([orbit[i], orbit[i + 1]], [orbit[i + 1], orbit[i+1]], '-', color = 'k')
    
    plt.tight_layout()
    plt.show()

x = np.random.random()
a = 2. #4.
N = 100
orbit = get_orbit(maps.logistic_map, a, x, N)
plot_graph(maps.logistic_map, a, orbit)

上記のコードでロジスティック写像における軌道をプロットした図が以下の図2です。
パラメータ$a=2$として99回反復しました。パラメータ$a=2$の場合はカオスは発生せず、軌道は1点に吸い込まれていくことが知られています³。ちなみにグラフ反復法で描いた図をしばしば蜘蛛の巣ダイアグラムと呼びますが²、図2の場合は蜘蛛の巣感がありません。

図2. ロジスティック写像の蜘蛛の巣ダイアグラム($a=2$)

軌道を時系列としても見てみましょう³。図3は図2で示した軌道を時系列にしたものです。点と点を結ぶ線はただの補助線です。時系列として見た場合でも、1点に吸い込まれたあとは点が動かない様子を観察できます。

図3. ロジスティック写像の軌道を時系列として見たもの($a=2$)

次にパラメータの値を変えて$a=4$にしてみましょう。初期点は図2と同じところにとって99回反復させた図が以下の図4です。図2のときと比べると複雑なグラフが描画されました。こちらはかなり蜘蛛の巣感があります。

図4. ロジスティック写像の蜘蛛の巣ダイアグラム($a=4$)

図4で示した軌道も時系列として見てみましょう。図5に図4の軌道を時系列にしたものを示します。図3の場合と異なり、1点に吸い込まれることなく振動しています。振動は規則的とは言えず、ランダムな様相を示しています。

図5. ロジスティック写像の軌道を時系列として見たもの($a=4$)

ロジスティック写像$f_a$はパラメータ$a=4$のときカオス力学系になることが知られています²。$a=4$のときのロジスティック写像$f_a$が生成する軌道は図4、5で示したように複雑なふるまいを示します。この複雑なふるまいは長期的な予測が不可能であり、カオスと呼ばれています。

以上で本記事の目的であるカオスを目で見てみることは一応、達成できたと思います。今後はロジスティック写像についてさらに詳しいことや2次元のカオス力学系などについて書きたいと思っています。

※タイトルとコードに誤りがあったので修正しました(2024.06.11)。

参考文献

引用はしていませんが本記事を書くにあたって振り返った文献を記載します。
本文中で引用した文献は脚注に載せています。

• 山口 昌哉, 『カオスとフラクタル』, ちくま学芸文庫 (2010).
• 国府 寛司, 『力学系の基礎』, 朝倉書店 (2000).
• 久保 泉, 『力学系1』, 岩波書店 (1997).

1. R. L. Devaney, 『新訂版 カオス力学系入門 第2版』, 共立出版 (2003). ↩

2. M. W. Hirsch, S. Smale and R. L. Devaney, 『力学系入門 原著第2版』, 共立出版 (2007). ↩ ↩² ↩³

3. 小室 元政, 『新版 基礎からの力学系』, サイエンス社 (2005). ↩ ↩²