まずはやってみよう ! - 数値計算の具体例
ここでは数値計算の実例として、ロジスティック・モデルを計算してみましょう。ロジスティック・モデルの解はカオスになる場合があることが知られています。カオスは数値計算によって発見された現象です。ロジスティック・モデルの例に示されるように、非常に簡単な法則から不規則で予測不可能な状態が出現するのです。いろいろな値とともにロジスティック・モデルの解の振る舞いを数値計算によって調べ、カオスを発見してみましょう。
ロジスティック・モデルとは？
ここでいうロジスティック（logistic）とは「論理学」という意味ではありません。もともとは、兵站（へいたん＝戦争の際の物資の補給）の意味でした。恐らく必要な補給物資の量を算出するためのモデルに端を発しているのでしょうが、普通は生物や人口の増減のモデルとして説明されています。生物の増減は、実際にはいろいろな複雑な要因が絡み合うのですが、非常に単純化できたと仮定しましょう。ここでは人口の増え方としてモデルを立ててみます。
まず、人口増加率は人口に比例するとしましょう。人口の調査は毎年１回行うとして、第n回目の調査時の人口をy(n)と表すします。すると、
y(n+1) = a y(n), a>0 … (1)
aは定数で、a>1ならば人口増加、a<1ならば人口減少になります。
一方、人口があまりに増えすぎると、人口増加率が低下するような要因が働く（食料不足とか）でしょう。この効果を表すために、人口増加率の低下分も人口に比例するとします。すなわち、
a → a (1- b y) … (2)
とaを置き換えることになります。(1)、(2)より、
y(n+1) = a (1 - b y(n)) y(n) … (3)
となります。ここで、両辺にbをかけて、b y(n) → x(n) と置き換えをすると
x(n+1) = a x(n) (1 - x(n)), a>0 … (4)
となります。(4)をロジスティック方程式と呼び、x(n) から x(n+1) へのロジスティック写像と呼びます。ロジスティック写像は図1のような放物線上のリターンマップとして表すことができます。
図1: ロジスティック写像のリターンマップの例
xの初期値としては 0<x(0)<1 を考えます。0<a<1の場合は、x(n) は単調に0に近づきます。また、a>1の場合は、人口増加と人口増加率の減少の効果が釣り合って、最終的には図2のような人口一定の状態になると予想できます。（この曲線をロジスティック曲線と呼びます。）しかし、実際には a>2 の場合にはまったく予想外の振る舞いを示すようになります。この予想外の振る舞いの中にカオスを見いだすことができます。この演習では a に対して x(n) がどのような振る舞いをするのかを調ることにしましょう。
図2: ロジスティック曲線の例
ロジスティック写像のプログラム（その1）
1 program logistic1 ! プログラム名の定義：program プログラム名
2 implicit none ! 暗黙の型宣言の無効化 テキスト4.1.2 必ず入れよう
3
4 integer,parameter:: nmax=100 ! 変数宣言
5 integer n
6 real a,x
7
8 a=1.6 ! 変数の値を与えている
9 x=0.2
10
11 do n=0,nmax
12 write(*,'(i5,f10.5)') n,x ! *:画面上に出力 (i5,f10.5):書き出す変数 = 今はn,x の型
13 x=(1.0-x)*x*a ! 計算式
14 end do
15
16 end program logistic1 ! 終りの宣言 endは必ず付けること。
上のリストを見てください。このプログラムはロジスティック写像をnmax回反復して、n回目の値を出力するプログラムになっています。nmaxは4行目で定義されていて、a の値は8行目で、x(0) の値は9行目で定義されています。まずはこのプログラムを入力して、logistic1.f90 というファイル名で保存してください。
保存が終了したら、
gfortran logistic1.f90
コンパイルをしてからプログラムを走らせて見てください。1カラム目に n、2カラム目に x(n) が表示されるのが確認できましたか？
ロジスティック写像のプログラム（その2） 〜 変数の入力
1 program logistic2
2 implicit none
3
4 integer,parameter:: nmax=100
5 integer n
6 real a,x
7
8 write(*,'(a,$)') 'input a,x0 : ' ← ★ ! 実行時に変数の値、入力可。
9 read(*,*) a,x ← ★ ! その入力した変数の読み込み
10
11 do n=0,nmax
12 write(*,'(i5,f10.5)') n,x
13 x=(1.0-x)*x*a
14 end do
15
16 end program logistic2
さて、最初のプログラムでは、a と x(0) の値を固定していました。しかし、これから私たちはいろいろな a と x(0) に対して結果を調べたいので、値を変えるたびにいちいちコンパイルし直すのは面倒です。プログラムの実行時に値を設定できるように改造してみましょう。8、9行目を上のプログラムのように変更して、logistic2.f90 として保存してください。logistic2.f90をコンパイルして実行すると、値を聞いてくるので、a と x(0) を順番に入力してください。最初のプログラムと同じような出力が得られるはずです。
ロジスティック写像のプログラム（その3） 〜 ファイルへの出力
1 program logistic3
2 implicit none
3
4 integer,parameter:: nmax=100 ! 変数宣言
5 integer n
6 real a,x
7 character(len=20) fname ! 文字の宣言 : len 以下は文字数,20 wordを超えられないので注意
8
9 write(*,'(a,$)') 'input a,x0 : '
10 read(*,*) a,x
11 write(*,'(a,$)') 'input file name : ' ← ★ ! file名も後で入力可に。
12 read(*,*) fname ← ★ ! そのfile名の読み込み
13
14 open(10,file=fname) ← ★ ! fileを開く
15
16 do n=0,nmax
17 write(10,'(i5,f10.5)') n,x ← ☆ ! *→10 : 出力場所を画面からfileに変更
18 x=(1.0-x)*x*a
19 end do
20
21 close(10) ← ★ ! fileを閉じるのも忘れずに
22
23 end program logistic3
ここまでのプログラムでは、折角の計算結果が画面に流れるだけで保存ができません。計算結果を好きなファイル名で保存できるように改造してみましょう。★の行が追加で、☆の行が変更です。12行目で入力したファイル名にしたがって、14行目でファイルをオープンし、17行目の計算結果をファイルに書き込みます。計算処理が終わったら21行目でファイルをクローズすることを忘れずに。このプログラムを logistic3.f90 として保存してコンパイル／実行してみましょう。
以上で好きな a と x(0) に対して、好きなファイル名のファイルの中に計算結果を保存できるようになりました。計算結果のデータを元にして gnuplot を用いて横軸に n、縦軸に x(n) をとった時系列データをプロットしてみましょう。a の値をいろいろと変えていくと、解の性質が変わることが見えるはずです。
ロジスティック写像のプログラム（その4） 〜 いろいろな a に対する計算
1 program logistic4
2 implicit none
3
4 integer,parameter:: nmax=200,nout=100
5 integer n,i,ia
6 real a,a1,a2,da,x,x0
7
8 write(*,'(a,$)') 'input x0 : ' ! 初期値を代入
9 read(*,*) x0 ! 入力した初期値を読み込む
10
11 open(10,file='chaos.dat') ! fileを開く
12
13 a1=0.0
14 a2=4.0
15 da=0.001
16 ia=int((a2-a1)/da)
17 do i=0,ia
18 a=a1+real(i)*da ! a1 から a2 の da 毎の a について
19 x=x0
20 do n=0,nmax ! 各 n 毎に
21 if (n>nout) write(10,'(2f10.5)') a,x ! 書き出し
22 x=(1.0-x)*x*a ! 計算する
23 end do
24 end do
25
26 close(10)
27
28 end program logistic4
いろいろな a に対して解の振る舞いが変わることが確認できましたか？aに対して解の振る舞いを整理する一つの方法として、横軸に a の値を、縦軸にある繰り返し回数以上の x(n) の値をプロットしてみましょう。そのためのデータファイルを作るプログラムは上のリストになります。logistic4.f90 として入力してください。実行結果は chaos.dat というファイルに保存されます。（何度も走らせると上書きされてしまうので注意！）
○ 課題
1. logistic4.f90 に示したプログラムを実行して、a の値によってロジスティック写像の振る舞いが変化することが見えたと思います。この結果をgnuplotでプロットしてください。次に、この図から読みとれるいくつかの特徴的な aを選んでください。それぞれのaの場合について、logistic3.f90を用いて n に対する x(n) の振る舞いを計算し、その結果をgnuplotでプロットしてください。それらの図を確認して、ロジスティック写像の性質についてまとめください。
2. カオスの特徴の1つとして、初期値に対する敏感な依存性があります。カオス的な振る舞いを示す a を適当に選び、同じ a のもとで非常に近い x(0) を選んで（例えば0.3 と 0.30001など）計算を行ってみましょう。ほんの少しだけ離れた x(0) であってもn が大きくなるにつれてx(n)が次第に異なった値を示す様子をlogistic3.f90を用いて計算し、結果を gnuplot でプロットしてください。
gnuplotは、関数のプロットや数値データのプロットなどを行う、非常に すぐれたフリーソフトです。多くの数学関数を標準 に持ち、解りやすいコマンドを使って対話的にグラフを作成することができます。簡単な使い方はこちら。もっと本格的に使いたい場合はこちら。