カオス・フラクタル 講義ノート #7
担当 : 井上 純一 (情報科学研究科棟 8-13)
URL : http://chaosweb.complex.eng.hokudai.ac.jp/˜j inoue/index.html

平成 23 年 6 月 14 日

目次
10 数値計算の準備: 常微分方程式の数値解法

60

10.1 テーラー展開による漸化式の導出 . . . . . . . . . . . . . . . . . . . . . . . . . . .
10.2 ルンゲ・クッタ法: 基本的アイデア . . . . . . . . . . . . . . . . . . . . . . . . . .
10.2.1 適用例 1: ロジスティック方程式 . . . . . . . . . . . . . . . . . . . . . . . .

61616210.2.2 適用例 2: レスラー方程式 . . . . . . . . . . . . . . . . . . . . . . . . . . .

6410

数値計算の準備: 常微分方程式の数値解法

今までこの講義では 1 次元写像の振る舞いを詳しくみてきた. 例えば, ロジスティック写像は次のロジスティック方程式:

dN
dt= λN (1 − µN )(139)を差分化した方程式であり, 本来解くべき方程式の「近似」である. しかし, 現実の自然現象は運動方程式など微分方程式で記述されるものが多く, 従って, これらの現象を正確に調べようとすれば,差分化された近似式ではなく, 微分方程式そのものを扱わなければならない. いくつかの常微分方程式は解析的に解くことができるが, 非線形項が入り, しかも, 変数が複数絡み合った連立微分方程式は多くの場合に手で解くことはできない. このような手で解くことのできない方程式の多くはカオスなどの複雑で多彩な挙動を示し, 従って, そのような現象の理解を深めることを目的の一つとする本講義では, そのような方程式を解くことを避けて通ることはできない (逆に言えば, 手できれいに解ける方程式は本講義の対象とならない).
そこで, 今回はやや本題を外れることになるが, 複雑な挙動を示す常微分方程式を解析するためのほぼ唯一の手段になる数値計算技法の典型例についてみていく (一部の人たちにとっては復習になるかもしれないが). 現在, 微分方程式を数値的に解くための様々なアプリケーションソフトやパッケージが利用できるが, はじめからこれらを「ブラックボックス」として用いるのではなく, 一度はその原理を学んでおくことは重要である.
具体的に我々が解こうとする方程式の典型例は次のようなものである.

dxdt
dy
dt

= −y − z(140)= x + ay(141)60

カオス・フラクタル 2011 #6

担当: 大学院情報科学研究科 井上 純一

dzdt= b + z(x − c)

(142)

この連立微分方程式はレスラー方程式と呼ばれる. ここに現れるパラメータ a, b, c を適切に選んだ
際, 我々が 1 次元写像でみたようなカオスが生じるかどうか, どのような力学的性質を持つかなどを調べることがここからの目標である. 今回の講義ではそのための準備を行う.

10.1

テーラー展開による漸化式の導出

話を簡単にするために次の 1 変数の微分方程式:

dydt= f (t, y)(143)

を考えよう. ただし, 関数 f は微分可能な関数であるとする. ここで h を十分小さい量としtn= t0 + nh, n = 1, 2, 3, · · ·(144)

で「時間」t を離散化すると (t0 は時間原点), tn+1 = t0 + (n + 1)h であるから, tn+1 − tn = h  1となることに注意しょう (あるいは tn+1 = tn + h). このとき, 小さな量 h に関するテーラー展開を思い出すとy(tn+1 )= y(tn + h) = y(tn ) +dydth + O(h2 )tn ,yn2= y(tn ) + f (tn , yn )h + O(h )(145)となる. 従って, y(tn ) = yn , y(tn+1 ) = yn+1 と略記することを約束すれば

yn+1
tn

= yn + f (tn , yn )h + O(h2 )(146)= t0 + nh(147)が得られる. この方程式 (146) を漸化式として解き, 各 n で (147) 式で与えられる tn を yn に対す
る時刻とみなすことで微分方程式を数値的に解く方法をオイラー法と呼ぶ. 上記から明らかなように, 各ステップでの所謂「局所誤差」は h2 のオーダである (O(h2 ) と表記した部分).

10.2

ルンゲ・クッタ法: 基本的アイデア

上記のテーラー展開で微分係数 dy/dt = f (t, y) を評価する点は (tn , yn ) であったが, これを(tn , yn ) と (tn+1 , yn+1 ) の中点である (tn+1/2 , yn+1/2 ) に選んだ場合, 上記のテーラー展開のプロセ
スと結果として得られる漸化式はそれぞれの区間でどうなるかを考えてみよう.
まず, (144) 式から tn+1/2 = t0 + (n + 1/2)h であるからtn+1= tn+1/2 + h/2, tn = tn+1/2 − h/2(148)であることに注意して

yn+1 = y(tn+1 ) =y(tn+1/2 + h/2)

= yn+1/2 + f (tn+1/2 , yn+1/2 )( )( )2h1 0h+ f (tn+1/2 , yn+1/2 )+ O(h3 )222(149)
ここは 61 ページ目

カオス・フラクタル 2011 #6

担当: 大学院情報科学研究科 井上 純一

およびyn = y(tn )= y(tn+1/2 − h/2)
= yn+1/2 − f (tn+1/2 , yn+1/2 )

( )
( )2
h1 0
h+ f (tn+1/2 , yn+1/2 )+ O(h3 )(150)
2220が得られる. ここで, f ≡ df /dt と略記したことに注意されたい. 従って, これらの方程式 (149)(150)を辺々引くと

yn+1 − yn= f (tn+1/2 , yn+1/2 )h + O(h3 )(151)となりyn+1/2= y(tn+1/2 ) = y(tn + h/2) = y(tn ) +hk1
f (tn , yn ) = yn +22(152)であるから, 問題の微分方程式 (143) をyn+1=yn + k2 + O(h3 )

(153)

k1

=

hf (tn , yn )

(154)

k2

=

hf (tn + h/2, yn + k1 /2)

(155)

tn

=

t0 + nh(156)の漸化式で書き直すことで, ステップごとの局所誤差をオイラー法での O(h2 ) から O(h3 ) へまで低減させることができることがわかる1 . これら (153)(154)(155)(156) 式を 2 次のルンゲ・クッタ法と呼ぶ.
このような方法でうまく中継点を選び, 上記の手続きを繰り返すことにより, 局所誤差を系統的に減らしていくことができることは容易に推測できるであろう. 導出が煩雑なので詳細は省略する
が, 局所誤差が O(h5 ) まで低減できるものとして, 次の 4 次のルンゲ・クッタ法が知られており, 多くの場面で用いられている (本講義でも数値計算では専らこの 4 次のルンゲ・クッタ法を用いる).

k1(k1 + 2k2 + 2k3 + k4 )+ O(h5 )6= hf (tn , yn )k2

= hf (tn + h/2, yn + k1 /2)

(159)

k3

= hf (tn + h/2, yn + k2 /2)

(160)

k4

= hf (tn + h, yn + k3 )

(161)

tn

= t0 + nh

(162)

yn+1

= yn +(157)(158)

次で具体的なコーディングについてみていこう.

10.2.1

適用例 1: ロジスティック方程式

以上が一般論であるが, 感じをつかむためにいくつかの具体的な方程式への適用例を見ておこう.

1h は十分に小さな量であったことを思い出す. 例えば h = 0.01 などとして, O(h2 ) と O(h3 ) を比較してみるとよい.

ここは 62 ページ目

カオス・フラクタル 2011 #6

担当: 大学院情報科学研究科 井上 純一

まずは最も簡単な例の一つとして, 上記の 4 次のルンゲ・クッタ法 (157)-(162) でロジスティック方程式:

dN
dt= λN (1 − µN )(163)を数値的に解いてみる. その結果を図 35 に載せよう. この図より, 解析解と極めて良好な一致を見せていることがわかる2 . また, そのプログラムを載せる.
20

N0 = 2, λ = 0.5, µ = 0.1

181/(0.1-(0.1-0.5)*exp(-0.5*x))N0 = 10, λ = 0.1, µ = 1

161/(1-(1-0.1)*exp(-0.1*x))N0 = 20, λ = 0.5, µ = 0.11/(0.1-(0.1-0.05)*exp(-0.5*x))

14N

121086420024t6810

図 35: ロジスティック方程式の 4 次のルンゲ・クッタ法に基づく数値解と解析解.

/* ルンゲ・クッタ法のロジスティック方程式への適用例 */#include <stdio.h>#include <math.h>#define h 0.01
#define mu 0.1#define lambda 0.5
double func(y)double y;{return (lambda*y*(1.0-mu*y));}/* 関数 f の定義 */main(){FILE *fpr;
int i,imax=1000;
double y,k1,k2,k3,k4,k;if((fpr = fopen("rk4.dat", "wt")) !=NULL){for(i = 0,y=20.0; i <= imax; i++){k1 = h*func(y);k2 = h*func(y+0.5*k1);
2 大域誤差 (真の解 y(t) と yn の差) についての評価は後の演習などで触れたいと思う.

ここは 63 ページ目

カオス・フラクタル 2011 #6

担当: 大学院情報科学研究科 井上 純一

k3 = h*func(y+0.5*k2);
k4 = h*func(y+k3);
k = (k1+2.0*k2+2.0*k3+k4)/6.0;
y = y + k;
fprintf(fpr,"%lf %lf\n",i*h,y);}}fclose(fpr);}

10.2.2

適用例 2: レスラー方程式

変数が 3 つに増えたレスラー方程式 (140)(141)(142) も上記 1 変数の場合のプログラムを拡張することで容易に解くことができる. 参考までにプログラミングのソースコードを載せておく.

/* ルンゲ・クッタ法のレスラー方程式への適用例 */#include <stdio.h>#include <math.h>#include <stdlib.h>#define h 0.01
#define a 0.5#define b 0.5#define c 5.7
double func1(x,y,z) /* 関数その 1 */double x;
double y;
double z;{return (-(y+z));}double func2(x,y,z) /* 関数その 2 */double x;
double y;
double z;{return (x+a*y);}double func3(x,y,z) /* 関数その 3 */double x;
double y;double z;{return (b+z*(x-c));}main(){FILE *fpr;
int i,imax=50000;

ここは 64 ページ目カオス・フラクタル 2011 #6

担当: 大学院情報科学研究科 井上 純一

double x,k1x,k2x,k3x,k4x,kx,y,k1y,k2y,k3y,k4y,ky,z,k1z,k2z,k3z,k4z,kz;
if((fpr = fopen("rossler.dat", "wt")) !=NULL){
for(i = 0,x=-0.01,y=-0.01,z=0.01; i <= imax; i++){k1x = h*func1(x,y,z);k1y = h*func2(x,y,z);
k1z = h*func3(x,y,z);k2x = h*func1(x+0.5*k1x,y+0.5*k1y,z+0.5*k1z);
k2y = h*func2(x+0.5*k1x,y+0.5*k1y,z+0.5*k1z);
k2z = h*func2(x+0.5*k1x,y+0.5*k1y,z+0.5*k1z);
k3x = h*func1(x+0.5*k2x,y+0.5*k2y,z+0.5*k2z);
k3y = h*func2(x+0.5*k2x,y+0.5*k2y,z+0.5*k2z);
k3z = h*func3(x+0.5*k2x,y+0.5*k2y,z+0.5*k2z);
k4x = h*func1(x+k3x,y+k3y,z+k3z);
k4y = h*func2(x+k3x,y+k3y,z+k3z);
k4z = h*func3(x+k3x,y+k3y,z+k3z);
kx = (k1x+2.0*k2x+2.0*k3x+k4x)/6.0;
ky = (k1y+2.0*k2y+2.0*k3y+k4y)/6.0;
kz = (k1z+2.0*k2z+2.0*k3z+k4z)/6.0;
x = x + kx;
y = y + ky;z = z + kz;
fprintf(fpr,"%lf %lf %lf %lf\n",i*h,x,y,z);}}
fclose(fpr);}
上記のプログラムを用いて, まずは, 3 つの変数 x, y, z の時間変化をプロットすると図 36 のようになる. この図 36 より, 1 次元写像に見たような複雑な挙動を x, y, z のそれぞれが示すことがわかる.
2060x(t)y(t)

15501040530020-510-100-15050100t

150200-10z(t)0100200t

300400500

図 36: レスラー方程式の 4 次のルンゲ・クッタ法に基づく数値解. x(t), y(t) および z(t) の振る舞い.
この図 36 を見方を変えて, xyz の 3 次元空間内にその軌道をプロットすると興味深い. これは上記
ここは 65 ページ目

カオス・フラクタル 2011 #6

担当: 大学院情報科学研究科 井上 純一

プログラムでファイルに書き出されたデータの第 2 軸を x 軸, 第 3 軸を y 軸, 第 4 軸を z 軸に選んだプロットに相当する.
ちなみに, 「計算機プログラミング I・同演習」や実験などで gnuplot を学んだものは gnuplot を起動したのちに

gnuplot> set parametricgnuplot> splot ’lossler.dat’ u 2:3:4 w lineとすればよい (先頭の gnuplot> は gnuplot のカーソルであることに注意). このようにして 3 次元プロットした結果を図 37 に載せる3 . 軌道の詳細を調べるため, 図 37 の 3 次元プロットを x-y 平

a = b = 0.5, c = 5.7z6050403020100-10208641520y10-2-45-6-80-10-12x-5-14 -10

図 37: レスラー方程式の数値解から得られる軌道の様子. パラメータは a = b = 0.5, c = 5.7 に選んである.
面, x-z 平面に射影して再プロットしたものが図 38 である.

以上で本講義の [カオス編] 後半で必要な数値計算上の道具立ては全て揃ったことになる. これらの図 37,38 の意味するところは次回以降に詳しくみていくことにしよう.

3 ここで出てきた gnuplot の使用法などは演習の時間に取り上げたいと思う.

ここは 66 ページ目

カオス・フラクタル 2011 #6

担当: 大学院情報科学研究科 井上 純一

86065042400y30-2z-420-610-8-100-12-14-10-505x101520-10-10-505x

101520

図 38: 図 37 の 3 次元プロットを x-y 平面, x-z 平面に射影してプロットしたもの
レポート課題 6
次の x に関する 2 階の常微分方程式:

d2 xdx− (1 − x2 )+x = 02dt
dt(164)は速度を dx/dt = v で定義すれば

dv
dtdxdt= (1 − x2 )v − x(165)= v(166)と x, v に関する 1 階の連立微分方程式に書き直すことができ4 , これは今回学んだルンゲ・クッタ法を用いて解くことができる. そのプログラムを書き (ソースコードの提出), いくつかの  の値, 初期条件に対する x, v の振る舞いをプロットせよ (gnuplot などの作図ツールを知らないものはソースコードのみでよい).

(連絡) : 次週 (6/21) は担当者 (井上) が午後 3 時より札幌北会場 WING 進学相談会（ダイヤ書房）
シャトレーゼガトーキングダムサッポロブースにて, 「北大工学部進学相談会」を実施するため, 通常より早めに終わる可能性があります. この講義で習得すべき内容が終わらない可能性が出てきた
場合, 1 回程度の補講を行う場合があることを了解ください.

4 多くの運動方程式はこの手の書き換えができることに注意しょう.

ここは 67 ページ目