反応拡散系に対する分岐計算ソフトウェア

RD-AUTO の使用マニュアル Ver.0.98
林悠帆（金沢大学大学院自然科学研究科）

田中直毅（金沢大学理工学域数物科学類）

樋口亮（金沢大学大学院自然科学研究科）

若井健（金沢大学大学院自然科学研究科）

長山雅晴（北海道大学電子科学研究所）
2018 年 12 月 17 日

マニュアル更新履歴
年月日

更新者

2010.03.20
2011.03.22

林悠帆
田中直毅

2011.10.17
2012.12.20
2013.08.10
2013.08.10
2014.12.13
今後の予定

樋口亮
若井健
長山雅晴
長山雅晴
長山雅晴

更新内容
マニュアルの大枠作成． デモ 1 作成．
第２章の詳しい説明の追加．デモ 2 作成．
鳥瞰図の節を作成．他．
デモ 1，デモ 2 の修正．デモ 3 作成．
デモの不具合を修正．Cross-Diﬀusion の定常解計算プログラム作成
本ソフトウエアの正式名称を「RD-AUTO」と決定
公開準備のためデモ 2 を削除
応用数学勉強会のために修正
Cross-diﬀusion のデモ作成
双安定反応拡散系のデモ作成
3 変数反応拡散系のデモ作成
空間 2 次元定常解とその安定性のデモ作成

1目次

第 1 章 はじめに

4

第 2 章 RD-AUTO のインストール
2.1 この章で扱う記号について . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 必要な環境 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.1 GLSC のインストール . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.2 GNUPLOT のインストール . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.3 MPICH のインストール . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.4 RSA 鍵を使った MPI のプロセス通信設定 . . . . . . . . . . . . . . . . . .
2.3 RD-AUTO の解凍 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3.1 Mac OSX を使って RD-AUTO を使う方法 . . . . . . . . . . . . . . . . .

5
5
5
5
5
7
9
10
11

第 3 章 発展方程式
3.1 使い方 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.1.1 時間発展方程式のプログラムに対するフローチャート . . . . . . . . . . . .
3.1.2 実行前の設定 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.1.3 実行コマンドの説明 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.1.4 出力されるデータファイル . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.1.5 gnuplot の設定 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.1.6 周期解の周期の求め方 . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.1.7 データカット . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 付録 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.1 プログラムで扱う変数について . . . . . . . . . . . . . . . . . . . . . . . .
3.2.2 反応拡散系の変数を変更する . . . . . . . . . . . . . . . . . . . . . . . . .

12
12
12
13
20
22
25
27
30
32
32
32

第 4 章 分岐計算プログラム
4.1 使い方 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.1.1 実行前の設定 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.1.2 実行コマンドの説明 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.1.3 解の枝の安定性を求める方法 . . . . . . . . . . . . . . . . . . . . . . . . .
4.1.4 出力されるデータファイル . . . . . . . . . . . . . . . . . . . . . . . . . . .

34
34
34
41
42
43

第 5 章 デモ 1
5.1 1 次元 Gray-Scott モデルの時間発展方程式 . . . . . . . . . . . . . . . . . . . . . .
5.2 1 次元 Gray-Scott モデルの分岐計算 . . . . . . . . . . . . . . . . . . . . . . . . . .

46
48
57

第 6 章 デモ 2
103
6.1 時間発展方程式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
6.2 1 次元発熱反応モデルの分岐計算 . . . . . . . . . . . . . . . . . . . . . . . . . . . 114

2

第 7 章 付録
153
7.1 鳥瞰図 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153
7.2 実行時に起きるエラーの対処法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 157

3

第1章

はじめに

本マニュアルは，次の空間 1 次元 n 変数反応拡散系に対する，定常解と周期解の分岐構造を数値計算するソフトウェア (RD-AUTO) の解説書です:


∂2u
∂u
= D 2 + f (u),∂t∂xt > 0,x ∈ (0, L),

ただし，


d1



D=


0
..

.

0



,




di ≥ 0(i = 1, 2, . . . , n),

dn







u1
 . 
n
. 
u=
 . ∈R ,

f1
.
n
.
f =
 . ∈R .

unfn※ di が全て 0 のとき方程式は常微分方程式になるが，ソフトウェアは動きます．ただし，常微分方程式を解く場合は AUTO の使用を推奨します．
RD-AUTO は Newton 法により定常解と周期解を求めるため，Newton 法を用いるための適切な初期値が必要です．その初期値は時間発展方程式を解いて求めるため，本マニュアルでは時間発展方程式を解くプログラムの使い方も説明します．自ら適切な初期値を与えることができる場合は，
本マニュアルで解説する時間発展方程式を使って初期値を与える必要はありませんが，RD-AUTOに受け渡すデータ形式に合わせる必要がありますので，注意してください．なお，データ形式についての解説は現時点1 で解説していません．RD-AUTO や時間発展方程式を解くプログラムの使い方を説明するために，例として 2 変数反応拡散系に現れるパルス波の分岐構造を求めいています．
そのためデモはすべて周期境界条件下で計算していますが，Dirichlet 境界条件や Neumann 境界
条件でも計算可能です．本マニュアルは次のような構成となっています．
第 2 章では，時間発展方程式と RD-AUTO のプログラムを動かすために必要な環境設定を解説します．
第 3 章では，時間発展方程式を解くプログラムの説明を行います．
第 4 章では，RD-AUTO の説明を行います．数値計算法の具体的な解説は後日公開します．
第 5 章，第 6 章では，具体的な方程式を用いて使い方の説明を行います．

1

2015 年 8 月 11 日

4

第2章

RD-AUTO のインストール

Intel 製の有料コンパイラである icc，ifort が使えることを前提としていますが，gcc, gfortan を用いて RD-AUTO を動かすことも可能です．周期解の数値計算では計算量が多くなるのである程度の計算機性能が必要です．ノート PC を使う場合は定常解まで計算できればよいと思ってくだ
さい．

この章で扱う記号について

2.1

最初に，この章で扱うコマンドプロンプトの記号の説明を行う．
$ のあとに書かれているコマンドは一般ユーザで入力する．
# のあとに書かれているコマンドはスーパーユーザで入力する．
斜体で version と書かれているところには，そのソフトウェアのバージョンの数字が入る．
枠で囲んであるコマンドは必ず入力するものである．
枠で囲んでいないコマンドは必ずしも入力する必要のないものである．

2.2

必要な環境

2.2.1

GLSC のインストール

鳥瞰図を書くソフトウエアとして GLSC を使うことにする場合はインストールしてください.

2.2.2

GNUPLOT のインストール

RD-AUTO ではバージョン 4.0 以降を必要とするので，Version が古い場合 Update してくだ
さい．
1. gnuplot のダウンロード先
以下のサイトから最新の gnuplot のバージョン (gnuplot- version.tar.gz ) を入手する．

http://www.gnuplot.info/
2. gnuplot のインストール/usr/local/gnuplot にインストールする場合を例とする．
gnuplot- version.tar.gz を展開する．
$ tar -zxvf gnuplot- version.tar.gz


展開後，gnuplot- version に移動し

5












$ ./configure -prefix=/usr/local/gnuplot$ make $ su  (スーパーユーザーのホームに移動した場合は，gnuplot- version まで移動
する．)
# make install # exit 







を行う．

3. パス1 を通す
使用しているシェル2 の任意の行に

setenv PATH ${PATH}:/usr/local/gnuplot/




と記述．
ただし，上記は tcsh を使用している場合の例であり，シェルによって書き方は異なる．
自分が使用しているシェルが分からない場合は，以下のコマンドで使用しているシェルが分かる．
$ echo $SHELL 
シェルを emacs で開きたいときは，自分のホームで以下のコマンドを使う．(シェルに tcshを使っている場合)
$ emacs .tcshrc注 1 または$ emacs .tcshrc.mine注 1 (注 1) 自分のホームに.tcshrc.mine が存在する場合はそれを編集する．
シェルに記述したら，以下のコマンドを入力する．
$ source .tcshrc注 2 (注 2).tcshrc.mine を編集した場合は$ source .tcshrc.mine 

シェルが正しく書かれているか確認するために，コマンドプロンプトを立ち上げなおす．
もし，コマンドプロンプトに図 2.2.1 のようなコメントが出たら，シェルの書き方が間違っている．
パスが通っているか確認するには以下のコマンドを入力すればよい．
$ which gnuplot 
このコマンドで gnuplot の場所が出れば，パスは通っている．
1
2

ディレクトリツリーの中のある場所からある場所への経路のこと．
ユーザと OS 間のコマンドのやりとりを仲介する．インタフェイスとしての役割を持つプログラムの総称．

6

図 2.2.1:

2.2.3

MPICH のインストール

1. MPICH のダウンロード
以下のコマンドを使い，MPICH(mpich.tar.gz) を入手する．
$ wget ftp://ftp.mcs.anl.gov/pub/mpi/mpich.tar.gz 

このコマンドを使うと，コマンドプロンプトに図 2.2.2 のようなメッセージが出て，ダウンロードできる．

図 2.2.2:

wget コマンドが使えない場合，インターネットブラウザの URL 欄に
ftp://ftp.mcs.anl.gov/pub/mpi/mpich.tar.gzと入力すればダウンロードできる．

2. MPICH のインストール/usr/local/mpich にインストールする場合を例とする．
mpich.tar.gz を展開する．
$ tar -zxvf mpich.tar.gz









展開後，mpich- version に移動 ($ cd mpich-version  ) し，次のコマンドでインストー

ルを行う．

7





$ ./configure -fc=/opt/∗ ∗ ∗注 3 /ifort -f90=/opt/∗ ∗ ∗注 3 /ifort-cc=/opt/∗ ∗ ∗注 3 /icc -rsh=ssh -prefix=/usr/local/mpich $ make $ su  (スーパーユーザーのホームに移動した場合は，mpich- version まで移動
する．)
# make install # exit (注 3)./configure の***の所には，ifort と icc のフルパス3 を入力する．フルパスが分からない場合は，以下のコマンドで調べる．
$ which ifort $ which icc 

Mac OSX の場合

$ ./configure --prefix=/opt/local CPP=’icc -E’ CC=icc CXX=icpc

FC=ifort F77=ifort 




とする．

3. パスを通す
使用しているシェルの任意の場所にsetenv PATH ${PATH}:/usr/local/mpich/bin




と記述．
ただし，上記は tcsh を使用している場合の例であり，シェルによって書き方は異なる．
自分が使用しているシェルが分からない場合は，以下のコマンドで使用しているシェルが分かる．
$ echo $SHELL 
シェルを emacs で開きたいときは，自分のホームで以下のコマンドを使う．(シェルに tcshを使っている場合)
$ emacs .tcshrc注 4 または$ emacs .tcshrc.mine注 4 (注 4) 自分のホームに.tcshrc.mine が存在する場合はそれを編集する．
シェルに記述したら，以下のコマンドを入力する．
3

ルートディレクトリ (/) を基点として，ファイルの場所を表わす文字列．

8






$ source .tcshrc注 5 (注 5).tcshrc.mine を編集した場合は$ source .tcshrc.mine 

シェルが正しく書かれているか確認するために，コマンドプロンプトを立ち上げなおす．
もし，コマンドプロンプトに図 2.2.1 のようなコメントが出たら，シェルの書き方が間違っている．
パスが通っているか確認するには，以下のコマンドを入力する．
$ which mpicc 

2.2.4

RSA 鍵を使った MPI のプロセス通信設定

MPICH の通信に ssh を使用すると，プログラム実行のたびにパスワードを入力しなければならなくなる．そのため，以下の操作によりそれを省略する．マルチコア4 のマシンを使用する場合，サーバーのマシン自体の RSA 鍵5 も登録する必要がある．
1. 使用するマシンの RSA 鍵を生成する
自分のホームで次のコマンドを入力する．







$ ssh-keygen -t rsa
Generating public/private rsa key pair.
Enter file in which to save the key (/home/***/.ssh/id_rsa):
Enter passphrase (empty for no passphrase): Enter same passphrase again: Your identification has been saved in /home/***/.ssh/id_rsa.
Your public key has been saved in /home/***/.ssh/id_rsa.pub.
The key fingerprint is:**:**:**:**:**:**:**:**:**:**:**:**:**:**:**:** ***@***








~/.ssh に id_rsa.pub が生成される．
2. サーバーとなるマシンに鍵を登録する
1 で生成した id_rsa.pub をサーバーとなるマシンに登録する．


自分のホームに移動 (cd $HOME  ) する．

.ssh の中に authorized_keys が存在するか確認する．($ ls .ssh/  )~/.ssh/authorized_keys が存在しない場合
$ cat .ssh/id_rsa.pub > .ssh/authorized_keys ~/.ssh/authorized_keys が存在する場合
4
5

複数のプロセッサコアを 1 つのパッケージに集約したマイクロプロセッサ．
公開鍵暗号の 1 つ

9









$ cat .ssh/id_rsa.pub >> .ssh/authorized_keys 

同様にして，使用する全てのマシンの RSA 鍵をサーバーとなるマシンに追加登録する．

3. 所有者に対してのみ読み込みと書き込みの権限を与える$ chmod 600 .ssh/authorized_keys 


2.3



RD-AUTO の解凍

1. RD-AUTO(分岐計算を行うソフトウエア) の入手先
長山研究室の〜から RD.tgz をダウンロードする．

2. 展開
適当なディレクトリで RD.tgz を展開する．展開するコマンドは$ tar -zxvf RD.tgz と実行する．このとき，RD.tgz は図 2.3.1 のように展開する．

RD
AS
TE
LIB
Useful̲ﬁles
Tools
DEMO
Exothermic-reaction̲Model
AS
TE
Gray-Scott̲ModelAS
TE

図 2.3.1: 展開図

RD/AS は分岐計算を行う RD-AUTO のディレクトリである．RD/TE は時間発展方程式を計算するディレクトリである．RD/LIB と RD/Useful_Files は後述するスクリプトファイル (datacut ,10

extremum など) のソースファイルが入っている．RD/Tools には mpich などが入っている．RD/DEMO
は本マニュアルで説明する分岐計算などを行うデモプログラムが入っている．

2.3.1

Mac OSX を使って RD-AUTO を使う方法

Mac 上で RD-AUTO を使うためには概ね上記の方法で動かせることができる．しかし，MPI を使って計算することから NotePC 等で固定 IP アドレスを使わないで Local.host として MPI を使う場合は注意が必要である．MPICH のインストールでは Mac Ports を使ってインストールした場合，gfortran を認識するれば問題ないが，gfortran を認識しない場合は，上記の方法で MPICHをインストールする必要があります．

11

第3章

発展方程式

3.1

使い方

3.1.1

時間発展方程式のプログラムに対するフローチャート
タᐃࡍࡿࣇ࢓࢖ࣝ
)XQFWLRQFJQXSORWF

6WDUWタᐃࡍࡿࣇ࢓࢖ࣝ
SDUDPHWHUGDWSDUDPHWHUGDW
LQLWLDOBVHWFࡲࡓࡣ,QLWLDOGDWࣃ࣓࣮ࣛࢱ࡜ึᮇ್ࡢタᐃฟຊࡉࢀࡿࣇ࢓࢖ࣝVROXWLRQBYDOXHGDWࣃ࣓࣮ࣛࢱ࡜ึᮇ್ࡢㄞࡳ㎸ࡳึᮇࢹ࣮ࢱࡢฟຊ

IRU 6WHS 6WHS 0D[VWHS6WHS᫬㛫Ⓨᒎ᪉⛬ᘧࢆゎࡃ

6WHS6WHS&XW ͤ

6WHS! 906 ͤ
6WHS*&9' ͤゎuࡸ᫬⣔ิࢹ࣮ࢱࢆJQXSORW࡛ฟຊฟຊࡉࢀࡿࣇ࢓࢖ࣝVROXWLRQBYDOXHGDW
B6WHSB7LPHB'[B3X69B3Y69GDW㏿ᗘࢹ࣮ࢱࢆJQXSORW࡛ฟຊ㏿ᗘࡢィ⟬ฟຊࡉࢀࡿࣇ࢓࢖ࣝ
B6WHSB7LPHB'LIIB'LVB9HOB$YHOGDWฟຊࡉࢀࡿࣇ࢓࢖ࣝ
ILQDOGDW᭱ᚋࡢࢹ࣮ࢱࢆฟຊ(QG※ 1 StepCut については parameter.dat の StepCut の項目を参照 (→ 16 ページ)．
※ 2 VMS については parameter.dat の Velocity_Measure_Start の項目を参照 (→ 16 ページ)．
※ 3 GCVD については parameter.dat の Gnuplot_CutN_Vel_Dis の項目を参照 (→ 16 ページ)．

12

3.1.2

実行前の設定

ここでは，次の Gray-Scott モデルを例に時間発展方程式を解くための設定を行う．なお，反応
拡散系の変数の個数を変更したい場合は，3.2.2 節を読んでください．



1

2

ut = ∆u + (−au + u v),

ε



(∗1)

t > 0,

0 < x < 2π,v = d∆v + h(v − v) − u2 v,
t

c





ただし，次のようにスケーリングを行って計算するものとする．



4π 2
1

2


ut = L2 ∆u + ε (−au + u v),(∗1)

t > 0,

0 < x < 2π,2vt = d 4π ∆v + h(vc − v) − u2 v,

L2





設定に必要なファイルは次の 5 つである:

1. 0parameter.dat
2. Function.c
3. parameter.dat
4. boundary.dat ：このファイルの境界条件数が反応拡散系の変数の数を決める
5. initial_set.c

1. 0parameter.dat
反応項及び拡散項で使われるパラメーターの設定．




1.0
3.0
1.0
0.07
1.0
0.018
500

;
;
;
;
;
;
;

du ; par[1] 方程式 u の拡散係数を第 1 行に記入
dv ; par[2] 方程式 v の拡散係数を第 2 行に記入
eps ; par[3]a; par[4]vc ; par[5]h; par[6]Lx ; par[7]

3 以降は Gray-Scott モデル (∗1) に現れるパラメータ (パラメータの順番は自由)．拡散係数は最初に入れなければならないので，例えば u, v, w の 3 変数の場合には第 3 行にパラメー
タ dw を入れる．
2. Function.c
反応項 (関数 Reaction ) 及び拡散項 (関数 Diffusion ) の設定．RD-AUTO は反応拡散系の空間離散化に対して差分法と Spectral 法のいずれかを用いており，Spectral 法を選択して数値計算する場合は，y = 2π
L x と変数変換を行った次の方程式を解くように設定する必要が
ある．

13





∂u
4π 2 ∂ 2 u= 2 D 2 + f (u),∂tL∂yt > 0,y ∈ [0, 2π]







#include <stdio.h>…………
void Reaction( TEMP TP, double *par, double *x, double *F ){int i;
int N = TP.N;/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */for( i = 1; i <= N; i++ ){/* ↓ u の反応項 f(u) を記述 ↓ */F[i ] = (-par[4]*x[i]+x[i]*x[i]*x[i+N]) / par[3];/* ↓ v の反応項 g(u) を記述 ↓ */F[i+N] = par[6]*(par[5]-x[i+N]) - x[i]*x[i]*x[i+N];/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */}}void Diffusion( TEMP TP, double *par ){/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - *//* ↓ 変数変換を行う記述 ↓ */par[TP.N_EP+1] = (4.0*M_PI*M_PI*par[1]) / (par[7]*par[7]);par[TP.N_EP+2] = (4.0*M_PI*M_PI*par[2]) / (par[7]*par[7]);/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */}

3. parameter.dat

14







256
; N
128
; N_Wave-1; T_SELECT-1; VEL_SELECT
1.0E-05; dt
6.283185307179586; Lx
100; Gnuplot_CutN_Vel_Dis
0.20; Velocity_Site
200000; Velocity_Measure_Start
3; BC_Select1; EE_Select5000; StepCut
300000; MaxStep
12; Gnuplot_Select2; Initial_Select


;
;
;
;
;
;
;
;
;
;
;
;
;
;
;

空間分割数
切断波数
⃝
3
4
⃝
時間分割幅
⃝
6
⃝
7
⃝
8
⃝
9
10
⃝
11
⃝
12
⃝
最大ステップ数
14
⃝
15
⃝



説明のないものは下記参照．

⃝
3 T_SELECT
T_SELECT の値に応じて，プログラムでは時間発展方程式の与え方が異なる．T_SELECT= 1 のとき，次の方程式を解く：
∂u
=T
∂t

(

)

∂2u
D 2 + f (u) ,∂x

ただし，T は周期 (Initial.dat から読み込む)．T_SELECT , 1 のとき，次の方程式を
解く：

∂u
∂2u
= D 2 + f (u).
∂t∂x⃝
4 VEL_SELECTVEL_SELECT の値に応じて，プログラムでは時間発展方程式の与え方が異なる．VEL_SELECT= 1 のとき，次の方程式を解く：
∂u
∂u
∂2u
=D 2 +c+ f (u),∂t∂x∂xただし，c は速度 (Initial.dat から読み込む)．VEL_SELECT , 1 のとき，次の方程式
を解く：

∂u
∂2u
= D 2 + f (u).
∂t
∂x




T_SELECT と VEL_SELECT の両方が 1 のとき，次の方程式を解く：
∂u
1
=
∂t
T

(

)

∂u
∂2u
D 2 +c+ f (u) .
∂x
∂x





15

⃝
6 Lx
空間 x に対する区間長．ただし，Spectral 法を用いて計算する場合は 2π を指定する．

⃝
7 Gnuplot_CutN_Vel_Dis
Gnuplot_CutN_Vel_Dis = n のとき，n ステップ1 おきに速度データをファイル
82_Step_Time_Diff_Dis_Vel_Avel.dat に出力する．
⃝
8 Velocity_Site
Velocity_Site は解の速度計算に必要とする．解 u の高さの位置を指定する (図 3.1.1)．

1

0.8

v

0.6

0.4

Velocity̲Site
0.2

u
0

0

1

2

3

4

5

6

図 3.1.1: Velocity_Site⃝
9 Velocity_Measure_Start
解の速度を測り始める Step 数を指定する．

10 BC_Select⃝
境界条件の選択

1:Dirichlet 境界条件
2:Neumann 境界条件
3:周期境界条件
11 EE_Select⃝離散スキームの選択

1:差分法 (陰解法)
2:Spectral 法 (Euler 法)
3:Spectral 法 (Runge-Kutta 法)
ただし，Spectral 法の場合は境界条件を周期境界条件として与えること．

12 StepCut⃝
StepCut = n のとき，n ステップおきにデータ (solution_value.dat (24 ページ)，
81_Step_Time_Sup_L1_L2.dat (24 ページ)) を出力する．
1

ステップの定義はフローチャートの step のこと

16

14 Gnuplot_Select⃝
時間発展方程式を解きながら計算途中の解情報を gnuplot を用いて表示する．デフォ
ルトでは 1 ∼ 6 の gnuplot 画面の設定がされている．Gnuplot_Select = 0 と指定したとき，6 つの gnuplot 画面が表示され，Gnuplot_Select = n1 n2 n3 . . . と指定すると
n1 , n2 , n3 , . . . で与えた数字の gnuplot 画面が表示されない．ただし，n1 , n2 , n3 , . . . は
1 ∼ 6 の値である．例えば Gnuplot_Select = 12 のとき 3,4,5,6 の gnuplot 画面が表示される．gnuplot 画面 1 ∼ 6 の意味は次の通りである．

1:初期値の出力
2:u(赤)，v(青) の出力
3:u だけ出力
4:Velocity_Site で指定した値を用いてパルス u のフロントとバックの距離 d を出力(Velocity_Site = 0.2 で与えたときの図 1 を参照)．
5:速度 (赤)，平均速度 (青) の出力
6:∥u∥∞ を時系列データとして出力 (∥u∥L1 ，∥u∥L2 を選択することも可能．)
ただし，u(xi ) = u(i) として∥u∥∞ := max(|u(1)|, . . . , |u(N )|)∥u∥L1 := |u(1)| + · · · + |u(N )|∥u∥L2 :=√

u(1)2 + · · · + u(N )2

gnuplot.c の項目 (→ 25 ページ) も参照のこと．
15 Initial_Select⃝
初期値と使用するパラメータの選択

1:初期値 u0 (x) は initial_set.c (19 ページ) で設定したものを用いる．パラメータは0parameter.dat のものを用いる．
2:初期値 u0 (x) は Initial.dat の 5 列目に入っているものを用いる．パラメータはInitial.dat の 3 列目に入っているものを用いる．
3:初期値 u0 (x) は Initial.dat の 5 列目に入っているものを用いる．パラメータは0parameter.dat のものを用いる．
4. boundary.dat

0.000000000;
0.000000000;
1.000000000;
1.000000000;

BC_000_U
BC_N+1_U
BC_000_V
BC_N+1_V; BC[1]; BC[2]; BC[3]; BC[4]

周期境界条件の場合は値を使うことはないが，変数の個数はここで判定しているために変数の個数 ×2 個の境界条件は設定しておかなければならない．

• Dirichlet 境界条件のとき
u(t, 0) =BC[1]，u(t, L) =BC[2]，v(t, 0) =BC[3]，v(t, L) =BC[4]．
• Neumann 境界条件のとき
ux (t, 0) =BC[1]，ux (t, L) =BC[2]，vx (t, 0) =BC[3]，vx (t, L) =BC[4]．
17• 周期境界条件のとき
u(t, 0) = u(t, L)，v(t, 0) = v(t, L)．

18

5. initial_set.c#include <stdio.h>…………
void initial_set2( TEMP *TP, double *par, double *BC, double *x ){…………
if( IS == 1 ){…………
/*ここに初期値を記述する*/}
else if( IS == 2 || IS == 3 ){…………
}else{
…………
}
…………
}
…………




初期値 u0 (x), v 0 (x) の設定は関数 initial_set2 内で行う．例えば図 3.1.2 のステップ関数は次のようになる．


if( IS == 1 ){for( i = 1; i <= N; i++ ){x[i ] = 0.000000000;
x[i+N] = 1.000000000;
if( N/10 <= i && N/5>= i )x[i] = 0.4;
if( 1 <= i && N/10>= i )x[i+N] = 0.4;
}
}

19



1

v

0.8

0.6

0.4

u

0.2

0

0

1

2

3

4

5

6

図 3.1.2: ステップ関数

3.1.3

実行コマンドの説明

$ ./aout  で時間発展方程式のプログラムを実行する．

aout はコンパイルを行い，コンパイルが成功したらプログラムを実行するまでのスクリプトである．aout は次のコマンドを行う．
$ make $ time␣./a.out 







aout のオプションについて• $ ./aout␣f  のコマンドを実行すると，$ cp␣final.dat␣Initial.dat $ ./aout 

と同等のことを行う．
• $ ./aout␣1  のコマンドを実行すると，オブジェクトファイルを再構成して実行する．す
なわち，





$ rm␣-rf␣ZZZZ_OBJ $ make $time␣./a.out と同等のことを行う．ただし，ZZZZ_OBJ はオブジェクトディレクトリのことで，全てのオブジェクトファイルはこのディレクトリの中にある．
• $ ./aout␣directoryname 
ディレクトリ (directoryname) を作り，その中に0parameter.dat
20

81_Step_Time_Sup_L1_L2.dat
82_Step_Time_Diff_Dis_Vel_Avel.dat
Initial.dat
boundary.dat
final.dat
parameter.dat
solution_value.dat
をコピーする．
• $ ./extremum ∥u∥∞ , ∥u∥L1 ，∥u∥L2 の時系列データ (81_Step_Time_Sup_L1_L2.dat ) の極大値 (もしくは極小値) を求め，周期解の周期を推測するときに用いるスクリプトファイル (→ 27 ページ)．
• $ ./datacut 
solution_value.dat を用いて時間発展方程式の解を描写する (→ 30 ページ)．

21

3.1.4出力されるデータファイル

1. final.dat(Initial.dat)
行数

1
2
3

Step
MaxStep..
.

パラメータ 1

パラメータ 2方向ベクトル

解の値N

par[1]

N_EP
Lx

par[2]
par[3]
..
.

0.0
..
.

x[1] = u[1]
..
.

4
5

dt
NM

6
7
8

PF
EES
NW

9
10
11
..
.

ME
REV_i_E
0.0
..
.

19
20
21
22
23
24
25
26
..
.

0.0
N_VAR
BC[1]
BC[2]
BC[3]
BC[4]
c_MP
0.0
..
.

..
.
par[n]
0.0..
.

Ns
Ns + 1..
.

..
.
x[Ns] = u[Ns]x[Ns+1] = v[1]
..
.

2N s
2N s + 1x[2Ns] = v[Ns]T

2N s + 2
2N s + 3..
.
MaxStep

..
.
0.0

..
.
0.0

• 1 列目最大ステップ数• 2 列目– N
空間分割数– N_EP
パラメータの個数
– Lx
区間
22

..
.
0.0

c
λ

– dt
時間分割幅
– NM
時間発展方程式では使用しない (分岐図の方で用いる)– PF
時間発展方程式では使用しない (分岐図の方で用いる)– EES
離散スキームの選択
1:差分法
2:スペクトル法 (Euler)
3:スペクトル法 (Runge-Kutta)– NW
スペクトル法を用いて計算する場合の切断波数 (差分法では意味のない値)– ME
時間発展方程式では使用しない (分岐図の方で用いる)– REV_i_E
時間発展方程式では使用しない (分岐図の方で用いる)– N_VAR
変数の個数– BC[1] ∼ BC[4]境界の値 (→ 17 ページ)
– c_MP
時間発展方程式では使用しない (分岐図の方で用いる)• 3 列目パラメータ par[1]∼par[N_EP] の値．ただし，N_EP は設定したパラメータの個数．

• 4 列目
分岐計算の場合は方向ベクトルが入る．(発展方程式の場合は 0 の値が入る)• 5 列目
上から順に，– 解の値 (u[1]. . . u[Ns],v[1]. . . v[Ns])– T …周期– c …解の平均速度– λ …コントロールパラメータの値が入る．

23

2. 81_Step_Time_Sup_L1_L2.dat
1 列目
ステップ数
..
.

2 列目
時間
..
.

3 列目
∥u∥∞
..
.

4 列目
∥u∥L1
..
.

5 列目
∥u∥L2
..
.

3. 82_Step_Time_Diff_Dis_Vel_Avel.dat
1 列目
ステップ数
..
.

2 列目
時間
..
.

3 列目
※
..
.

4 列目
あとから
..
.

5 列目
解の速度
..
.

6 列目
解の平均速度
..
.

3 列目の ※ は Velocity_Site (→ 16 ページ) の高さでのパルスのフロントとバックの距離d(図 3.1.3)

1

v

0.8

0.6

0.4

Velocity̲Site
0.2

u

0

フロントとバックの距離d
0

1

2

3

4

5

図 3.1.3: distance

4. solution_value.dat
時間 t

位置 x

解 u の値

解 v の値

t1
..
.
..
.
t1
t2
..
.
..
.
t2
..
.

x1

u[1] の値

v[1] の値x2..
.

u[2] の値..
.

v[2] の値
..
.

xN s
x1

u[Ns] の値
u[1] の値

v[Ns] の値
v[1] の値x2..
.

u[2] の値..
.

v[2] の値
..
.

xN s
..
.

u[Ns] の値..
.

v[Ns] の値
..
.

24

6

3.1.5gnuplot の設定

1. gnuplot.c#include <stdio.h>…………
void Gnuplot_InitialCommand( TEMP *TP, int *c_G ){…………
/*- - - - ↓解の描写範囲を指定する↓- - - - - - - - - - - -*/fprintf( GP2, "set yr[-0.1:1.1]\n" );fprintf( GP3, "set yr[-0.1:1.1]\n" );/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/…………
}void Gnuplot_Draw( TEMP TP, int *c_G ){…………
}void Gnuplot_Close( TEMP *TP ){…………
}void Gnuplot_DrawCommand( TEMP TP, FILE *GP, int SELECT ){…………
switch( SELECT ){…………
case 6:/*↓ u 2:3 だと sup, u 2:4 だと L1, u 2:5 だと L2 の norm を出力する*/
fprintf( GP, "plot ’%s’ u 2:3 pt 7 ps 2.0 lt 1\n", G3 );
…………
}
…………
}





関数 Gnuplot_InitialCommand では gnuplot の初期設定を行う．解の描写範囲もここで指定する．
関数 Gnuplot_Draw では関数 Gnuplot_DrawCommand を呼び出し gnuplot による画面出力を行う．
関数 Gnuplot_Close では gnuplot を閉じる．

25

関数 Gnuplot_DrawCommand では各 gnuplot 画面にどのような出力表示を与えるか指定する．

26

3.1.6周期解の周期の求め方

ここでは 81_Step_Time_Sup_L1_L2.dat の 5 列目のデータ ∥u(T̂ )∥L2 を用いて周期を近似的に求める方法を説明する．具体的な使い方は，第 4 章 2 節の Appendix1.1(→ 96 ページ) で説明する．∥u∥L2 の時系列データが次の図 3.1.4 のように得られたとき，extremum を実行すると図 3.1.5の時系列の極大値 ∥u(T1 )∥L2 , ∥u(T2 )∥L2 , . . . , ∥u(Tn )∥L2 が求められ，Ti+1√ − Ti とその平均の値

n −T1
Tave = Tn−1も求められる．ただし n は極大値の個数．∥u(T, x)∥L2 :=値も求められる．

||u||L∫L2−L u dx．同様に極小

2

2.645

"81_Step_Time_Sup_L1_L2.dat" u 2:5

2.64

2.635

2.63

2.625

2.62

2.615
0

100

200

300

400

500

600

700

800

900

1000

d

図 3.1.4:||u||L2

2.645’81_Step_Time_Sup_L1_L2.dat’ u 2:5’./EXTREMUM_DATA/1_LocalMaxData.dat’ u 2:5’./EXTREMUM_DATA/5_LocalMinData.dat’ u 2:5|| u(T1)||L

2

2.64

2.635

2.63

2.625

^

2.62

|| u(T1)||L

2

^

T1

2.615
0

100

^

T1
200

T2
300

T2
400

500

^

T3

T3

600

700

^

T4
800

T4
900

1000

d

図 3.1.5:$ ./extremum  と実行する．
解の極大値もしくは極小値の間隔と周期，及びその平均値が求められ，ディレクトリ EXTREMUM_DATAにファイルに出力される．出力されるファイルについては後述．ここで求めた周期を使って周期解の分岐計算を行う．
27

• extremum の設定について#!/bin/sh#----------------------------------------------------------------Index_Comparison_Data=5
Index_Period_Data=2#- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - N_Comparison_Max=10
N_Comparison_Min=10#- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - # Gnuplot_Command0="set term dumb"
Gnuplot_Command="u 2:5"#- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - …………
– Index_Comparison_Data
81_Step_Time_Sup_L1_L2.dat の norm データを用いて周期を求める．
例 Index_Comparison_Data=5 の場合，81_Step_Time_Sup_L1_L2.dat の左から 5 列目のデータ (つまり ∥u∥L2 ) を使って周期を求める．

– N_Comparison_Max ，N_Comparison_Min
時系列データの点を順番に y1 , y2 , . . . と表記すると，yi が極大値であるか調べるには
yi−l , . . . , yi−1 < yi かつ yi > yi+1 , . . . , yi+l となればよい．ただし，l は正の定数でl の値を N_Comparison_Max に指定する．同様に，yi が極小値であるか調べるには
yi−r , . . . , yi−1 > yi かつ yi < yi+1 , . . . , yi+r となればよいので，r の値を N_Comparison_Minに指定する．
例 N_Comparison_Min が 10 だったとすると，ある点の 10 個前と 10 個後全てがその点より大きかった場合，その点を極小値とする (図 3.1.6)．

– Gnuplot_Command
Gnuplot_Command は"u n:m"と与える．ただし，n と m は n = Index_Period_Data，m = Index_Comparison_Data である．
• 出力されるファイルについて実行後、ディレクトリ EXTREMUM_DATA が作成される．

– 1_LocalMaxData.dat
すべての極大値の値 (∥u(T1 )∥L2 , ∥u(T2 )∥L2 , . . . , ∥u(Tn )∥L2 )– 2_LocalMaxPeriod.dat
すべての極大値に対する時間間隔 (T2 − T1 , T3 − T2 , . . . , Tn − Tn−1 )– 3_LocalMaxAveData.dat∥u(T1 )∥L2 + ∥u(T2 )∥L2 + · · · + ∥u(Tn )∥L2
極大値の平均値 ()
n28

2.645’81_Step_Time_Sup_L1_L2.dat’ u 2:5’./EXTREMUM_DATA/1_LocalMaxData.dat’ u 2:5’./EXTREMUM_DATA/5_LocalMinData.dat’ u 2:5
2.64

2.635

2.63

2.625

2.62

例.矢印の範囲まで比較して極小値を出している

2.615
0

100

200

300

400

500

600

700

800

900

1000

図 3.1.6:– 4_LocalMaxAvePeriod.dat
極大値の時間間隔に対する平均値 (Tave =Tn − T1)
n−1– 5_LocalMinData.dat
すべての極小値の値 (∥u(T̂1 )∥L2 , ∥u(T̂2 )∥L2 , . . . , ∥u(T̂n )∥L2 )– 6_LocalMinPeriod.dat
すべての極小値に対する時間間隔 (T̂2 − T̂1 , T̂3 − T̂2 , . . . , T̂n − T̂n−1 )– 7_LocalMinAveData.dat∥u(T̂1 )∥L2 + ∥u(T̂2 )∥L2 + · · · + ∥u(T̂n )∥L2
極小値の平均値 ()
n– 8_LocalMinAvePeriod.dat
T̂n − T̂1
極小値の時間間隔に対する平均値 (T̂ave =
)
n−1

29

3.1.7

データカットsolution_value.dat のデータを用いて解 u(x) の描写を行うとき，
$ ./datacut と実行する．
• 設定#!/usr/bin/perl -w#--------------------------------------------------------------------$InitialSelect= 3;$OutPut_NumData = 1;
$sleep= 0.0;$command0= "set yr[-0.1:1.1]";$command11= "";$command12= "u 1:2 w l, u 1:3 w l";$command13= "";$command10= "plot";$DivsionPoint= 1;$DeletingData= "1";#- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - …………


– デフォルトのオプション指定$InitialSelect の値によって描画の表示が変化する．
$InitialSelect = 1 ならば，Return Key を押すことで 1 枚ずつ描写する．
$InitialSelect = 2 ならば，Return Key を押すことで最初から最後まで自動で描写する．描写が終わったら，もう一度 Return Key を押すことで終了する．
$InitialSelect = 3 ならば，実行すると最初から最後まで自動で描写する．描写が終わったら自動的に終了する．
– 縦軸の描写範囲を変えたい場合$command0="set yr[-0.1:1.1]"; の数字を変える．

• 出力されるファイルについて実行後、ディレクトリ cut が作成される．
データは全てその中に入る．

• オプション
– $ ./datacut␣i1 
Return Key を押すことで 1 枚ずつ描写する ($InitialSelect = 1 で実行するのと同じ)．
– $ ./datacut␣i2 
Return Key を押すことで最初から最後まで自動で描写する．描写が終わったら，もう一度 Return Key を押すことで終了する ($InitialSelect = 2 で実行するのと同じ)．
30

– $ ./datacut␣i3 
実行すると最初から最後まで自動で描写する．描写が終わったら自動的に終了する($InitialSelect = 3 で実行するのと同じ)．
– $ ./datacut␣on 
画面を n 回置きに表示する．
– $ ./datacut␣sx 
画面を x 秒置きに表示する ($InitialSelect = 2, 3 のときのみ)．

31

3.2

付録

3.2.1

プログラムで扱う変数について

1. N と N s の定義
N の値は境界条件によって異なる．空間分割数を N̂ と与えたすると，
Dirichlet 境界条件のとき，N = N̂ − 1
Neumann 境界条件のとき，N = N̂ + 1
周期境界条件のとき，N = N̂となる．また，N s = N̂ + 1

2. x と t_x の定義
L
区間 L を間隔 ∆x = N
で離散化して，区間の格子点を xj = (j − 1)∆x, (1, . . . , N + 1) と表す．このとき，定常解 u(x) の格子点上の u(xi ) の値を配列 u[i] に保存し，周期解の場合
は u(0, xi ) の格子上の値を u[i] に保存する．

v についても同様．
配列 x の 1 ∼ N には u[i] の値が，N + 1 ∼ 2N には v[i] の値が入っている，配列 t_x の
1 ∼ N s には u[i] の値が，N s + 1 ∼ 2N s には v[i] の値が入っている:




u[1]
 . 
 . 
 . 



u[N]

x=,v[1]

 . 
 .. 



v[N]

3.2.2





u[1]
 . 
 . 
 . 



u[Ns]
t_x =  v[1]  . 
 .. 


v[Ns]反応拡散系の変数を変更する

RD-AUTO 全体を通じて，反応拡散系の変数を変更するためには，boundary.dat の中の境界条件数を増減することによって実現している．1 変数反応拡散方程式ならば 2 行
0.000000000;
0.000000000;

BC_000_U
BC_N+1_U; BC[1]; BC[2]

2 変数反応拡散系ならば 4 行
0.000000000;
0.000000000;
1.000000000;
1.000000000;

BC_000_U
BC_N+1_U
BC_000_V
BC_N+1_V; BC[1]; BC[2]; BC[3]; BC[4]

3 変数反応拡散系ならば 6 行

32





0.000000000;
0.000000000;
1.000000000;
1.000000000;
2.000000000;
2.000000000;

BC_000_U
BC_N+1_U
BC_000_V
BC_N+1_V
BC_000_W
BC_N+1_W; BC[1]; BC[2]; BC[3]; BC[4]; BC[5]; BC[6]

以下どうように変数を増やすことは可能である．ただし，変数を増やす毎に，0parameter.dat
内におてい，par[1], par[2], par[3],· · · のように順々にパラメータの先頭から拡散係数を決める必要がある．それに応じて，par[TP.N_EP+1], par[TP.N_EP+2], par[TP.N_EP+3], と変数分だけ変数変換も変更する必要がある．区間長のパラメータ番号が変更になったら，もちろん
par[7] を変更するように気をつけてください．
void Diffusion( TEMP TP, double *par ){/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - *//* ↓ 変数変換を行う記述 ↓ */par[TP.N_EP+1] = (4.0*M_PI*M_PI*par[1]) / (par[7]*par[7]);par[TP.N_EP+2] = (4.0*M_PI*M_PI*par[2]) / (par[7]*par[7]);par[TP.N_EP+3] = (4.0*M_PI*M_PI*par[2]) / (par[7]*par[7]);/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
}



33



第4章

分岐計算プログラム

4.1

使い方

4.1.1

実行前の設定1. aout
RD-AUTO は並列計算 (MPICH) を用いるため，最初に並列計算を行うための設定を行う．
MPICH で使用するプロセス数を変更するには，AS/aout を開き，$NumProcs の数を変えることで変更することができる．次の例は，4CPU を持つ cherry09 のマシン環境での設定である．使用するマシンの設定は$HeadName_Machine と@i_Machine を組み合わせて行い，$MachineFileName で指定されたファイルに使用するマシン名が出力される．$NumProcs を1 にすると，1 プロセスでの計算になる．(ホスト名に数字がない場合等の設定は未完成)#!/usr/bin/perl -w#------------------------------------------------------------------------$NumProcs= 4;@i_Machine = (9,9,9,9);#- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - …………
#- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - $MachineFileName= "9_MachineFileName";$HeadName_Machine = ’cherry0’;#------------------------------------------------------------------------…………


この場合 9_MachineFileName に次のようにマシン名が書き込まれ (ただし，() 内は注釈である．)，
cherry09cherry09cherry09cherry09(の 1 個目の CPU)(の 2 個目の CPU)(の 3 個目の CPU)(の 4 個目の CPU)

上から順に$NumProcs の数だけ使用する．
また次の例は，4CPU を持つ peach201〜peach204 を繋げたマシン環境での設定である．

34





#!/usr/bin/perl -w#------------------------------------------------------------------------$NumProcs= 16;@i_Machine = (1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4);#- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - …………
#- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - $MachineFileName= "9_MachineFileName";$HeadName_Machine = ’peach20’;#------------------------------------------------------------------------…………


この場合も 9_MachineFileName に次のようにマシン名が書き込まれる (ただし，() 内は注釈である．)．
peach201
peach202
peach203
peach204
peach201
peach202..
.



(の 1 個目の CPU)(の 1 個目の CPU)(の 1 個目の CPU)(の 1 個目の CPU)(の 2 個目の CPU)(の 2 個目の CPU)

peach203 (の 4 個目の CPU)
peach204 (の 4 個目の CPU)

この場合も上から順に$NumProcs の数だけ使用する．

2. Function.c
反応項 (関数 Reaction ) 及び拡散項 (関数 Diffusion ) の設定を行う．この設定は時間発展方程式の Function.c で設定したものと同じである．

35





#include <stdio.h>…………
void Reaction( TEMP TP, double *par, double *x, double *F ){…………
for( i = 1; i <= N; i++ ){/* ↓ u の反応項 f(u) を記述 ↓ */F[i ] = (-par[4]*x[i]+x[i]*x[i]*x[i+N]) / par[3];/* ↓ v の反応項 g(u) を記述 ↓ */F[i+N] = par[6]*(par[5]-x[i+N]) - x[i]*x[i]*x[i+N];}}
void PartialDiff_Function( TEMP TP, double *par, ………… ){…………
for( i = 1; i <= N; i++ ){…………
}}void i_Diffusion( int *i_par ){…………
}void Diffusion( TEMP TP, double *par ){…………
}void PartialDiff_Function_Di( TEMP TP, double *par, double *parD )
{
…………
}





3. Initial.dat
時間発展方程式によって求めた final.dat や，分岐計算によって求めた final.dat ，middle_final.dat，9middle_final.dat ，固有関数のデータなどを Initial.dat にコピーして使用する．どのような場合にどの初期値を使うのかはデモを参照すること．

4. parameter.dat

36





0.000000
0.0001-130
1.0E-03
1
-1
2
1
3
6
-1
-1
1.0E-08
1
1
0
1

;
;
;
;
;
;
;
;
;
;
;
;
;
;
;
;
;
;



ds
PseudoNewton
KellerEqn
MaxIterdt
N_MSM
PeriodFix
ConPar
SolSelect
BC_Select
TE_Select
EigenValue
EigenVector
MachineEpsion
StepCut
MaxStep
Gnuplot
Initial; パラメータの変化量
; ⃝
2
; ⃝
3
; ⃝
4
; ⃝
5
; ⃝
6
; ⃝
7
; ⃝
8
; ⃝
9
; ⃝
10
; ⃝
11
; ⃝
12
; ⃝
13
; ⃝
14
; ⃝
15
; 最大 step 数
; ⃝
17
; ⃝
18

⃝
2 PseudoNewton
Newton 法の反復 n 回目に対する近似解の norm 値 (CC_Norm ) が PseudoNewton で指定した値 (0.0001) より小さいとき，反復 n 回目は Newton 法から擬似 Newton 法へ変
更する．
⃝
3 KellerEqn
KellerEqn = 1 なら疑似弧長法を用いる．
KellerEqn , 1 なら疑似弧長法を用いない．
⃝
4 MaxIter
MaxIter は Newton 法による最大反復回数
最大反復回数を越えた場合は，解が収束しなかったとしてプログラムを終了する．

⃝
5 dt
時間分割幅 (周期解の分岐計算に必要とする)⃝
6 N_MSM
N_MSM = 1 のとき，single shooting 法を用いる．
N_MSM ≥ 2 のとき，multiple 法を用いる．multiple 法は周期軌道上の n 個の初期値を用意する．
⃝
7 PeriodFix
周期をコントロールパラメータとして扱うかどうかの設定．

PeriodFix = 1 のとき，周期をコントロールパラメータとする．
PeriodFix , 1 のとき，ConPar で指定したパラメータをコントロールパラメータとする．
⃝
8 ConPar37

コントロールパラメータとして扱うパラメータ番号の設定．

ConPar = n のとき，Initial.dat の左から 3 列目，上から n 行目のパラメータをコ
ントロールパラメータとして与える．
⃝
9 SolSelect
求めるべき解を指定する．

SolSelect = 1 :定常解．
SolSelect = 2 :進行解．
SolSelect = 3 :脈動定常解．
SolSelect = 4 :脈動進行解．
10 BC_Select⃝
境界条件の選択．

1:Dirichlet 境界条件．
2:Neumann 境界条件．
3:周期境界条件．
11 TE_Select⃝離散スキームの選択．

1:差分法を用いる．
2:Spectral 法による波数計算を Euler 法で解く．
3:Spectral 法による波数計算を 4 次の Runge-Kutta 法で解く．
4:差分法を用いる．(一周期分の各時刻の解をメモリに保存する．メモリを大量に使用するが，1 より高速に解ける．)．
5:Spectral 法による波数計算を Euler 法で解く (一周期分の各時刻の解をメモリに保存する．メモリを大量に使用するが，2 より高速に解ける．)．
6:Spectral 法による波数計算を 4 次の Runge-Kutta 法で解く (一周期分の各時刻の解をメモリに保存する．メモリを大量に使用するが，3 より高速に解ける．)．
しかし 4，5，6 を選ぶと非常にメモリーを使うので、エラーが出たら 1 か 2 か 3 にする．

12 EigenValue⃝
EigenValue = 1 なら解 (x(∗) ) を求めた後，線形化行列 A(x(∗) ) の固有値計算を行う．
13 EigenVector⃝
EigenVector = 1 なら解 (x(∗) ) を求めた後，線形化行列 A(x(∗) ) の固有値に対する固有関数の計算をする．
14 MachineEpsion⃝
Newton 法の反復 n 回目に対する近似解の norm 値 (CC_Norm ) が MachineEpsilon で指定した値 (10−8 ) より小さいとき，反復 n 回目の近似解を Newton 法により解が収束したと判定する．
15 StepCut⃝
StepCut = n のとき，n ステップおきにデータを出力する．
17 Gnuplot (未完成)⃝
Initial (後述), 9 のとき，Newton 法を解きながら計算途中の解情報を gnuplot を用いて表示する．デフォルトでは 1 ∼ 7 または 1 ∼ 9 の gnuplot 画面の設定がさ38れている．Gnuplot_Select = 0 と指定したとき，n つの gnuplot 画面が表示され，Gnuplot_Select = n1 n2 n3 . . . と指定すると n1 , n2 , n3 , . . . で与えた数字の gnuplot 画
面が表示されない．ただし，n1 , n2 , n3 , . . . は 1 ∼ n の値である．例えば Gnuplot_Select= 12 のとき 3,4,5,6 の gnuplot 画面が表示される．gnuplot 画面 1 ∼ n の意味は次の通りである．
定常解，進行解のとき

1:初期値の出力 (赤が u，青が v)
2:u の解と v の解の出力
3:u の解4:横軸にパラメータ，縦軸に ∥u∥∞ を与えた分岐図
5:横軸にパラメータ，縦軸に固有値の実部の値
6:5 の拡大図
7:横軸にパラメータ，縦軸に ∥u∥L2 を与えた分岐図 (Newton 法の反復途中の解も分岐図中に表示する)
脈動定常解，脈動進行解のとき

1:初期値の出力 (赤が u，青が v)
2:u の解と v の解の出力
3:u の解4:横軸にパラメータ，縦軸に ∥u∥∞ を与えた分岐図
5:横軸にパラメータ，縦軸に固有値の実部の値
6:5 の拡大図
7:横軸にパラメータ，縦軸に ∥u∥L2 を与えた分岐図 (Newton 法の反復途中の解も分岐図中に表示する)
8:横軸に時間，縦軸に ∥u∥∞ の値
9:横軸に時間，縦軸に ∥u∥L2 の値
18 Initial⃝
1:解構造を求める．
2:分岐点における固有ベクトルのデータを用いて分岐点から伸びる解の枝を求める．
7:既に求めたデータの固有値を求める
8:既に求めたデータの固有ベクトルを求める
9:既に求めたデータをディスプレイに表示する．また，midle_final.dat のステップ
数を修正する (→ 69 ページ)．
5. 0_BranchChangeParameter.dat
parameter.dat の Initial = 2 のとき使用する．


0.01000
0.00
0.0
0.00000000
0.000

;
;
;
;
;

eps
period_delta
velocity_delta
par_deltatime; uini の構成に必要; 周期の微小変量; 速度の微小変量; パラメータの微小変量; uini の構成に必要

390_BranchChangeParameter.dat で指定したパラメータを用いて，分岐点から伸びる解の枝を求めるための初期値を構成する．
• 進行解の場合
uini = u0 + εΦR ,cini = δc ,pini = δp ,ただし，u0 は定常解の pitch-fork 分岐点における解 u，ΦR は分岐点における零固有値に対する固有関数の実部，cini は解の速度，pini はコントロールパラメータ，ε, δc , δp はそれぞれ 0_BranchChangeParameter.dat における eps , velocity_delta , par_deltaである.

• 脈動定常解の場合
2π+ δT , pini = δp ,β
ただし，u0 は定常解の Hopf 分岐点における解 u，ΦR , ΦI はそれぞれ分岐点における零固有値に対する固有関数の実部と虚部，Tini は解の周期，pini はコントロールパ
ラメータ，β は定常解の Hopf 分岐点における固有値虚部の値，ε, δT , δp , t はそれぞれ
0_BranchChangeParameter.dat における eps , period_delta , par_delta , time である. β:定常解の Hopf 分岐点における固有値虚部の値
uini = u0 + ε(sin(2πt)ΦR + cos(2πt)ΦI ), Tini =• 脈動進行解の場合
2π+δT , cini = c0 +δc , pini = δp ,β
ただし，u0 は進行解の Hopf 分岐点における解 u，ΦR , ΦI はそれぞれ分岐点における零固有値に対する固有関数の実部と虚部，cini は解の速度，Tini は解の周期，pini はコン
トロールパラメータ，c0 は進行解の Hopf 分岐点における速度，ε, δT , δc , δp , t はそれぞ
れ 0_BranchChangeParameter.dat における eps , period_delta , velocity_delta, par_delta , time である.
uini = u0 +ε(sin(2πt)ΦR +cos(2πt)ΦI ),40

Tini =

4.1.2

実行コマンドの説明

$ ./aout  で分岐計算のプログラムを実行する．



aout はコンパイルを行い，コンパイルが成功したらプログラムを実行するまでのスクリプトである．aout は次のコマンドを行う．
$ make $ time␣mpirun␣-np␣$NumProcs␣-machinefile␣$MachineFileName␣./a.out 
ただし，$NumProcs は使用するプロセス数，$MachineFileName は使用するマシン名の記述されたファイル (→ 34 ページ参照)．
aout のオプションについて• $ ./aout␣f $ cp␣final.dat␣Initial.dat $ ./aout と同等のことを行う．
• $ ./aout␣1  のコマンドを実行すると，オブジェクトファイルを再構成して実行する．
すなわち，





$ rm␣-rf␣オブジェクトファイル $ make $ ./aout と同等のことを行う．
• $ ./aout␣directoryname 
ディレクトリ (directoryname) を作り，その中に0_BranchChangeParameter.dat
81_CP_TV_NR.dat
82_S_CP_TV_NR_REV_IMV_ABV_NUEV.dat
83_S_CP_TV_NR_REV_IMV_ABV_NUEV.dat
9midle_final.dat
Initial.dat
midle_final.dat
multiple_points.dat
parameter.dat
をコピーする．(Initial.dat は final.dat と中身は同じ)

41

4.1.3解の枝の安定性を求める方法

ある解の枝を求めたとして，その安定性を調べたいとする．このとき以下の作業を行う．
ただし，middle_final.dat には解の枝の情報が入っている．

1. middle_final.dat を Initial.dat に変更する．
cp␣middle_final.dat␣Initial.dat 
2. parameter.dat の Initial の値を 8 にする．
3. 実行
$ ./aout 42

4.1.4出力されるデータファイル

1. middle_final.dat
6 列目，7 列目は場合に応じて出力される．

行数

パラメータ

方向ベクトル

解の値

固有値実部

固有値虚部

N

par[1]

N_EP
Lx

par[2]
par[3]
..
.

u̇[1]
..
.

x[1]= u[1]..
.

ReΦ[1]..
.

ImΦ[1]..
.

Ns
Ns + 1..
.

..
.
u̇[Ns]
v̇[1]
..
.

..
.
x[Ns]= u[Ns]x[Ns+1]= v[1]
..
.

2N s
2N s + 1v̇[Ns]
Ṫx[2Ns]= v[Ns]T..
.
ReΦ[2Ns]
0.0..
.
ImΦ[2Ns]
0.0

ċ
λ̇
u̇(1)
..
.

c
λ
x[1]
..
.

0.0
0.0
ReΦ[1]
..
.

0.0
0.0
ImΦ[1]..
.

1
2
3

Step 数
1
..
.

4
5

dt
NM

6
7
8

PF
EES
NW

9
10
11
..
.

ME
REV_i_E
0.0
..
.

19
20
21
22
23
24
25
26
..
.

0.0
N_VAR
BC[1]
BC[2]
BC[3]
BC[4]
c_MP
0.0
..
.

2N s + 2
2N s + 3
2N s + 4
2N s + 5

..
.
1
2
..
.

..
.
par[n]
0.0
..
.

..
.
0.0
N

..
.
0.0
par[1]

N_EP

par[2]• 1 列目
Step 数• 2 列目
– NM

43

– PF
parameter.dat (→ 36 ページ) の PeriodFix の値．
– ME
parameter.dat (→ 36 ページ) の MachineEpsion の値．
– REV_i_E
39 ページの β の値．
– c_MP
特に気にする必要はない (削除される予定)．
時間発展方程式の final.dat (→ 22 ページ) の 2 列目と同様．

• 3 列目時間発展方程式の final.dat の 3 列目と同様．

• 4 列目方向ベクトル• 5 列目時間発展方程式の final.dat の 5 列目と同様．

• 6 列目
固有値実部• 7 列目
固有値虚部

2. 9middle_final.dat
middle_final.dat の第 1 列から第 5 列目までの情報．
9middle_final.dat のデータが存在する場合，9middle_final.dat の必要なステップ部
分以外のデータを全て消して Initial.dat に変えて，初期値として使うことができる．


viを用いて 9middle_final.dat を編集する方法．
 
d G のコマンドでカーソル位置からファイルの末尾までを削除する．
d 1 G のコマンドでカーソル位置からファイルの先頭までを削除する．



3. final.dat
9middle_final.dat の最後のステップのデータ．
4. 81_CP_TV_NR.dat
1 列目
コントロールパラメータ
..
.

2 列目
解の周期
..
.

5. 82_S_CP_TV_NR_REV_IMV_ABV_NUEV.dat
固有値計算をしているときのみ出力される．

44

3 列目
解の速度
..
.

4 列目
∥u∥∞
..
.

5 列目
∥u∥L1
..
.

6 列目
∥u∥L2
..
.



1 列目
ステップ番号..
.

2 列目
コントロールパラメータ
..
.

8 列目
固有値の実部..
.

9 列目
固有値の虚部
..
.

3 列目
解の周期
..
.

10 列目
固有値の絶対値
..
.

4 列目
解の速度
..
.

2 列目
位置 x
..
.

45

6 列目
∥u∥L1
..
.

11 列目
不安定な固有値の数 (未完成)..
.

6. solution_value.dat
1 列目
時間
..
.

5 列目
∥u∥∞
..
.

3 列目
u の値
..
.

4 列目
v の値
..
.

7 列目
∥u∥L2
..
.

第5章

デモ 1

この章では，RD-AUTO を用いて１次元 Gray-Scott モデルの分岐図を描く工程を説明する．分岐点の分岐の種類についても調べる．
このデモに必要なデータは DEMO ディレクトリの中に図 5.0.1 のように入っている．
DEMO

Gray-Scott_ModelAS

01_INI

02_SPE

͐

TE

01_SP

02_TP

͐

図 5.0.1: DEMO の中身

DEMO 以下には計算に必要な設定ファイルだけがある．実行する場合は，TE（発展方程式計算）あるいは AS（分岐構造計算）にコピーして使用すること1 ．

1

問題ごとに TE 内あるいは AS 内の全てのファイル46

デモプログラム

内容

実行後の保存先

01_SP
02_TP
03_TB1
04_TB2
01_INI
02_SPE
03_SP1
04_SP2
05_SP306_TP1
07_TP208_TP3
09_TB1
10_AP1
11_AP2
12_SB1
13_SB2
14_SB3定常解 (TE)SP
TP

進行解 (TE)
脈動進行解 (TE)
脈動進行解 (TE)
定常解 (AS)

TB
INI

定常解 (AS)
定常解 (AS)
定常解 (AS)
定常解 (AS)
進行解 (AS)
定常解 (AS)
定常解 (AS)

SP1
SP2
SP
TP_INI
TP
TP

脈動進行解 (AS)
脈動進行解 (AS)
脈動進行解 (AS)
脈動定常解 (AS)
脈動定常解 (AS)
脈動定常解 (AS)

SB_INI
SB

47

5.1

1 次元 Gray-Scott モデルの時間発展方程式

ここでは，1 次元 Gray-Scott モデルに現れるパルス波を求める．カレントディレクトリは時間
発展方程式を解くところに移動 ($ cd␣RD/TE  ) して，最初に次の方程式の設定を行う．




































1
4π 2
uxx + (−au + u2 v),
2
L
ε
4π 2
vt = dv 2 vxx + h(vc − v) − u2 v,
L
ut = duu(t, 0) = u(t, 2π),v(t, 0) = v(t, 2π),u(0, x) = u0 (x),v(0, x) = v0 (x),t > 0,

ux (t, 0) = ux (t, 2π), t > 0,vx (t, 0) = vx (t, 2π), t > 0,
0 ≤ x ≤ 2π,
0 ≤ x ≤ 2π
0 < x < 2π, 






(∗1)

}

(∗2)

}

(∗3)





1 方程式 (∗1) の設定⃝
方程式の拡散項と反応項の設定は Function.c で行う．拡散項は関数 Diffusion をvoid Diffusion( TEMP TP, double *par ){/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */par[TP.N_EP+1] = (4.0*M_PI*M_PI*par[1]) / (par[7]*par[7]);par[TP.N_EP+2] = (4.0*M_PI*M_PI*par[2]) / (par[7]*par[7]);/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */}と記述して，反応項は関数 Reaction をvoid Reaction( TEMP TP, double *par, double *x, double *F ){/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */int i;
int N = TP.N;
for( i = 1; i <= N; i++ ){F[i ] = (-par[4]*x[i]+x[i]*x[i]*x[i+N]) / par[3];F[i+N] = par[6]*(par[5]-x[i+N]) - x[i]*x[i]*x[i+N];}/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */}と記述する．

⃝
2 境界条件 (∗2) の設定
境界条件の設定は parameter.dat を開き，BC_Select を

48







3

;

…………
BC_Select
…………

;

OP[9]





と与える．BC_Select = 3 は周期境界条件を表す．

⃝
3 初期値 (∗3) の設定
時間発展方程式で用いる初期値は求めるパルス波に依存して異なるので，各パルス波を求めるところで説明する．
⃝
4⃝
1 〜⃝
3 以外の設定
パルス波を求めるには ⃝1 〜⃝
3 以外にも設定する箇所 (parameter.dat や 0parameter.datなど) があるが，これらの設定は求めるパルス波に対して異なるので，各パルス波を求める
ところで説明する．
Step1 定常パルス波を時間発展方程式を用いて求める
ステップ関数を初期値として進行パルス波を求める．
1. 実行前の設定(a) 0parameter.dat
Function.c 内のパラメータの値を設定する．今回は以下のように設定する．


1.0
2.0
1.1
0.07
1.0
0.018
500

;
;
;
;
;
;
;



du ; par[1]dv ; par[2]eps ; par[3]a; par[4]vc ; par[5]h; par[6]Lx ; par[7](b) Function.c
反応項はf (u, v) = −au + u2 v
g(u, v) = h(vc − v) − u2 vであるので、関数 Reaction 内でfor( i = 1 ; i <= N ; i++ ){F[i ] = ( -par[4]*x[i] + x[i]*x[i]*x[i+N] ) / par[3];F[i+N] = par[6]*( par[5] - x[i+N] ) - x[i]*x[i]*x[i+N];
}



と記述．

(c) parameter.dat
今回は以下のように設定する．

49









256
; N
; 空間分割数は 256
128; N_Wave; 切断波数は 128-1; T_SELECT; (∗1); (∗1)-1; VEL_SELECT
1.0E-01; dt; 時間の刻み幅は 10−1
6.283185307179586; Lx; 区間は 2π100
; Gnuplot_CutN_Vel_Dis; (∗2)
0.2; Velocity_Site; (∗3)
10000; Velocity_Measure_Start ; step10000 から速度を測る
3; BC_Select; 周期境界条件3; redEE_Select; (∗4)
500; StepCut; 500 回おきにデータを出力
100000; MaxStep; 最大 step 数0; Gnuplot_Select; 全ての gnuplot 画面を表示; 初期値はステップ関数
1; Initial_Select(∗1) 時間発展方程式は ut = uxx + f (u), v t = v xx + g(v) を解く(∗2) 速度データを 100 回おきに出力する(∗3) u のパルス波の高さ 0.2 のところで速度を測る(∗4)Spectral 法 (Runge-Kutta 法)(d) boundary.dat
今回は周期境界条件のため，BC[1] ∼ BC[4] は特に設定する必要はない．
(e) initial_set.c
初期値の設定を行う．⃝
3 初期設定 (∗3) の初期値は図 5.1.1 のステップ関数を与える．

1

v

0.8

0.6

0.4

u

0.2

0

0

1

2

3

4

5

6

図 5.1.1:
プログラムでは次のように設定する2 ．
2プログラムで x[i] = 0 , x[i+N] = 1 と与えたのは，Gray-Scott モデルの拡散項を除いた常微分方程式の平衡点がu = 0, v = 1 だからである．

50





void initial_set2( TEMP *TP, double *par, ………… ){…………
if( IS == 1 ){for( i = 1; i <= N; i++ ){x[i ] = 0.000000000;
x[i+N] = 1.000000000;
if( N/10 <= i && N/5>= i )x[i] = 0.4;
if( 1 <= i && N/10>= i )x[i+N] = 0.4;}}
else if( IS == 2 || IS == 3 ){
…………
}
…………
}



(f) gnuplot.c
gnuplot で描写する領域を指定する．図 5.1.1 の縦軸の範囲を −0.1 ∼ 1.1 で設定する．
void Gnuplot_InitialCommand( TEMP *TP, int *c_G ){…………
/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/fprintf( GP2, "set yr[-0.1:1.1]\n" );fprintf( GP3, "set yr[-0.1:1.1]\n" );/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
…………
}



2. 実行



$ ./aout␣1 で実行する．

51











デモプログラム 1
Step1 の設定から実行までの操作は
$ cp␣../DEMO/Gray-Scott_Model/TE/01_SP/*␣./ $ ./aout␣1 
で行える．

RD/DEMO/Gray-Scott_Model/TE/01_SP には Step1 の設定で必要なファイル．
• 0parameter.dat• Function.c• parameter.dat• initial_set.c• gnuplot.cが存在する．


3. 実行結果
定常パルス波が出れば成功．
1

v

0.8

0.6

0.4

0.2

u

0

0

1

2

3

4

5

6

図 5.1.2: 定常パルス波

4. データの保存
ここで求めた定常パルス波のデータを SP というディレクトリに保存する．次のコマンドを実行することでデータをディレクトリ SP に保存できる．

$ ./aout␣SP




52





一般に，データをディレクトリ directoryname に保存する場合は次のコマンドを実行する．
$ ./aout␣directoryname ただし，directoryname には f 及び数字から始まる文字列を指定することはでき
ない．




Step2 進行パルス波
ステップ関数を初期値として，進行パルス波を求める．
1. 実行前の設定(a) 0parameter.dat




…………
; dv ; par[2]
…………

3.0




と変更する．

(b) parameter.dat
変更点なし．
(c) initial_set.c
変更点なし．
2. 実行

$ ./aout␣1  で実行する．


デモプログラム 2
Step2 の設定から実行までの操作は$ cp␣../DEMO/Gray-Scott_Model/TE/02_TP/*␣./





$ ./aout␣1 
で行える．

RD/DEMO/Gray-Scott_Model/TE/02_TP には Step2 の設定で必要なファイル．
• 0parameter.dat• Function.c• parameter.dat• initial_set.c• gnuplot.cが存在する．


3. 結果
進行パルス波が出れば成功．

4. データの保存$ ./aout␣TP 
53

1

v

0.8

0.6

0.4

0.2

u

0

0

1

2

3

4

5

6

図 5.1.3: 進行パルス波
でディレクトリ TP にデータを保存する．

Step3.1 脈動進行パルス波
Step2 の結果を初期値として，脈動進行パルス波を求める．
1. 実行前の設定(a) 0parameter.dat


3.38



…………
; dv ; par[2]
…………





と変更する．

(b) parameter.dat




…………
3
; Initial_Select; (∗1)…………
300000; MaxStep; 最大 step 数…………
(∗1) 初期値 u0 (x) は Initial.dat に入っているものを用いる．パラメータは0parameter.dat のものを用いる．
と変更する．脈動進行解を安定させるために最大 step 数を大きくとる．

(c) Initial.dat
Step2 で出力された final.dat を Initial.dat として使用する．
2. 実行$ cp␣TP/final.dat␣Initial.dat





$ ./aout␣1 で実行する．
54

デモプログラム 3.1
Step3.1 の設定から実行までの操作は
$ cp␣../DEMO/Gray-Scott_Model/TE/03_TB1/*␣./ $ ./aout␣1 
で行える．

RD/DEMO/Gray-Scott_Model/TE/03_TB1 には Step3.1 の設定で必要なファ
イル．
• 0parameter.dat• Function.c• parameter.dat• Initial.dat• gnuplot.cが存在する．


3. 結果
脈動進行パルス波が出れば成功．
1

1

0.8

v

0.8

v

0.6

0.6

0.4

0.4

0.2

0.2

u

u
0

0

0

1

2

3

4

5

6

0

図 5.1.4: 脈動進行パルス波 1

1

2

3

4

5

6

図 5.1.5: 脈動進行パルス波 2

Step3.2 脈動進行パルス波
安定した脈動解を 30000step ほど計算する．
1. 実行前の設定(a) 0parameter.datと変更点なし．
(b) parameter.dat


30000



;

…………
MaxStep
…………



;


と変更する．

(c) Initial.dat
55

Step3.1 で出力された final.dat を Initial.dat として使用する．
2. 実行$ cp␣final.dat␣Initial.dat





$ ./aout␣1 で実行する．


デモプログラム 3.2
Step3.2 の設定から実行までの操作は
$ cp␣../DEMO/Gray-Scott_Model/TE/04_TB2/*␣./ $ ./aout␣1 
で行える．

RD/DEMO/Gray-Scott_Model/TE/04_TB2 には Step3.2 の設定で必要なファ
イル．
• 0parameter.dat• Function.c• parameter.dat• Initial.dat• gnuplot.cが存在する．


3. データの保存

$ ./aout␣TB 
でディレクトリ TB にデータを保存する．

56

5.2

1 次元 Gray-Scott モデルの分岐計算

ここでは，5.1 節で求めたパルス波の分岐図を求める．ただし，d を分岐パラメータとして与える．

カレントディレクトリは分岐図を求めるところに移動 ($ cd␣RD/AS  ) して，最初に Gray-Scott
モデル (48 ページの (∗1)) の拡散項と反応項の設定を行う．これらの設定は時間発展方程式と同じであり，拡散項は Function.c の関数 Diffusion をvoid Diffusion( TEMP TP, double *par ){/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */par[TP.N_EP+1] = (4.0*M_PI*M_PI*par[1]) / (par[7]*par[7]);par[TP.N_EP+2] = (4.0*M_PI*M_PI*par[2]) / (par[7]*par[7]);/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */}と記述して，反応項は Function.c の関数 Reaction をvoid Reaction( TEMP TP, double *par, double *x, double *F ){/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */int i;
int N = TP.N;
for( i = 1; i <= N; i++ ){F[i ] = (-par[4]*x[i]+x[i]*x[i]*x[i+N]) / par[3];F[i+N] = par[6]*(par[5]-x[i+N]) - x[i]*x[i]*x[i+N];}/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */}と記述する．次に周期境界条件 (48 ページの (∗2)) の設定を行う．境界条件の設定は parameter.dat
を開き，BC_Select を




3

;

…………
BC_Select
…………

;

OP[9]





と与える．ただし，BC_Select = 3 は周期境界条件を表す．
Arrow Stay の分岐計算ソフトウェアは Newton 法を用いて，解構造を数値的に求める．解構造の計算は図 5.2.1(a) のように，各 step ごとにパラメータを固定しながら ũ1 → ũ2 → . . . と順次に求める方法や図 (b) のように，パラメータを動かしながら û1 → û2 → . . . と順次に求める方法が
ある．

57

u^

u~

u^

2

3

2

ᗙᶆ

u~

u^

1

1

ΔS
S

S

S

S

図 5.2.1: (a)図 5.2.2: (b)(a)，(b) の選択およびパラメータの変化量 ∆s は parameter.dat 内で設定する．このとき，Newton 法で用いる初期値は求めるパルス波に依存して変わるため，各パルス波の分岐図を求める
ところで説明する．分岐計算を行うには，他にも parameter.dat 内のパラメータ調整や必要に応じて 0_BranchChangeParameter の設定を必要とするが，これらの設定はデモプログラム内で説明する．
Step1.1 定常解の分岐計算
5.1 節で求めた Step1 の解を用いて，定常解の分岐図を求める．ただし，Step1.1 では定常解の安定性を計算しないとする．定常解の選択およびその安定性計算による設定は parameter.dat
を開き，次のように設定する．




1

;

-1
-1

;
;

…………
SolSelect…………
EigenValue
EigenVector
…………

;

OP[8]

;
;

OP[11]
OP[12]



ただし，SolSelect = 1 は定常解を選択し，EigenValue と EigenVector を 1 以外の値に選択すると解の安定性を計算しない．

1. 実行前の設定(a) Initial.dat
5.1 節の Step1 で求めた解を Newton 法の初期値に与える．
$ cp␣../TE/SP/final.dat␣Initial.dat (b) parameter.dat
今回は以下のように設定する．

58





0.000000; ds; パラメータの変化量なし
0.0001; PseudoNewton; (∗1)-1; KellerEqn; 擬似弧長法の設定なし (図 4.2.1(a)); Newton 法の最大反復回数30; MaxIter
1.0E-03; dt; 定常解の分岐計算では不必要な設定
1; N_MSM; 定常解の分岐計算では不必要な設定-1; PeriodFix; 定常解の分岐計算では不必要な設定
2; ConPar; d を分岐パラメータに選択1; SolSelect; 定常解を選択
3; BC_Select; 周期境界条件を選択
6; TE_Select; 定常解の分岐計算では不必要な設定-1; EigenValue; 固有値の計算をしない-1; EigenVector; 固有ベクトルの計算をしない
1.0E-08; MachineEpsion ; (∗2)
; 毎 step ごとにデータを出力する
1; StepCut
1; MaxStep; 最大 step 回数0; Gnuplot; 全ての gnuplot 画面を出力する
1; Initial; 分岐図を求めることを選択(∗1)Newton 法 の 反 復 n 回 目 に 対 す る 近 似 解 の norm 値 (CC_Norm ) がPseudoNewton で指定した値 (0.0001) より小さいとき，反復 n 回目は Newton 法から擬似 Newton 法へ変更する．
(∗2)Newton 法 の 反 復 n 回 目 に 対 す る 近 似 解 の norm 値 (CC_Norm ) がMachineEps で指定した値 (10−8 ) より小さいとき，反復 n 回目の近似解をNewton 法により解が収束したと判定する．


分岐計算のプログラムを実行すると，標準出力に CC_Norm の値などが出力される．
例えば，実行後次の画面が出力されたとすると，

59




Newton 法の反復回数

step 数
最大 step 回数

コントロールパラメータ d の値
収束判定 norm 値
iter 回数目の norm 値
CC̲Norm > PseudoNewton(0.001)
であるため通常のNewton法を用いる

CC̲Norm < PseudoNewton(0.001)
であるため擬似Newton法を用いる

CC̲Norm < MachineEpsilon(0.00000001)であるためNewton法で解が収束したと判定する



反復 1 回目で CC_Norm < PseudoNewton (0.0001) であるため，Newton 法から擬似 Newton 法になる．
反復 2 回目で CC_Norm < MachineEps (10−8 ) であるため，解が収束したと判定
する．

2. 実行
上記の設定が終わったら，$ ./aout␣1 で実行する．


デモプログラム 1.1
Step1.1 の設定から実行までの操作は
$ cp␣../DEMO/Gray-Scott_Model/AS/01_INI/*␣./ $ ./aout␣1 
で行える．

RD/DEMO/Gray-Scott_Model/AS/01_INI には Step1.1 の設定で必要なファイル• Function.c• Initial.dat• parameter.datが存在する．


3. 実行結果
以下のように「Convergence !!」と表示されれば，Newton 法を用いて解が収束したと60

判定する．

step 数

Newton 法の反復回数

最大 step 回数

コントロールパラメータ d の値
収束判定 norm 値
iter 回数目の norm 値

4. データの保存
Step1.1 で求めた定常解のデータを INI というディレクトリに保存する．
次のコマンドを実行することでデータをディレクトリ INI に保存できる．

$ ./aout␣INI







時間発展方程式のデモプログラムで説明したように，一般には次のコマンドを実行することでデータをディレクトリ directoryname に保存できる．
$ ./aout␣directoryname ただし，directoryname には f 及び数字から始まる文字列を指定することはでき
ない．




Step1.2 定常解の分岐計算
次に，Step1.1 で求めた解を用いてパラメータを変化させながら定常解の枝を求める．
1. 実行前の設定(a) Initial.dat
初期値は Step1.1 で求めた解を用いる．
(b) parameter.dat




0.005

;

10
1000

;
;

; パラメータの変化量 0.005

ds…………
StepCut
MaxStep…………

; 10step ごとにデータを出力する; 最大 step 回数と変更する．

61

2. 実行$ cp␣INI/9middle_final.dat␣Initial.dat





$ ./aout␣1 

デモプログラム 1.2
Step1.2 の設定から実行までの操作は
$ cp␣../DEMO/Gray-Scott_Model/AS/02_SPE/*␣./ $ ./aout␣1 
で行える．

RD/DEMO/Gray-Scott_Model/AS/02_SPE には Step1.2 の設定で必要なファイル• Function.c• Initial.dat• parameter.datが存在する．


3. 実行結果
…………
Step / MaxStep
iter / MaxIter
Control Parameter
MachineEps
CC_Norm



=
=
=
=
=

73 / 1000
29 / 30
2.364999999999992
0.000000010000000
0.000006462960947

Step / MaxStep=
73 / 1000
iter / MaxIter=
30 / 30
Control Parameter=
2.364999999999992
MachineEps=
0.000000010000000
CC_Norm=
0.000006461318227--------------------------------------------------------------Not Convergence !!
-------------------------------------------------------------------------------------------------------------------------------Iteration Over !! (30)
Iteration Over !! (30)
Iteration Over !! (30)
Iteration Over !! (30)
Iteration Over !! (30)-----------------------------------------------------------------…………




62

73step において，Newton 法による反復回数が最大反復回数 (この場合は 30 回) を越えてしまい，73step の定常解が収束しなかったと判定される．
Newton 法により解が収束しない原因として，• saddle-node 分岐により解が存在しないところで計算した• パラメータの変化量を大きくしすぎたなどが考えられる．
ここでは，73step 近傍で saddle-node 分岐が起きるため，解が収束しない (図 5.2.3 参照)．図5.2.3 の横軸はコントロールパラメータ d，縦軸は ∥u∥∞ である．

||u||∞
0.235

10step20step0.23"81_CP_TV_NR.dat" u 1:4

30step

求まった解

40step50step60step70step

0.225

0.22

0.215

求められなかった解

0.21

0.205

0.2
2

2.05

2.1

2.15

2.2

2.25

2.3

2.35

2.4

2.45

2.5

d

図 5.2.3: 定常パルス波の分岐図

Step1.3 定常解の分岐計算
Step1.2 の 73step 以降を求めるには擬似弧長法を用いる．ここでは，擬似弧長法を用いてStep1.1 の結果から再び定常解の分岐計算を行う．
1. 実行前の設定(a) Initial.dat
初期値は step1.1 で求めた 9middle_final.dat を用いる．
(b) parameter.dat-0.005;1;; パラメータの変化量 −0.005

ds…………
KellerEqn…………

; 擬似弧長法の設定あり

と変更する3 ．
3ds は負の値に変更したが，擬似弧長法を用いると自分の進みたい方向と逆方向に分岐の解の枝が進むことがおきる．この場合，例えば ds の正負の値を変更して自分の進みたい方向に解の枝が進むよう調整する．

63

2. 実行$ cp␣INI/9middle_final.dat␣Initial.dat $ ./aout␣1 

デモプログラム 1.3
Step1.3 の設定から実行までの操作は
$ cp␣../DEMO/Gray-Scott_Model/AS/03_SP1/*␣./ $ ./aout␣1 
で行える．

RD/DEMO/Gray-Scott_Model/AS/03_SP1 には Step1.3 の設定で必要なファイル• Function.c• Initial.dat• parameter.datが存在する．


3. 実行結果


…………
Step / MaxStep
iter / MaxIter
Control Parameter
MachineEps
CC_Norm=====

1000 / 1000
0 / 304.087204836309247
0.000000010000000
0.005000000000000

Step / MaxStep
iter / MaxIter
Control Parameter
MachineEps
CC_Norm=====

1000 / 1000
1 / 304.092168559232388
0.000000010000000
0.000000011138078

Step / MaxStep=
1000 / 1000
iter / MaxIter=
2 / 30Control Parameter=4.092168623444494
MachineEps=
0.000000010000000
CC_Norm=
0.000000000011126--------------------------------------------------------------Convergence !!
--------------------------------------------------------------…………


実行結果の分岐図を見るには，gnuplot を開き$ plot "81_CP_TV_NR.dat" u 1:4
64

と実行する4 ．分岐図は図 5.2.4 となり，saddle-node 分岐点が 2 つ存在する．点 S は
step1.3 で求めた最初の定常解であり，点 S から点 E1 まで矢印の方向に定常解の枝が求まる．

||u||∞
0.24"81_CP_TV_NR.dat" u 1:4

S

0.23
0.22
0.21
0.2
0.19
0.18
0.17
0.16
0.15

E1

0.14
0.13
1.5

2

2.5

3

3.5

4

4.5

d

図 5.2.4: 定常パルス波の分岐図

4. データの保存
step1.3 で求めた分岐図のデータを SP1 のディレクトリに保存する．
$ ./aout␣SP1 
Step1.4 定常解の分岐計算
次に，図 5.2.4 の点 S から左側の定常解の分岐計算をする．
1. 実行前の設定(a) Initial.dat
点 S を初期値として設定するために，SP1 の 9middle_final.dat を Initial.dat
にコピーし ($ cp␣SP1/9middle_final.dat␣Initial.dat  )， Initial.datの最初のステップ以外のデータを全て削除する．この削除方法に vi を用いる. 以下に vi を用いてファイルを編集する方法を説明する.

4
81_CP_TV_NR.dat は 1 列目にコントロールパラメータ，2 列目に解の周期，3 列目に解の速度，4 列目に ∥u∥∞ ，
5 列目に ∥u∥L1 ，6 列目に ∥u∥L2 が入っている．gnuplot で u 1:4 と指定すると，横軸にコントロールパラメータ，縦軸に ∥u∥∞ のグラフが表示する．

65





i. vi で Initial.dat を開く．
$ vi␣Initial.dat ii. 削除したい部分の最初の行にカーソルを移動する．
$ /␣␣␣␣2  と入力すると 2step 目のデータの最初の行にカーソルが移動する．
iii. カーソル位置からファイルの末尾までを削除する．

 
d G のコマンドでカーソル位置からファイルの末尾までを削除することができる．

iv. ファイルを保存して終了する．
$ :wq  と入力することでデータを保存して終了する．

(b) parameter.dat


0.005




;

; パラメータの変化量 0.005

ds
…………





と変更する．

2. 実行

$ ./aout␣1 で実行する．


デモプログラム 1.4
Step1.4 の設定から実行までの操作は
$ cp␣../DEMO/Gray-Scott_Model/AS/04_SP2/*␣./ $ ./aout␣1 
で行える．

RD/DEMO/Gray-Scott_Model/AS/04_SP2 には Step1.4 の設定で必要なファ
イル．
• Function.c• Initial.dat• parameter.datが存在する．


3. 実行結果

66





…………
Step / MaxStep
iter / MaxIter
Control Parameter
MachineEps
CC_Norm=====

1000 / 1000
0 / 304.258176296652051
0.000000010000000
0.005000000000000

Step / MaxStep
iter / MaxIter
Control Parameter
MachineEps
CC_Norm=====

1000 / 1000
1 / 304.263162309425487
0.000000010000000
0.000000011191342

Step / MaxStep=
1000 / 1000
iter / MaxIter=
2 / 30Control Parameter=
4.263162334360556
MachineEps=
0.000000010000000
CC_Norm=
0.000000000011498--------------------------------------------------------------Convergence !!
--------------------------------------------------------------…………


Step1.4 で求めた定常解の分岐図を Step1.3 と同様な方法で gnuplot で表示すると図5.2.5 の結果になる．

||u||∞
0.24

S"81_CP_TV_NR.dat" u 1:4

0.22

0.2

0.18

0.16

0.14

E2

0.12
1

1.5

2

2.5

3

3.5

図 5.2.5: 定常パルス波の分岐図

4. データの保存

67

4

4.5

d

step1.4 で求めた分岐図のデータを SP2 というディレクトリに保存する．
$ ./aout␣SP2



Step1.5 定常解の分岐計算
Step1.3 と Step1.4 で求めたデータを 1 つのデータに結合する．データを結合するだけなら$ cat␣SP1/middle_final.dat␣SP2/middle_final.dat > newf ilenameと実行することでデータをまとめることができる．しかし，今回の場合は図 5.2.4 と図 5.2.5で分岐計算を進めた方向が違うため，Step1.3 か Step1.4 で求めた middle_final.dat のどち
らかの step の進み方を逆にする必要がある．ここでは Step1.4 で求めた middle_final.datの step の進み方を逆にする．

1. Step1.4 で求めた middle_final.dat の step の進み方を逆にする

Sort は midle_final.dat の step の進み方を逆に並び替えるスクリプトである．
$ ./Sort と実行すると，middle_final.dat の step の進み方が逆になる．


midle_final.dat の中身が以下のようであったとする．
step 数
1
..
.
1
2
..
.
2
..
.
n
..
.
n


...
..
.
..
.
..
.
..
.
..
.
..
.
..
.
..
.
..
.
..
.

...
..
.
..
.
..
.
..
.
..
.
..
.
..
.
..
.
..
.
..
.

...
..
.
..
.
..
.
..
.
..
.
..
.
..
.
..
.
..
.
..
.

$ ./Sort  と実行することで，以下のように並び替わる．

68

step 数
n
..
.
n
..
.
2
..
.
2
1
..
.
1

...
..
.
..
.
..
.
..
.
..
.
..
.
..
.
..
.
..
.
..
.

...
..
.
..
.
..
.
..
.
..
.
..
.
..
.
..
.
..
.
..
.

...
..
.
..
.
..
.
..
.
..
.
..
.
..
.
..
.
..
.
..
.

ここでは Step1.4 のデータを逆の方向に並べ替えたいため，まず SP2 に保存した middle_final.dat

を持ってくる ($ cp␣SP2/middle_final.dat␣middle_final.dat  )．


そして$ ./Sort  を実行する．

2. データを結合する方法
一般に，f ilename1 と f ilename2 を結合して newf ilename というファイルに保存する場合には以下のコマンドを実行する．
$ cat␣f ilename1␣f ilename2 > newf ilename 

3 つのデータを結合させる場合


$ cat␣f ilename1␣f ilename2␣f ilename3 > newf ilename 

このコマンドを使い，順番を並べ替えた SP2 の middle_final.dat と，SP1 の middle_final.datのデータを結合し適当なファイル (ここでは tmp.dat ) に保存するには，$ cat␣middle_final.dat␣SP1/middle_final.dat > tmp.dat と実行する．

3. 結合したデータのステップ数を修正する方法
4. parameter.dat




…………
; StepCut…………
; Initial1

9; 読み飛ばすデータステップ数;

(a) 結合したデータの間引き間隔を parameter.dat の StepCut で変更する．1 とするとデータ間引きを行わない．ここでは 1 にしておく．
(b) 結合したデータのファイル名を middle_final.dat に変更する．
$ mv␣tmp.dat␣middle_final.dat (c) parameter.dat の Initial の値を 9 に変更する．
(d) $ ./aout␣1  と実行する．
69

実行後，結合した middle_final.dat に対応する• 81_CP_TV_NR.dat• 82_S_CP_TV_NR_REV_IMV_ABV_NUEV.dat• solution_value.dat• 1• 11が出力される．ただし，1 はステップ数が修正された middle_final.dat ，11 はステッ
プ数が修正された 9middle_final.dat である．
結合した分岐図を Step1.3 と同様な方法で gnuplot で表示すると図 5.2.6 のようになる．

||u||∞
0.24"81_CP_TV_NR.dat" u 1:4

0.22

0.2

0.18

0.16

E1

0.14

E2
0.12
1

1.5

2

2.5

3

3.5

4

4.5

d

図 5.2.6: 定常パルス波の分岐図

縦軸を ||u||L2 にすると図 5.2.7 のようになる．

||W|| L2
2.47

E2

2.46

2.45

2.44

2.43

E1

2.42

2.41

2.4

2.39

2.38

2.37
1

1.5

2

2.5

3

3.5

図 5.2.7: 定常パルス波の分岐図

70

4

4.5

d

1 ，11 をそれぞれ middle_final.dat ，9middle_final.dat に上書きし，データを SP というディレクトリに保存する．
$ ./aout␣SP

1 つのデータにまとめたので SP1 と SP2 のデータは$ rm␣-rf␣SP1␣SP2 で削除する．
Step1.6 定常解の分岐計算
ここでは Step1.5 で求めた定常解の安定性を求める方法について説明する．解の安定性は，⃝
1 解の枝を求めながら安定性を求める方法と⃝2 解の枝を求め終わった後、それらの解の安定性を調べる方法の 2 通りある．ここでは，定常解の枝を Step1.5 で既に求めたので 2 の方法を用いる

1. 実行前の設定(a) Initial.dat
Step1.5 で保存したディレクトリ SP の中から安定性計算に必要とする 9middle_final.dat
を Initial.dat に変更して計算するディレクトリに持ってくる．
$ cp␣SP/9middle_final.dat␣Initial.dat(b) parameter.dat




…………
8
; Initial; (∗1)(∗1) 各 step における解の固有値と固有関数の計算を行う
と変更する．

2. 実行

$ ./aout␣1  で実行する．
(42 ページにも解説がある．)



デモプログラム 1.5
Step1.6 の設定から実行までの操作は
$ cp␣../DEMO/Gray-Scott_Model/AS/05_SP3/*␣./ $ ./aout␣1 
で行える．

RD/DEMO/Gray-Scott_Model/AS/05_SP3 には Step1.6 の設定で必要なファ
イル．
• Function.c• Initial.dat• parameter.datが存在する．


3. 実行結果安定性の結果は 82_S_CP_TV_NR_REV_IMV_ABV_NUEV.dat のデータを参照する．gnuplot を開き$ gnuplot > plot␣"82_S_CP_TV_NR_REV_IMV_ABV_NUEV.dat"␣u␣1:8 と71

実行することで横軸にステップ数，縦軸に固有値実部を与えたグラフが表示される (図5.2.8)．
固有値実部
0.07"82_S_CP_TV_NR_REV_IMV_ABV_NUEV.dat" u 1:8
0.06

0.05

0.04

0.03

0.02

0.01

0

-0.01

-0.02
0

20

40

60

80

100

120

140

160

180

200

step数

図 5.2.8:

step の変化と共に固有値の実部が 0 を横切るとき，分岐が起きる．この場合，6 つの分岐が
存在する．
固有値実部
0.07
"82_S_CP_TV_NR_REV_IMV_ABV_NUEV.dat" u 1:8
0.06

0.05

0.04

0.03

0.02

4

2 3

1

0.01

0

-0.01

5

6

120

140

-0.02
0

20

40

60

80

100

160

180

200

step数

図 5.2.9:

図 5.2.9 の 1〜6 の分岐点の種類について述べる．

1 の分岐点について
1 の分岐点は step= 70 の近くで起こり，82_S_CP_TV_NR_REV_IMV_ABV_NUEV.dat のstep= 70 の近傍では次の結果となる．

72





左から 1 列目の step 数，8 列目の固有値実部，9 列目の固有値虚部の値．

69
69
69
69
69
70
70
70
70
70

……………………
……………………
……………………
……………………
……………………
……………………
……………………
……………………
……………………
……………………

…………
2.72392992E-02 0.00000000E+00 …………
6.77088438E-09 0.00000000E+00 …………
-1.38528671E-03 0.00000000E+00 …………
-1.62637187E-02 0.00000000E+00 …………
-1.80527589E-02 0.00000000E+00 …………
2.38602696E-02 0.00000000E+00 …………
7.67131139E-04 0.00000000E+00 …………
6.04405127E-09 0.00000000E+00 …………
-1.57752778E-02 0.00000000E+00 …………
-1.80527990E-02 0.00000000E+00 …………
…………





step= 70 近くで固有値実部が 1 つ 0 を横切ることがわかる．ただし，定常解の周期境界条件の影響でシフトの自由度があるため，固有値のひとつは零固有値を持つ．この零固有値は安定性の条件に影響を与えない．固有値実部の 1 つが 0 を横切るとき，次の 3 つのタイプの分岐が現れる:
saddle-node 分岐
pitch-fork 分岐transcritical 分岐
このうち step の変化と共に分岐点近傍でコントロールパラメータの進む向きが変わる
とき，saddle-node 分岐が起きる．gnuplot を開き，

$ gnuplot > plot␣"82_S_CP_TV_NR_REV_IMV_ABV_NUEV.dat"␣u␣1:2と実行して横軸に step 数，縦軸にコントロールパラメータ d を表示させる (図 5.2.10)．

d4.5"82_S_CP_TV_NR_REV_IMV_ABV_NUEV.dat" u 1:2

4

3.5

3

2.5

2

1.5

1
0

20

40

60

80

100

120

140

160

180

200

step数

図 5.2.10:
このとき step= 70, 110, 141 でコントロールパラメータの向きが変化することがわかる．
すなわち，step= 70, 110, 141 で saddle-node 分岐が起きる．従って，図 5.2.9 の 1 はsaddle-node 分岐点である．同様に，図 5.2.9 の 4 と 6 も saddle-node 分岐点であること

73

がわかる (図 5.2.11)．

||u||∞
0.24"81_CP_TV_NR.dat" u 1:4

110step

0.22

0.2

141step

0.18

0.16

70step
0.14

0.12
1

1.5

2

2.5

3

3.5

4

4.5

d

図 5.2.11:2 の分岐点について
2 の分岐点は step= 96 の近くで起こり，82_S_CP_TV_NR_REV_IMV_ABV_NUEV.dat のstep= 96 をみると




左から 1 列目の step 数，8 列目の固有値実部，9 列目の固有値虚部の値．

95
95
95
95
95
96
96
96
96
96

……………………
……………………
……………………
……………………
……………………
……………………
……………………
……………………
……………………
……………………

…………
3.39914803E-04 2.62426635E-02 …………
3.39914803E-04 -2.62426635E-02 …………
1.89058850E-08 0.00000000E+00 …………
-5.28527110E-03 0.00000000E+00 …………
-1.80952017E-02 0.00000000E+00 …………
2.04794297E-08 0.00000000E+00 …………
-4.20137537E-05 -2.64192278E-02 …………
-4.20137537E-05 2.64192278E-02 …………
-4.86485952E-03 0.00000000E+00 …………
-1.80979881E-02 0.00000000E+00 …………
…………





となる．step の変化と共に 2 つの固有値が固有値虚部を持ちながら複素平面の虚軸を横切ることがわかる．この場合，Hopf 分岐が起きる．同様に，図 5.2.9 の 5 の場合も
82_S_CP_TV_NR_REV_IMV_ABV_NUEV.dat の step= 115 をみると

74





左から 1 列目の step 数，8 列目の固有値実部，9 列目の固有値虚部の値．

114
114
114
114
114
115
115
115
115
115

……………………
……………………
……………………
……………………
……………………
……………………
……………………
……………………
……………………
……………………

…………
9.09383226E-03 0.00000000E+00 …………
6.06194605E-03 0.00000000E+00 …………
-7.51816727E-09 0.00000000E+00 …………
-3.51838009E-05 -2.54489667E-02 …………
-3.51838009E-05 2.54489667E-02 …………
9.93006414E-03 0.00000000E+00 …………
7.15965142E-03 0.00000000E+00 …………
6.13616555E-04 2.47780506E-02 …………
6.13616555E-04 -2.47780506E-02 …………
-5.09436851E-09 0.00000000E+00 …………
…………





同じことが起きるため，Hopf 分岐が起きている．

3 の分岐点について
3 の分岐点は step= 107 の近くで起こり，82_S_CP_TV_NR_REV_IMV_ABV_NUEV.dat を
みると，




左から 1 列目の step 数，8 列目の固有値実部，9 列目の固有値虚部の値．

106
106
106
106
106
107
107
107
107
107

……………………
……………………
……………………
……………………
……………………
……………………
……………………
……………………
……………………
……………………

…………
7.19556478E-04 0.00000000E+00 …………
-1.35210633E-08 0.00000000E+00 …………
-3.01568890E-03 2.77459724E-02 …………
-3.01568890E-03 -2.77459724E-02 …………
-1.26064364E-02 4.66843397E-03 …………
1.55106188E-03 0.00000000E+00 …………
7.65809216E-09 0.00000000E+00 …………
-3.06796344E-03 2.78066326E-02 …………
-3.06796344E-03 -2.78066326E-02 …………
-1.17306086E-02 1.92013563E-03 …………
…………





となる．step の変化と共に 1 つの固有値が 0 を横切ることがわかる．この場合は，先程と同様に 3 つの分岐が現れるが，図 4.2.9 の情報から saddle-node 分岐ではないことがわかるため，pitch-fork 分岐か transcritical 分岐である．これ以上の判別は固有値の情報から理解できない．これを判別するには，3 の分岐点から伸びる解の枝を求める必要がある．詳しい説明は省略するが，3 は pitch-fork 分岐点であることが説明できる．
以上の結果から，図 5.2.9 の 1〜6 では次の 3 つの種類の分岐が現れることがわかる．

• saddle-node 分岐 (1,4,6)• Hopf 分岐 (2,5)• pitch-fork 分岐 (3)
データの保存75

Step1.6 で求めた分岐図のデータを SP というディレクトリに上書きする．
$ ./aout␣SP

このとき，ディレクトリを上書きしますかと尋ねられるので，$ yes と入力する．

Step2.1 進行解の数値計算
ここでは Step1.6 で求めた定常解の pitch-fork 分岐点における固有関数の情報を用いて進行解を求める．用いる固有関数の情報は，固有値の実部が正である step のものを用いる．Step1.6
でも示したように，pitch-fork 分岐点での固有値の情報は以下のようになっている．


左から 1 列目の step 数，8 列目の固有値実部，9 列目の固有値虚部の値．

106
106
106
106
106
107
107
107
107
107

……………………
……………………
……………………
……………………
……………………
……………………
……………………
……………………
……………………
……………………

…………
7.19556478E-04 0.00000000E+00 …………
-1.35210633E-08 0.00000000E+00 …………
-3.01568890E-03 2.77459724E-02 …………
-3.01568890E-03 -2.77459724E-02 …………
-1.26064364E-02 4.66843397E-03 …………
1.55106188E-03 0.00000000E+00 …………
7.65809216E-09 0.00000000E+00 …………
-3.06796344E-03 2.78066326E-02 …………
-3.06796344E-03 -2.78066326E-02 …………
-1.17306086E-02 1.92013563E-03 …………
…………





step106 では負の値だった 2 行目の固有値実部が，step107 では正の値になっている．したがって，ここでは step107 の固有関数の情報を用いる．
1. 実行前の設定(a) Initial.dat
pitch-folk 分岐点の固有関数を Initial.dat として使用する．固有ベクトルのデー
タは 6_EigenVector_Data に入っている．
6_EigenVector_Data/ 107_CP . . . 5 /o1 . . . 6 のデータを Initial.dat に変更して持ってくる．
(b) parameter.dat

5
6

ここにはコントロールパラメータの値が入っている
ここには固有値実部の値が入っている

76





ds
; パラメータの変化量なし
KellerEqn; 擬似弧長法の設定なし…………
; 進行解を選択
2; SolSelect…………
1; MaxStep; 最大 step 回数…………
2; Initial; (∗1)(∗1) 分岐点から伸びる解の枝を求める
0.000000
-1

;
;





と変更する．

(c) 0_BranchChangeParameter.dat
分岐点における零固有値に対する固有関数の情報と 0_BranchChangeParameter.datのパラメータを用いて Newton 法に用いる初期値 (uini , cini , pini ) を構成する.
• 初期値の構成 uini = u0 + εΦR ,• 速度の構成 cini = δc ,• パラメータの構成 pini = δp ,ただし，u0 は定常解の pitch-fork 分岐点における解 u，ΦR は分岐点における零固有値に対する固有関数の実部，ε, δc , δp はそれぞれ 0_BranchChangeParameter.dat
における eps , velocity_delta , par_delta である.
ここでは，






0.08
0.00
0.0
0.00000000
0.000

;
;
;
;
;

eps
period_delta
velocity_delta
par_deltatimeと変更する．

2. 実行

$ ./aout␣1  で実行する

77; uini の構成に必要; ここでは不必要な設定; 速度の微小変量; パラメータの微小変量; ここでは不必要な設定







デモプログラム 2.1
Step2.1 の設定から実行までの操作は
$ cp␣../DEMO/Gray-Scott_Model/AS/06_TP1/*␣./ $ ./aout␣1 
で行える．

RD/DEMO/Gray-Scott_Model/AS/06_TP1 には Step2.1 の設定で必要なファ
イル．
• Function.c• Initial.dat• parameter.dat• 0_BranchChangeParameterが存在する．


3. 実行結果


…………

Step / MaxStep
iter / MaxIter
Control Parameter
Velocity
MachineEps
CC_Norm

=
=
=
=
=
=

1 / 1
6 / 302.300901915266015
0.000337775542350
0.000000010000000
0.000000043803793

Step / MaxStep
iter / MaxIter
Control Parameter
Velocity
MachineEps
CC_Norm

=
=
=
=
=
=

1 / 1
7 / 302.300901915266015
0.000337643953297
0.000000010000000
0.000000015729003

Step / MaxStep=
1 / 1
iter / MaxIter=
8 / 30Control Parameter=
2.300901915266015
Velocity=
0.000337596714033
MachineEps=
0.000000010000000
CC_Norm=
0.000000005653990--------------------------------------------------------------Convergence !!
---------------------------------------------------------------



78



収束した解は速度を持っているため，進行解であるとわかる．

Newton 法による計算が発散したり，進行解と異なる別の解が求まったときは，0_BranchChangeParameterの各パラメータを調整 (eps を大きくするなど) して，進行解が求まるまで再実行する．
データの保存Step2.1 で求めた分岐図のデータを TP_INI というディレクトリに保存する．


$ ./aout␣TP_INI

Step2.2 進行解の分岐計算
Step2.1 で収束した解を用いてパラメータを変化させながら進行解の枝を求める．
1. 実行前の設定(a) Initial.dat
初期値は Step2.1 で収束させた解を用いる．
(b) parameter.dat




0.01

;

200
1

;
;

1

;



; パラメータの変化量 0.01

ds…………
MaxStep
StepCut…………
Initial; 最大 step 回数; 1step ごとにデータを出力する; 分岐図を求めることを選択と変更する．

2. 実行$ cp␣TP_INI/middle_final.dat␣Initial.dat





$ ./aout␣1 

デモプログラム 2.2
Step2.2 の設定から実行までの操作は
$ cp␣../DEMO/Gray-Scott_Model/AS/07_TP2/*␣./ $ ./aout␣1 
で行える．

RD/DEMO/Gray-Scott_Model/AS/07_TP2 には Step2.2 の設定で必要なファイル• Function.c• Initial.dat• parameter.datが存在する．


3. 実行結果

79





…………
Step / MaxStep
iter / MaxIter
Control Parameter
Velocity
MachineEps
CC_Norm

=
=
=
=
=
=

200 / 200
0 / 304.300901915265973
0.001219755984295
0.000000010000000
0.000036275013261

Step / MaxStep
iter / MaxIter
Control Parameter
Velocity
MachineEps
CC_Norm

=
=
=
=
=
=

200 / 200
1 / 304.300901915265973
0.001220913746472
0.000000010000000
0.000000033525231

Step / MaxStep=
200 / 200
iter / MaxIter=
2 / 30Control Parameter=
4.300901915265973
Velocity=
0.001220904601074
MachineEps=
0.000000010000000
CC_Norm=
0.000000000157812--------------------------------------------------------------Convergence !!
---------------------------------------------------------------



進行解の分岐図を表示すると図 5.2.12 のようになる．

||u||∞
0.36’TP3/81_CP_TV_NR.dat’ u 1:4

0.34

0.32

0.3

進行解の分岐図

0.28

0.26

S
0.24
2.2

2.4

2.6

2.8

3

3.2

3.4

3.6

3.8

4

4.2

4.4

d

図 5.2.12: 進行解の分岐図
図 5.2.12 を定常解とあわせて表示すると次のようになる (※はデモでは求めていない部
分の進行解の枝)．

80

||u||∞
0.35’SP/81_CP_TV_NR.dat’ u 1:4’TP3/81_CP_TV_NR.dat’ u 1:4

0.3

進行解の分岐図

0.25

※

0.2

定常解の分岐図

0.15

0.1
1

1.5

2

2.5

3

3.5

4

4.5

d

図 5.2.13: 進行解の分岐図

Step1.6 と Step2.3 で求める安定性を考慮し，縦軸を ||u||L2 にすると図 5.2.14 のようになる．

||W|| L2
2.46

定常解の分岐図

2.44

2.42

進行解の分岐図

2.4

2.38

2.36

d

2.34
1

1.5

2

2.5

3

3.5

4

4.5

図 5.2.14: 進行解の分岐図

4. データの保存
Step2.2 で求めた分岐図のデータを TP というディレクトリに保存する．
$ ./aout␣TP

Step2.3 進行解の分岐計算
ここでは Step2.2 で求めた進行解の安定性を求める方法について説明する．
1. 実行前の設定(a) Initial.dat
まず，Step2.2 で保存したディレクトリ TP の中から固有値計算に必要なデータ
9middle_final.dat を Initial.dat に変更して計算するディレクトリに持ってくる．
81

$ cp␣TP/9middle_final.dat␣Initial.dat(b) parameter.dat


8





;

…………
Initial; 安定性の計算を選択と変更する．

2. 実行

$ ./aout␣1  で実行する．


デモプログラム 2.3
Step2.3 の設定から実行までの操作は
$ cp␣../DEMO/Gray-Scott_Model/AS/08_TP3/*␣./ $ ./aout␣1 
で行える．

RD/DEMO/Gray-Scott_Model/AS/08_TP3 には Step2.3 の設定で必要なファイル• Function.c• Initial.dat• parameter.datが存在する．


3. 実行結果
82_S_CP_TV_NR_REV_IMV_ABV_NUEV.dat のデータから横軸にステップ数，縦軸に固有値実部をグラフで示すと図 5.2.15 となる．
固有値実部
0.005
’TP3/82_S_CP_TV_NR_REV_IMV_ABV_NUEV.dat’ u 1:8

pitch-folk分岐点
0

Hopf分岐点-0.005-0.01-0.015

-0.02
0

20

40

60

80

100

120

図 5.2.15:

82

140

160

180

200

step数

進行解の安定性は 128 ページの定常解の安定性と同様に調べる．デモプログラムを使用した場合は，106step で 82_S_CP_TV_NR_REV_IMV_ABV_NUEV.dat の固有値の実部が2 つ 0 を横切るため Hopf 分岐が存在する．Hopf 分岐点からは脈動進行解の枝が存在
する．




左から 1 列目の step 数，8 列目の固有値実部，9 列目の固有値虚部の値．

105
105
105
105
105
106
106
106
106
106

……………………
……………………
……………………
……………………
……………………
……………………
……………………
……………………
……………………
……………………

…………
1.02512634E-08 0.00000000E+00 …………
-7.69879710E-06 -1.61517239E-02 …………
-7.69879710E-06 1.61517239E-02 …………
-1.53054324E-02 -3.86482927E-02 …………
-1.53054324E-02 3.86482927E-02 …………
3.08429993E-05 1.61766085E-02 …………
3.08429993E-05 -1.61766085E-02 …………
1.02497315E-08 0.00000000E+00 …………
-1.53711572E-02 3.87068421E-02 …………
-1.53711572E-02 -3.87068421E-02 …………
…………





4. データの保存
Step2.3 で求めた分岐図のデータを TP というディレクトリに上書きする．
$ ./aout␣TP

このとき，ディレクトリを上書きしますかと尋ねられるので，$ yes と入力する．

Step3.1 脈動進行解の分岐計算
ここでは Step2.3 で求めた進行解の Hopf 分岐点の固有関数の情報を用いて脈動進行解を求める．
1. 実行前の設定(a) Initial.dat
Hopf 分岐点の固有関数を Initial.dat として使用する．固有ベクトルのデータ
は 6_EigenVector_Data に入っている．6_EigenVector_Data/ 106_CP . . . 7 /o1. . . 8 のデータを Initial.dat に変更して持ってくる．
(b) parameter.dat

78

ここにはコントロールパラメータの値が入っている
ここには固有値実部の値が入っている

83





ds
; パラメータの変化量なし…………
4; SolSelect; 脈動進行解を選択…………
1; StepCut; 毎 step ごとにデータを出力する
1; MaxStep; 最大 step 回数…………
2; Initial; (∗1)(∗1) 分岐点から伸びる解の枝を求める．
0.000000;と変更する．

(c) 0_BranchChangeParameter.dat
分岐点における零固有値に対する固有関数の情報と 0_BranchChangeParameter.datのパラメータを用いて Newton 法に用いる初期値 (uini , Tini , cini , pini ) を構成する.
• 初期値の構成 uini = u0 + ε(sin(2πt)ΦR + cos(2πt)ΦI ),2π• 周期の構成 Tini =+ δTβ• 速度の構成 cini = δc ,• パラメータの構成 pini = δp ,ただし，u0 は定常解の pitch-fork 分岐点における解 u，ΦR , ΦI はそれぞれ分岐点における零固有値に対する固有関数の実部と虚部，c は解の速度，T は解の周期 p はコン
トロールパラメータ，c0 は進行解の Hopf 分岐点における速度，ε, δT , δc , δp , t はそれ
ぞれ 0_BranchChangeParameter.dat における eps , period_delta , velocity_delta, par_delta , time である.
ここでは，




0.05
0.00
0.0
0.00000000
0.000

;
;
;
;
;

eps
period_delta
velocity_delta
par_delta
time



SS_2
SS_2
SS_2

&
&
&
&
&

SS_3
SS_3
SS_3
SS_3

&
&
&
&
&

SS_4
SS_4
SS_4
SS_4
SS_4


と変更する．

2. 実行

;
;
;
;
;



$ ./aout␣1  で実行する

84





デモプログラム 3.1
Step3.1 の設定から実行までの操作は
$ cp␣../DEMO/Gray-Scott_Model/AS/09_TB1/*␣./ $ ./aout␣1 
で行える．

RD/DEMO/Gray-Scott_Model/AS/09_TB1 には Step3.1 の設定で必要なファ
イル．
• Function.c• Initial.dat• parameter.dat• 0_BranchChangeParameterが存在する．


3. 実行結果

85





…………

Step / MaxStep
iter / MaxIter
Control Parameter
PeriodVelocity
MachineEps
CC_Norm

=
=
=
=
=
=
=

1 / 1
7 / 303.360901915265992
394.759258820564980
0.001064648438453
0.000000010000000
0.000000060381561

Step / MaxStep
iter / MaxIter
Control Parameter
PeriodVelocity
MachineEps
CC_Norm

=
=
=
=
=
=
=

1 / 1
8 / 303.360901915265992
394.759140363937661
0.001064648372300
0.000000010000000
0.000000013427036

Step / MaxStep=
1 / 1
iter / MaxIter=
9 / 30Control Parameter=
3.360901915265992
Period=394.759166707443910
Velocity=
0.001064648387011
MachineEps=
0.000000010000000
CC_Norm=
0.000000002985457--------------------------------------------------------------Convergence !!
--------------------------------------------------------------…………


ただし，プログラムのバグで「Convergence !!」と表示されていても，周期 (Period) が負になっていたり，極端に大きな値になっていたりする場合は，求めたい解に収束していないことに注意する．

Step3.2 脈動進行解の分岐計算
定常解，進行解のときと同様に脈動進行解の分岐図を求める．
脈動進行解の枝を追うのは非常に難しい．小まめにデータを保存しながら少しずつ進めてい
く必要がある．そのためデモプログラムは省略する．
分岐図は以下のようになる．

Step4.1 脈動定常解の分岐計算
ここでは，Step1.6 で求めた定常解の Hopf 分岐点における固有関数の情報を用いて脈動定常解を求める．方法は Step2.1 とほぼ同様であるが，固有関数は保存されていないため求めな
86

||u||∞
0.314"81_CP_TV_NR.dat" u 1:4

0.313

0.312

0.311

脈動進行解の分岐図

0.31

0.309

0.308

0.307
3.355

3.36

3.365

3.37

3.375

3.38

3.385

3.39

3.395

3.4

d

4.4

d

図 5.2.16:||u||∞
0.36"81_CP_TV_NR.dat" u 1:4"../8_TP3/81_CP_TV_NR.dat" u 1:4
0.34

脈動進行解の分岐図
0.32

0.3

0.28

進行解の分岐図

0.26

0.24

0.22
2.2

2.4

2.6

2.8

3

3.2

3.4

3.6

3.8

4

4.2

図 5.2.17:

おす必要がある．Step4.1 では，Hopf 分岐点における固有関数を求めなおす．まず，Step1.6で保存したデータを AS にもってくる．
$ cp␣SP/*␣./ Step1.6 でも示したように，最初に現れる Hopf 分岐点での固有値の情報は以下のようになっ
ている．

87





左から 1 列目の step 数，8 列目の固有値実部，9 列目の固有値虚部の値．

95
95
95
95
95
96
96
96
96
96

……………………
……………………
……………………
……………………
……………………
……………………
……………………
……………………
……………………
……………………

…………
3.39914803E-04 2.62426635E-02 …………
3.39914803E-04 -2.62426635E-02 …………
1.89058850E-08 0.00000000E+00 …………
-5.28527110E-03 0.00000000E+00 …………
-1.80952017E-02 0.00000000E+00 …………
2.04794297E-08 0.00000000E+00 …………
-4.20137537E-05 -2.64192278E-02 …………
-4.20137537E-05 2.64192278E-02 …………
-4.86485952E-03 0.00000000E+00 …………
-1.80979881E-02 0.00000000E+00 …………
…………





95step では正の値であった 1 列目と 2 列目の固有値実部の値が，96step では 2 列目と 3 列目に負の値として記録されていることがわかる．したがって，ここでは 95step の固有関数を求める．
1. 実行前の設定(a) Initial.dat
Hopf 分岐点の 9middle_final.dat を Initial.dat として使用する．9middle_final.dat
を任意のテキストエディタで開き，95step 以外の全ての値を消去し，保存する．その後，9middle_final.dat を Initial.dat に変更する ($ cp␣9middle_final.dat
␣Initial.dat  )．
(b) parameter.dat




…………
8
; Initial; (∗1)(∗1) 各 step における解の固有値と固有関数の計算を行う
と変更する．

2. 実行

$ ./aout␣1  で実行する

88





デモプログラム 4.1
Step4.1 の設定から実行までの操作は
$ cp␣../DEMO/Gray-Scott_Model/AS/11_SB1/*␣./ $ ./aout␣1 
で行える．

RD/DEMO/Gray-Scott_Model/AS/11_SB1 には Step4.1 の設定で必要なファ
イル．
• Function.c• Initial.dat• parameter.datが存在する．


3. 実行結果…………

---------------------------------------------------------------Step / MaxStep= 1 / 1
Control Parameter=
1.869650891473
Nearly Most Zero EigenValue =0.000000018906- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - REV[1] = 0.000339914803 IEV[1] = 0.026242663526
RI[1] = 0.026244864851
REV[2] = 0.000339914803 IEV[2] = -0.026242663526 RI[2] = 0.026244864851
REV[3] = 0.000000018906 IEV[3] = 0.000000000000
RI[3] = 0.000000018906----------------------------------------------------------------



Step4.2 脈動定常解の分岐計算
ここでは，Step4.1 で求めた定常解の Hopf 分岐点における固有関数の情報を用いて脈動定常解を求める．
1. 実行前の設定(a) Initial.dat
Hopf 分岐点の固有関数を Initial.dat として使用する．固有ベクトルのデータ
は 6_EigenVector_Data に入っている．
6_EigenVector_Data/ 1_CP . . . 9 /o1 . . . 10 のデータを Initial.dat に変更して持ってくる．
(b) parameter.dat

910

ここにはコントロールパラメータの値が入っている
ここには固有値実部の値が入っている

89





ds
; パラメータの変化量なし
KellerEqn; 擬似弧長法の設定あり…………
; 脈動定常解を選択
3; SolSelect…………
1; MaxStep; 最大 step 回数…………
2; Initial; (∗1)(∗1) 分岐点から伸びる解の枝を求める
0.000000
1

;
;





と変更する．

(c) 0_BranchChangeParameter.dat
分岐点における零固有値に対する固有関数の情報と 0_BranchChangeParameter.datのパラメータを用いて Newton 法に用いる初期値 (uini , pini ) を構成する.
• 初期値の構成 uini = u0 + ε(sin(2πt)ΦR + cos(2πt)ΦI ),• 周期の構成 Tini = 2π
β + δT ,• パラメータの構成 pini = δp ,ただし，u0 は定常解の Hopf 分岐点における解 u，ΦR , ΦI はそれぞれ分岐点における零固有値に対する固有関数の実部と虚部，T は解の周期 p はコントロールパ
ラメータ，β は定常解の Hopf 分岐点における固有値虚部の値，ε, δT , δp , t はそれぞ
れ 0_BranchChangeParameter.dat における eps , period_delta , par_delta ,time である. β:定常解の Hopf 分岐点における固有値虚部の値
ここでは，






0.05
0.00
0.0
0.00000000
0.000

;
;
;
;
;

eps
period_delta
velocity_delta
par_deltatimeと変更する．

2. 実行

$ ./aout␣1  で実行する

90; uini の構成に必要; 周期の微小変量; ここでは不必要な設定; パラメータの微小変量; uini の構成に必要







デモプログラム 4.2
Step4.2 の設定から実行までの操作は
$ cp␣../DEMO/Gray-Scott_Model/AS/12_SB2/*␣./ $ ./aout␣1 
で行える．

RD/DEMO/Gray-Scott_Model/AS/12_SB2 には Step4.2 の設定で必要なファ
イル．
• Function.c• Initial.dat• parameter.dat• 0_BranchChangeParameterが存在する．


3. 実行結果


…………

Step / MaxStep
iter / MaxIter
Control Parameter
PeriodMachineEps
CC_Norm

=
=
=
=
=
=

1 / 1
4 / 301.994435415511887
259.665896228138479
0.000000010000000
0.000004054062691

Step / MaxStep
iter / MaxIter
Control Parameter
PeriodMachineEps
CC_Norm

=
=
=
=
=
=

1 / 1
5 / 301.994440834482065
259.670083565151799
0.000000010000000
0.000000050642485

Step / MaxStep=
1 / 1
iter / MaxIter=
6 / 30Control Parameter=
1.994440898624708
Period=259.670129810254650
MachineEps=
0.000000010000000
CC_Norm=
0.000000000454148--------------------------------------------------------------Convergence !!
---------------------------------------------------------------



収束した解は周期を持っているため，脈動解であるとわかる．

91今回のように，脈動解の場合は擬似弧長法の設定を入れると解を収束させやすい．
データの保存Step4.2 で求めた分岐図のデータを SB_INI というディレクトリに保存する．


$ ./aout␣SB_INI

Step4.3 脈動定常解の分岐計算
Step4.2 で収束した解を用いてパラメータを変化させながら脈動定常解の枝を求める．
1. 実行前の設定(a) Initial.dat
初期値は Step4.2 で収束させた解を用いる．
(b) parameter.dat-0.00005

;

1
1

;
;

1
100

;
;

1

;



ds
…………
EigenValue
EigenVector…………
StepCut
MaxStep…………
Initial; パラメータの変化量 −0.00005; 固有値の計算あり; 固有ベクトル計算あり; 1step ごとにデータを出力する; 最大 step 回数; 分岐図を求めることを選択と変更する．

2. 実行$ cp␣SB_INI/middle_final.dat␣Initial.dat





$ ./aout␣1 

デモプログラム 4.3
Step4.3 の設定から実行までの操作は
$ cp␣../DEMO/Gray-Scott_Model/AS/13_SB3/*␣./ $ ./aout␣1 
で行える．

RD/DEMO/Gray-Scott_Model/AS/13_SB3 には Step4.3 の設定で必要なファイル• Function.c• Initial.dat• parameter.datが存在する．


3. 実行結果
途中で計算が止まってしまう．

92





…………

--------------------------------------------------------------= 50 / 100
Step / MaxStep
Control Parameter=
2.190080118610
Nearly Most One EigenValue=
0.961863982654- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - REV[1] = -1776.464390405427 IEV[1] = 0.000000000000
RI[1] = 1776.464390405427
REV[2] = 11.796518082181 IEV[2] = 0.000000000000
RI[2] = 11.796518082181
REV[3] = 2.095899744166 IEV[3] = 0.000000000000
RI[3] = 2.095899744166--------------------------------------------------------------Step / MaxStep
iter / MaxIter
Control Parameter
PeriodMachineEps
CC_Norm

=
=
=
=
=
=

51 / 100
0 / 502.190080118609975
469.503536462149782
0.000000010000000
0.000121647676807

Step / MaxStep
iter / MaxIter
Control Parameter
PeriodMachineEps
CC_Norm

=
=
=
=
=
=

51 / 100
1 / 502.184345104892049
481.258229677027430
0.000000010000000
0.060014529587429

Step / MaxStep
iter / MaxIter
Control Parameter
PeriodMachineEps
CC_Norm

=
=
=
=
=
=

51 / 100
2 / 502.171010326160396
489.541509841286711
0.000000010000000
0.235433907099954

Step / MaxStep=
51 / 100iter / MaxIter=
3 / 50Control Parameter=
1.698490166970823= -438.571793196803696
PeriodMachineEps=
0.000000010000000
CC_Norm=
0.000000000000000--------------------------------------------------------------Convergence !!
---------------------------------------------------------------



「Convergence !!」と表示されてはいるが，周期の値が負になっている．横軸に step 数，
縦軸に固有値の実部をとる (図 5.2.18) と，ひとつの固有値の実部が負の方向に大きく93

なりすぎているために計算が止まったのだとわかる．
固有値実部
200
"82_S_CP_TV_NR_REV_IMV_ABV_NUEV.dat" u 1:8
0
-200
-400
-600
-800
-1000
-1200
-1400
-1600
-1800
0

5

10

15

20

25

30

35

40

45

50

step数

図 5.2.18: 固有値実部
脈動定常解の分岐図は以下のようになる．

||u||∞
0.23"81_CP_TV_NR.dat" u 1:4

脈動定常解の分岐図
0.225

0.22

0.215

0.21

0.205
1.95

2

2.05

2.1

2.15

2.2

図 5.2.19: 定常解と脈動定常解の分岐図

94

2.25

d

図 5.2.19 を定常解とあわせて表示すると次のようになる (※はデモでは求めていない部
分の脈動定常解の枝)．

||u||∞
0.24"81_CP_TV_NR.dat" u 1:4"SP3/81_CP_TV_NR.dat" u 1:4

0.22

※

脈動定常解の分岐図

0.2

0.18

定常解の分岐図

0.16

0.14

0.12
1

1.5

2

2.5

3

3.5

4

4.5

d

図 5.2.20: 定常解と脈動定常解の分岐図

4. データの保存
Step4.3 で求めた分岐図のデータを SB というディレクトリに保存する．
$ ./aout␣SB

Appendix1.1 脈動進行解の分岐計算
Appendix1.1 では，5.1 節で求めた Step3.2 の解を用いて，脈動進行解を求める．
時間発展解で求めた脈動進行解を Newton 法の初期値として与える場合は脈動進行解の速度と周期を求める必要がある．速度は適切な Velocity_Site と Velocity_Measure_Start(16 ページ参照) を与えて求め，周期は 27 ページで説明した方法を用いる．Appendix1.1 では，この一連の作業について説明する．


まず，カレントディレクトリを時間発展方程式を求めるところに移動 ($ cd␣RD/TE  )し，3.2 節の Step3.2 で保存した全てのデータ TB/*.dat をカレントディレクトリに持ってくる ($ cp␣TB/*␣./ )．これによって，3.2 節の Step3.2 を実行したあとの状態が復元できる．その後 final.dat 開くと，下から 5 行の 1 列目と 5 列目は以下のようになっていることがわかる．




…………

周期
速度





1. 解の速度について
適切な Velocity_Site と Velocity_Measure_Start を選択してあれば final.dat の
左から 5 列目，下から 2 番目のところに速度の値が入っている．

95

(a) Velocity_Site
脈動進行解は解が上下に振動するため，パルスの高さに注意して Velocity_Site
はパルス波のフロントとバックが常に存在する位置に設定する．デモでは 0.2 と与えている．
(b) Velocity_Measure_Start
Velocity_Measure_Start は MaxStep より小さな値を与える．
2. 解の周期について
final.dat の左から 5 列目，下から 3 番目には脈動進行解の周期を入れる必要がある．
(a) extremum を使う
$ ./extremum と実行すると，以下のような画面が表示される．

||u||L2

2.592"81_Step_Time_Sup_L1_L2.dat" u 2:5"EXTREMUM_DATA/1_LocalMaxData.dat" u 2:5"EXTREMUM_DATA/5_LocalMinData.dat" u 2:5

2.59
2.588
2.586
2.584
2.582
2.58
2.578
2.576
2.574
2.572
0

500

1000

1500

2000

2500

3000

d

図 5.2.21:
しかし，この図では極大値と極小値が求められていない．これは StepCut が大きいゆえに，時系列データが粗くて極大値と極小値が求められていない (脈動解の周期が求められていない)．したがって，StepCut を小さくし，3.2 節の Step3.2 を再計算する．StepCut を小さくすると gnuplot 画面の表示が追いつかないので，gnuplot
画面の表示を変更する．

(b) StepCut を調整して再計算する
parameter.dat を

96





…………
200
; Gnuplot_CutN_Vel_Dis…………
200; StepCut…………
1245; Gnuplot_Select…………
(∗1) 速度データを 200 回おきに出力する; (∗1); 200 回おきにデータを出力; gnuplot 画面の 3,6 番を表示と変更する．Initial.dat は 3.2 節の Step3.2 で使用した Initial.dat を使用
する．
$ ./aout␣1 で実行する．


デモプログラム A1.1
Appendix1.1(b) の設定から実行までの操作は$ cp␣../DEMO/Gray-Scott_Model/TE/05_TB3/*␣./





$ ./aout␣1 
で行える．

RD/DEMO/Gray-Scott_Model/TE/05_TP3 には Appendix1.1(b) の設定
で必要なファイル．
• 0parameter.dat• Function.c• parameter.dat• Initial.dat• gnuplot.cが存在する．
(c) extremum を使う
再び extremum を使うと以下の図が表示される．

97

||u||L2

2.595"81_Step_Time_Sup_L1_L2.dat" u 2:5"EXTREMUM_DATA/1_LocalMaxData.dat" u 2:5"EXTREMUM_DATA/5_LocalMinData.dat" u 2:5
2.59

2.585

2.58

2.575

2.57

2.565
0

500

1000

1500

2000

2500

3000

d

図 5.2.22:

4_LocalMaxAvePeriod.dat を開くと極大値の時間間隔値が，

423.3333333333



とわかるので，この値を脈動解の周期に用いる．極小値の時間間隔値 8_LocalMinAvePeriod.dat
を用いることもできる．

(d) final.dat に入力
final.dat の左から 5 列目，下から 3 番目の 0.0 を周期の値に書き換える．




…………
30000 ………… ………… ………… 9.999979890995608E-01
30000 ………… ………… ………… 9.999978794605785E-01
30000 ………… ………… ………… 423.3333333333
30000 ………… ………… ………… 1.074709848399322E-03
30000 ………… ………… ………… 0.000000000000000E+00


Appendix1.2 脈動進行解の分岐計算
Appendix1.1 で求めた final.dat を Newton 法の初期値に与え脈動進行解を求める．カレントディレクトリを分岐図を求めるところに移動 ($ cd␣../AS  ) する．
1. Initial.dat
時間発展方程式で求めた final.dat を Initial.dat として使用する.

2. parameter.dat

98





0.000000

;

4

;

1
1

;
;

1

;



; パラメータの変化量なし

ds…………
SolSelect…………
StepCut
MaxStep…………
Initial; 脈動進行解を選択; 毎 step ごとにデータを出力する; 最大 step 回数; 分岐図を求めることを選択と変更する．

3. 実行$ cp␣../TE/final.dat␣Initial.dat





$ ./aout␣1 で実行する．


デモプログラム A1.2
Appendix1.2 の設定から実行までの操作は
$ cp␣../DEMO/Gray-Scott_Model/AS/11_AP1/*␣./ $ ./aout␣1 
で行える．

RD/DEMO/Gray-Scott_Model/AS/11_AP1 には Appendix1.2 の設定で必要なファイル．
• Function.c• Initial.dat• parameter.datが存在する．


4. 実行結果

99





…………
Step / MaxStep
iter / MaxIter
Control Parameter
PeriodVelocity
MachineEps
CC_Norm

=
=
=
=
=
=
=

1 / 1
0 / 303.380000000000000
423.333333333300004
0.001074709848399
0.000000010000000
0.010293873403380

Step / MaxStep=
1 / 1
iter / MaxIter=
1 / 30Control Parameter=
3.380000000000000
Period= 3495.413967540044723
Velocity=
0.001649339719272=
0.000000010000000
MachineEps
CC_Norm=
0.000000000141766--------------------------------------------------------------Convergence !!
--------------------------------------------------------------

「Convergence !!」と表示されているが，周期が明らかに大きすぎる値になっているので脈動進行解の枝に収束していないことがわかる．

Appendix1.3 脈動進行解の分岐計算
Appendix1.2 と同じ Initial.dat を初期値として再計算を行う．Appendix1.2 のようにうまく収束しなかった場合は，擬似弧長法を入れて試してみるとよい．
1. Initial.dat
Appendix1.2 と同じものを用いる．
2. parameter.dat


1



;

…………
KellerEqn…………

; 擬似弧長法の設定あり

と変更する．

3. 実行$ ./aout␣1 で実行する．

100

デモプログラム A1.3
Appendix1.3 の設定から実行までの操作は
$ cp␣../DEMO/Gray-Scott_Model/AS/12_AP2/*␣./ $ ./aout␣1 
で行える．

RD/DEMO/Gray-Scott_Model/AS/12_AP2 には Appendix1.2 の設定で必要なファイル．
• Function.c• Initial.dat• parameter.datが存在する．


4. 実行結果

101





…………

Step / MaxStep
iter / MaxIter
Control Parameter
PeriodVelocity
MachineEps
CC_Norm

=
=
=
=
=
=
=

1 / 1
2 / 303.372113877407729
423.332759369808628
0.001076007147629
0.000000010000000
0.000115980494426

Step / MaxStep
iter / MaxIter
Control Parameter
PeriodVelocity
MachineEps
CC_Norm

=
=
=
=
=
=
=

1 / 1
3 / 303.372366125499470
423.332638193281355
0.001076063757467
0.000000010000000
0.000000136682410

Step / MaxStep=
1 / 1
iter / MaxIter=
4 / 30Control Parameter=
3.372366652005462
Period=423.332638008723507
Velocity=
0.001076063893677
MachineEps=
0.000000010000000
CC_Norm=
0.000000000135279--------------------------------------------------------------Convergence !!
--------------------------------------------------------------



時間発展方程式の解を初期値として，脈動進行解の枝を求めることができた．脈動解の枝を求めるときも同様である．

102

第6章

デモ 2

この章では RD-AUTO を用いて 1 次元発熱反応モデルの分岐図を描く工程を説明する．分岐の種類についても調べる．
このデモに必要なデータは DEMO ディレクトリの中に図 6.0.1 のように入っている．
DEMO
Exothermic-reaction̲Model

AS

01̲SP̲INI

02̲SP̲01

…

TE
01̲SP

02̲TP

…

図 6.0.1: DEMO の中身

103

デモプログラム

内容

実行後の保存先

01_SP
02_TP
03_TB1
04_TB2
01_SPINI
02_SP01
03_SP02
04_SP03
05_SP04
06_SP05
07_SP06
08_SP0709_TP1_01
10_TP2_INI
11_TP2_01
12_TP2_02
13_TB3_01
14_TB3_02
15_TB3_03
16_TB3_04
17_TB3_05
18_TB_INI
19_TB_01
20_TB_02
21_TB_03
22_TB_04
23_TB_05

定常解 (TE)SP
TP

進行解 (TE)
脈動進行解 (TE)
脈動進行解 (TE)
定常解 (AS)
定常解 (AS)
定常解 (AS)
定常解 (AS)
定常解 (AS)
定常解 (AS)
定常解 (AS)
定常解 (AS)
進行解 (AS)
進行解 (AS)
進行解 (AS)
進行解 (AS)TB
SP_INI
SP_01
SP_02
SP_03
SP_04
SP_05
SP_06
SP_F
TP1_F
TP2_INI
TP2_F
TP2_F

進行解 (AS)
進行解 (AS)
進行解 (AS)
進行解 (AS)
進行解 (AS)
脈動進行解 (AS)
脈動進行解 (AS)
脈動進行解 (AS)
脈動進行解 (AS)
脈動進行解 (AS)
脈動進行解 (AS)

TP3_F
TB_INI
TB_01
TB_02
TB_03
TB_04
TB_05

104

時間発展方程式

6.1








































(

) )

(

1
u
4π 2
−au + expv ,ut = 2 du ∆u +Lε
1 + u/c
(
)
u
4π 2
v,
vt = 2 dv ∆v + h(1 − v) − expL
1 + u/cu(t, 0) = u(t, 2π),v(t, 0) = v(t, 2π),u(0, x) = u0 (x),v(0, x) = v0 (x),ux (t, 0) = ux (t, 2π), t > 0,vx (t, 0) = vx (t, 2π), t > 0,
0 ≤ x ≤ 2π,
0 ≤ x ≤ 2π

t > 0,
0 < x < 2π, 






(∗1)

}

(∗2)

}

(∗3)





この節では，1 次元発熱反応モデルに現れるパルス波を求める．カレントディレクトリは時間発
展方程式を解くところに移動 ($ cd␣RD/TE  ) して，最初に方程式と境界条件の設定を行う．

⃝
1 方程式 (∗1) の設定
方程式の拡散項と反応項の設定は Function.c で行う．拡散項は関数 Diffusion をvoid Diffusion( TEMP TP, double *par ){/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */par[TP.N_EP+1] = (4.0*M_PI*M_PI*par[1]) / (par[7]*par[7]);par[TP.N_EP+2] = (4.0*M_PI*M_PI*par[2]) / (par[7]*par[7]);/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */}と記述して，反応項は関数 Reaction をvoid Reaction( TEMP TP, double *par, double *x, double *F ){/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */int i;
int N = TP.N;
for( i = 1; i <= N; i++ ){F[i] = ((-par[4]*x[i])+exp(x[i]/(1.0+(x[i]/par[5])))*x[i+N])/par[3];F[i+N] = (par[6]*(1.0-x[i+N]))-exp(x[i]/(1.0+(x[i]+par[5]))*x[i+N]);}/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */}と記述する．

⃝
2 境界条件の設定
境界条件の設定は parameter.dat を開き，BC_Select を

105







3

;

…………
BC_Select
…………

;

OP[9]





と与える．BC_Select = 3 は周期境界条件を表す．

⃝
3⃝
1 〜⃝
2 以外の設定⃝
1 〜⃝
2 以外の設定はパルス波の種類によって異なるため，それらの設定はパルス波を求める
ところで逐次説明していく．
Step1 定常パルス波を時間発展方程式を用いて求める
ステップ関数を初期値として定常パルス波を求める．
1. 実行前の設定(a) 0parameter.dat
Function.c 内のパラメータの値を設定する．今回は以下のように設定する．


1.0
5.0
1.1
0.07
1.0
45.0
10

;
;
;
;
;
;
;



du ; par[1]dv ; par[2]eps ; par[3]a; par[4]c; par[5]h; par[6]Lx ; par[7](b) Function.c
反応項は
1f (u, v) = (−au + k(u)v)
ϵ(k(u) = exp (u))
1 + u/cg(u, v) = h(1 − v) − k(u)vであるので，関数 Reaction 内で

for( i = 1 ; i <= N ; i++ ){F[i]=((-par[4]*x[i])+exp(x[i]/(1.0+(x[i]/par[5])))*x[i+N])/par[3];F[i+N]=(par[6]*(1.0-x[i+N]))-exp(x[i]/(1.0+(x[i]+par[5]))*x[i+N]);
}



と記述．

(c) parameter.dat
今回は以下のように設定する．

106





256
; N
; 空間分割数は 256
128; N_Wave; 切断波数は 128-1; T_SELECT; (∗1); (∗1)-1; VEL_SELECT
1.0E-01; dt; 時間の刻み幅は 10−1
6.283185307179586; Lx; 区間は 2π100
; Gnuplot_CutN_Vel_Dis; (∗2)
0.2; Velocity_Site; (∗3)
10000; Velocity_Measure_Start ; step10000 か ら 速 度 を 測る．
3; BC_Select; 周期境界条件3; redEE_Select; (∗4)
500; StepCut; 500 回おきにデータを出力
100000; MaxStep; 最大 step 数; gnuplot 画面をすべて表示
0; Gnuplot_Select1; Initial_Select; 初期値はステップ関数(∗1) 発展方程式は ut = uxx + f (u), v t = v xx + g(v) を解く(∗2) 速度データを 100 回おきに出力する(∗3) u のパルス波の高さ 0.2 のところで高さを測る(∗4)Spectral 法 (Runge-Kutta 法)(d) boundary.dat
今回は周期境界条件のため，BC[1] ∼ BC[4] は特に変更する必要はない．
(e) initial_set.c
初期値の設定を行う．105 ページの (∗3) の初期値は図 6.1.1 のステップ関数を与
える．
3

u
v

2.5

2

1.5

1

0

1

2

3

図 6.1.1:

107

4

5

6

よって，プログラムでは次のように設定する．
void initial_set2( TEMP *TP, double *par, ……… ){…………
if( IS == 1 ){for( i = 1; i <= N; i++ ){x[i ] = 1.3992660000;
x[i+N] = 0.9378105400;
if( N/10 <= i && N/5>= i )x[i] = 5.0;
if( 1 <= i && N/10>= i )x[i+N] = 0.4;}}
else if( IS == 2 || IS == 3 ){
…………
}
…………
}



(f) gnuplot.c
gnuplot で描写する領域を指定する．図 6.1.1 の縦軸の範囲を −1.0 ∼ 20.0 で設定
する．
void Gnuplot_InitialCommand( TEMP *TP, int *c_G ){…………
/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/fprintf( GP2, "set yr[-1.0:20.0]\n" );fprintf( GP3, "set yr[-1.0:20.0]\n" );/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
…………
}



2. 実行


3. $ ./aout␣1 で実行する．

108













デモプログラム 1
Step1 の設定から実行までの操作は
$ cp␣../DEMO/Exothermic-reaction_Model/TE/01_SP/*␣./ $ ./aout␣1 
で行える．

RD/DEMO/Exothermic-reaction_Model/TE/01_SP には Step1 の設定で必要なファイル• 0parameter.dat• Function.c• parameter.dat• initial_set.c• gnuplot.cが存在する．


4. 実行結果
定常パルス波がでれば成功．

14

u
v

12
10
8
6
4
2
0
0

1

2

3

4

5

6

図 6.1.2: 定常パルス波

5. データの保存
ここで求めた定常パルス波のデータを SP というディレクトリの保存する．次のコマンドを実行することでデータをディレクトリ SP に保存できる．

109




$ ./aout␣SP

Step2 進行パルス波
ステップ関数を初期値として，進行パルス波を求める．
1. 実行前の設定(a) 0parameter.dat




…………
; dv ; par[2]
…………

3.0




と変更する．

(b) parameter.dat
変更点なし．
(c) initial_set.c
変更点なし．
2. 実行



3. $ ./aout␣1 で実行する．


デモプログラム 2
Step2 の設定から実行までの操作は$ cp␣../DEMO/Exothermic-reaction_Model/TE/02_TP/*␣./





$ ./aout␣1 
で行える．

RD/DEMO/Exothermic-reaction_Model/TE/02_TP には Step2 の設定で必要なファイル• 0parameter.dat• Function.c• parameter.dat• initial_set.c• gnuplot.cが存在する．


4. 結果
進行パルス波がでれば成功．

5. データの保存$ ./aout␣TP 
でディレクトリ TP にデータを保存する．

Step3.1 脈動進行パルス波
Step2 の結果を初期値として，脈動進行パルス波を求める．
110

18

u
v

16
14
12
10
8
6
4
2
0
0

1

2

3

4

5

6

図 6.1.3: 進行パルス波

1. 実行前の設定(a) 0parameter.dat


3.38



…………
; dv ; par[2]
…………





と変更する．

(b) parameter.dat




…………
3
; Initial_Select; (∗1)…………
100000; MaxStep; 最大 step 数…………
(∗1) 初期値 u0 (x) は Initial.dat に入っているものを用いる．パラメータは0parameter.dat のものを用いる．
と変更する．

(c) Initial.dat
Step2 で出力された final.dat を Initial.dat として使用する．
2. 実行$ cp␣TP/final.dat␣Initial.dat




111



$ ./aout␣1 で実行する．


デモプログラム 3.1
Step3.1 の設定から実行までの操作は
$ cp␣../DEMO/Exothermic-reaction_Model/TE/03_TB1/*␣./ $ ./aout␣1 
で行える．

RD/DEMO/Exothermic-reaction_Model/TE/03_TB1 には Step3.1 の設定で必要なファイル．
• 0parameter.dat• Function.c• parameter.dat• Initial.dat• gnuplot.cが存在する．


3. 結果
脈動進行パルス波が出れば成功．

14

u
v

12
10
8
6
4
2
0
0

1

2

3

4

図 6.1.4: 脈動進行パルス波

Step3.2 脈動進行パルス波
112

5

6

安定した脈動解を 100000step ほど計算する．

1. 実行前の設定(a) 0parameter.dat
変更点なし．
(b) parameter.dat
変更点なし．
(c) Initial.dat
Step3.1 で出力された final.dat を Initial.dat として使用する．
2. 実行$ cp␣final.dat␣Initial.dat





$ ./aout␣1 で実行する．


デモプログラム 3.2
Step3.2 の設定から実行までの操作は
$ cp␣../DEMO/Exothermic-reaction_Model/TE/04_TB2/*␣./ $ ./aout␣1 
で行える．

RD/DEMO/Exothermic-reaction_Model/TE/04_TB2 には Step3.2 の設定で必要なファイル．
• 0parameter.dat• Function.c• parameter.dat• Initial.dat• gnuplot.cが存在する．


3. データの保存

$ ./aout␣TB 
でディレクトリ TB にデータを保存する．

113

6.2

1 次元発熱反応モデルの分岐計算

ここでは，7.1 節で求めたパルス波の分岐図を求める．ただし，d を分岐パラメータとして与え
る．カレントディレクトリは分岐図を求めるところに移動 ($ cd␣RD/AS  ) して，最初に発熱反応モデル (105 ページの (∗1)) の拡散項と反応項の設定を行う．これらの設定は時間発展方程式と同じであり，拡散項は Function.c の関数 Diffusion をvoid Diffusion( TEMP TP, double *par ){/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */par[TP.N_EP+1] = (4.0*M_PI*M_PI*par[1]) / (par[7]*par[7]);par[TP.N_EP+2] = (4.0*M_PI*M_PI*par[2]) / (par[7]*par[7]);/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */}と記述して，反応項は Function.c の関数 Reaction をvoid Reaction( TEMP TP, double *par, double *x, double *F ){/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */int i;
int N = TP.N;
for( i = 1; i <= N; i++ ){F[i ] = ((-par[4]*x[i])+exp(x[i]/(1.0+(x[i]/par[5])))*x[i+N])/par[3];F[i+N] = (par[6]*(1.0-x[i+N]))-exp(x[i]/(1.0+(x[i]+par[5]))*x[i+N]);}/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */}と記述する．次に周期境界条件 (105 ページの (∗2)) の設定を行う．境界条件の設定は parameter.dat
を開き，BC_Select を




3

;

…………
BC_Select
…………

;

OP[9]





と与える．ただし，BC_Select = 3 は周期境界条件を表す．

Step1.1 定常解の分岐計算
7.1 節で求めた Step1 の解を用いて，定常解の分岐図を求める．ただし，Step1.1 では定常解の安定性を計算しないとする．定常解の選択およびその安定性計算による設定は parameter.dat
を開き，次のように設定する．

114





1

;

-1
-1

;
;

…………
SolSelect…………
EigenValue
EigenVector
…………

;

OP[8]

;
;

OP[11]
OP[12]



ただし，SolSelect = 1 は定常解を選択し，EigenValue と EigenVector を 1 以外の値に選択すると解の安定性を計算しない．

1. 実行前の設定(a) Initial.dat
7.1 節の Step1 で求めた解を Newton 法の初期値に与える．
$ cp␣../TE/SP/final.dat␣Initial.dat (b) parameter.dat

0.000000
0.000101-1200
1.0E-03
1
-1
2
1
3
6
-1
-1
1.0E-08
1
1
0
1



;
;
;
;
;
;
;
;
;
;
;
;
;
;
;
;
;
;

ds
PseudoNewton
KellerEqn
MaxIterdt
N_MSM
PeriodFix
ConPar
SolSelect
BC_Select
EE_Select
EigenValue
EigenVector
MachineEpsion
StepCut
MaxStep
Gnuplot
Initial2. 実行
上記の設定が終わったら，$ ./aout␣1 で実行する．

115

;
;
;
;
;
;
;
;
;
;
;
;
;
;
;
;
;
;

OP[1]
OP[2]
OP[3]
OP[4]
OP[5]
OP[6]
OP[7]
OP[7]
OP[8]
OP[9]
OP[10]
OP[11]
OP[12]
OP[13]
OP[14]
OP[15]
OP[16]
OP[17]



デモプログラム 1.1
Step1.1 の設定から実行までの操作は
$ cp␣../DEMO/Exothermic-reaction_Model/AS/01_SPINI/*␣./ $ ./aout␣1 
で行える．

RD/DEMO/Exothermic-reaction_Model/AS/01_SPINI には Step1.1 の設定で必要なファイル
• Function.c• Initial.dat• parameter.datが存在する．


3. 実行結果


…………

Step / MaxStep
iter / MaxIter
Control Parameter
MachineEps
CC_Norm

=
=
=
=
=

1 / 1
8 / 200
5.000000000000000
0.000000010000000
0.000270406062100

Step / MaxStep
iter / MaxIter
Control Parameter
MachineEps
CC_Norm

=
=
=
=
=

1 / 1
9 / 200
5.000000000000000
0.000000010000000
0.000000012730538

Step / MaxStep=
1 / 1
iter / MaxIter=
10 / 200
Control Parameter=
5.000000000000000
MachineEps=
0.000000010000000
CC_Norm=
0.000000000264208--------------------------------------------------------------Convergence !!
--------------------------------------------------------------



4. データの保存
Step1.1 で求めた定常解のデータを SP_INI というディレクトリに保存する．
次のコマンドを実行することでデータをディレクトリ SP_INI に保存できる．

116

$ ./aout␣SP_INI

Step1.2 定常解の分岐計算
次に，Step1.1 で求めた解を用いてパラメータを変化させながら定常解の枝を求める．
1. 実行前の設定(a) Initial.dat
初期値は Step1.1 で求めた解を用いる．
(b) parameter.dat




0.005

;

10
400

;
;

; パラメータの変化量 0.005

ds…………
StepCut
MaxStep…………

; 10step ごとにデータを出力する; 最大 step 回数と変更する．

2. 実行$ cp␣SP_INI/middle_final.dat␣Initial.dat





$ ./aout␣1 

デモプログラム 1.2
Step1.2 の設定から実行までの操作は
$ cp␣../DEMO/Exothermic-reaction_Model/AS/02_SP01/*␣./ $ ./aout␣1 
で行える．

RD/DEMO/Exothermic-reaction_Model/AS/02_SP01 には Step1.2 の設定で必要なファイル• Function.c• Initial.dat• parameter.datが存在する．


3. 実行結果

117





…………
Step / MaxStep
iter / MaxIter
Control Parameter
MachineEps
CC_Norm=====

400 / 400
3 / 200
6.999999999999957
0.000000010000000
0.000000123511199

Step / MaxStep
iter / MaxIter
Control Parameter
MachineEps
CC_Norm=====

400 / 400
4 / 200
6.999999999999957
0.000000010000000
0.000000015726073

Step / MaxStep=
400 / 400
iter / MaxIter=
5 / 200
Control Parameter=6.999999999999957
MachineEps=
0.000000010000000
CC_Norm=
0.000000001747026--------------------------------------------------------------Convergence !!
--------------------------------------------------------------…………


4. データの保存
step1.2 で求めた分岐図のデータを SP_01 のディレクトリに保存する．
$ ./aout␣SP_01 
Step1.3 定常解の分岐計算
Step1.2 で求めた枝の続きを求める．その際，ニュートン法では解が収束しなくなる．
実際，初期値に step1.2 で求めた final.dat を用いて実行すると，

118





…………
Step / MaxStep
iter / MaxIter
Control Parameter
MachineEps
CC_Norm=====

8 / 400
115 / 200
7.039999999999957
0.000000010000000
1252903.503148945979774

Step / MaxStep=
8 / 400
iter / MaxIter=
116 / 200
Control Parameter=7.039999999999957
MachineEps=
0.000000010000000
CC_Norm= 252355404632210472960.000000000000000--------------------------------------------------------------Not Convergence !!
-------------------------------------------------------------------------------------------------------------------------------CC_Norm is not Converge !! (2.523554046322105E+20)
CC_Norm is not Converge !! (2.523554046322105E+20)
CC_Norm is not Converge !! (2.523554046322105E+20)
CC_Norm is not Converge !! (2.523554046322105E+20)
CC_Norm is not Converge !! (2.523554046322105E+20)-----------------------------------------------------------------…………


となってしまう．これは saddle-node 分岐により解が存在しないところで計算してしまっているためである．なので、ニュートン法から疑似弧長法に切り替えて計算を行う．

1. 実行前の設定(a) Initial.dat
初期値は step1.2 で求めた final.dat を用いる．
(b) parameter.dat


1
4000



…………
; KellerEqn…………
; MaxStep…………

; 疑似弧長法; 最大ステップ数と変更する．

2. 実行$ cp␣SP_01/middle_final.dat␣Initial.dat $ ./aout␣1 
119





デモプログラム 1.3
Step1.3 の設定から実行までの操作は
$ cp␣../DEMO/Exothermic-reaction_Model/AS/03_SP02/*␣./ $ ./aout␣1 
で行える．

RD/DEMO/Exothermic-reaction_Model/AS/03_SP02 には Step1.3 の設定で必要なファイル• Function.c• Initial.dat• parameter.datが存在する．


3. 実行結果


…………

Step / MaxStep
iter / MaxIter
Control Parameter
MachineEps
CC_Norm=====4000 / 4000
1 / 200
2.230401984086159
0.000000010000000
0.000096142430039

Step / MaxStep
iter / MaxIter
Control Parameter
MachineEps
CC_Norm=====4000 / 4000
2 / 200
2.230402182048909
0.000000010000000
0.000000023149141

Step / MaxStep=
4000 / 4000
iter / MaxIter=
3 / 200
Control Parameter=
2.230402182084287
MachineEps=
0.000000010000000
CC_Norm=
0.000000000013074--------------------------------------------------------------Convergence !!
--------------------------------------------------------------…………


4. データの保存
step1.3 で求めた分岐図のデータを SP_02 のディレクトリに保存する．

120



$ ./aout␣SP_02 
Step1.4 定常解の分岐計算
parameter.dat の値を調節しながらニュートン法 (KellerEqn=-1) で計算を行い，ニュートン法で解が収束しない場合は疑似弧長法 (KellerEqn=1) で計算する，という作業を繰り返す．
そのため詳しい解説は省略し，デモと操作方法のみ記しておく．
step1.4 のデモと操作方法
1. デモの読み込み$ cp␣../DEMO/Exothermic-reaction_Model/AS/(*1)/* ␣./ 
2. 実行
$ ./aout␣1 
3. 保存
$ ./aout␣(*2) 
上記の 1\UTF{FF5E}3 の操作を 2 回繰り返す．(*1)，(*2) は，
1 回目 (*1):04_SP03(*2):SP_03
2 回目 (*1):05_SP04(*2):SP_04と変更して実行する．

Step1.2 から 1.4 で求めた分岐図を見るには，gnuplot を開き，
$plot "SP_01/81_CP_TV_NR.dat" u 1:4，"SP_02/81_CP_TV_NR.dat"u 1:4，

"SP_03/81_CP_TV_NR.dat" u 1:4，"SP_04/81_CP_TV_NR.dat" u 1:4 と実行する．分岐図は図 6.2.1 のようになり，saddle-node 分岐点が 2 つ存在する．
Step1.5 定常解の分岐計算
次に，図 6.2.1 の点 S から左側の分岐計算をする．
1. 実行前の設定(a) Initial.dat
初期値として SP_INI の middle_final.dat を用いる．
$ cp␣SP_INI/middle_final.dat␣Initial.dat (b) parameter.dat-0.005840

; パラメータの変化量 −0.005

;

ds
…………
; MaxStep…………

; 最大ステップ数と変更する．

2. 実行

$ ./aout␣1 で実行する．

121

u

∞

13

S

12
11
10
9

ステップの進む方向

87"SP_01/81_CP_TV_NR.dat" u 1:4"SP_02/81_CP_TV_NR.dat" u 1:4"SP_03/81_CP_TV_NR.dat" u 1:4"SP_04/81_CP_TV_NR.dat" u 1:4

6
5
4
1

2

3

4

5

6

7

8

9

10

d

図 6.2.1: 定常パルス波の分岐図

デモプログラム 1.5
Step1.5 の設定から実行までの操作は
$ cp␣../DEMO/Exothermic-reaction_Model/AS/06_SP05/*␣./ $ ./aout␣1 
で行える．

RD/DEMO/Exothermic-reaction_Model/AS/06_SP05 には Step1.5 の設定で必要なファイル．
• Function.c• Initial.dat• parameter.datが存在する．


3. 実行結果

122





…………

Step / MaxStep
iter / MaxIter
Control Parameter
MachineEps
CC_Norm=====

840 / 840
2 / 200
0.800000000000085
0.000000010000000
0.000974439979814

Step / MaxStep
iter / MaxIter
Control Parameter
MachineEps
CC_Norm=====

840 / 840
3 / 200
0.800000000000085
0.000000010000000
0.000000018493665

Step / MaxStep=
840 / 840
iter / MaxIter=
4 / 200
Control Parameter=
0.800000000000085
MachineEps=
0.000000010000000
CC_Norm=
0.000000000937121--------------------------------------------------------------Convergence !!
--------------------------------------------------------------…………


4. データの保存
step1.5 で求めた分岐図のデータを SP_05 というディレクトリに保存する．
$ ./aout␣SP_05



Step1.6 定常解の分岐計算
step1.5 で求めた枝の続きを求める．
1. 実行前の設定(a) Initial.dat
初期値として SP_05 の final.dat を用いる．
( $ cp␣final.dat␣Initial.dat  )(b) parameter.dat-0.005
1
3000



;

ds
…………
; KellerEqu…………
; MaxStep…………

; パラメータの変化量 −0.005; 疑似弧長法; 最大ステップ数

123

と変更する．

2. 実行

$ ./aout␣1 で実行する．


デモプログラム 1.6
Step1.6 の設定から実行までの操作は
$ cp␣../DEMO/Exothermic-reaction_Model/AS/07_SP06/*␣./ $ ./aout␣1 
で行える．

RD/DEMO/Exothermic-reaction_Model/AS/07_SP06 には Step1.6 の設定で必要なファイル．
• Function.c• Initial.dat• parameter.datが存在する．


3. 実行結果


…………
Step / MaxStep
iter / MaxIter
Control Parameter
MachineEps
CC_Norm=====

3000 / 3000
0 / 6000
9.999085259217248
0.000000010000000
0.005000000000000

Step / MaxStep
iter / MaxIter
Control Parameter
MachineEps
CC_Norm=====

3000 / 3000
1 / 6000
10.004076269062283
0.000000010000000
0.000001169675227

Step / MaxStep=
3000 / 3000
iter / MaxIter=
2 / 6000
Control Parameter=10.004076275823657
MachineEps=
0.000000010000000
CC_Norm=
0.000000000584759--------------------------------------------------------------Convergence !!
--------------------------------------------------------------…………


4. データの保存
124

step1.6 で求めた分岐図のデータを SP_06 というディレクトリに保存する．


$ ./aout␣SP_06

Step1.5 と 1.6 で求めた分岐図を見るには，gnuplot を開き，
$plot "SP_05/81_CP_TV_NR.dat" u 1:4，"SP_06/81_CP_TV_NR.dat"u 1:4 と実行する．分岐図は図 6.2.2 のようになる．

u

∞

13

S

12
11

"06_SP05/81_CP_TV_NR.dat" u 1:4"07_SP06/81_CP_TV_NR.dat" u 1:4

10
9
8
7

ステップの進む方向

6
5
4
3
0

1

2

3

4

5

6

7

8

9

10

11

d

図 6.2.2: 定常パルス波の分岐図

Step1.7 定常解の分岐計算
step1.2 から step1.6 までのデータを結合する．その際，図 6.2.1 と図 6.2.2 をみるとわかるよ
うに分岐計算をすすめた向きがそれぞれ異なる．なので，ここでは step1.5 と step1.6 で求めたデータの向きを逆にして結合する．

1. step1.5 と step1.6 のデータを結合する．
次のコマンドで，step1.5 と step1.6 の middle_final.dat のデータを結合し，middle_final.datに保存する．
$cat␣SP_05/middle_final.dat␣SP_06/middle_final.dat > middle_final.dat 
2. 結合したデータのステップ数を逆にする．

125



$./Sort (詳しい説明はデモ 1 の分岐計算 step1.5 を参照)
3. step1.2 から step1.4 までのデータとステップ数を逆にしたデータを結合し，middle_final.datに保存する．
$cat␣middle_final.dat␣SP_01/middle_final.dat␣SP_02/middle_final.dat␣
SP_03/middle_final.dat␣SP_04/middle_final.dat > tmp.dat $mv␣tmp.dat␣middle_final.dat 
4. 結合したデータのステップ数を修正する
parameter.dat を次のように変更する．




9

…………
; Initial;OP[17]



$ ./aout␣1 で実行する．
実行後，結合した middle_final.dat に対応する• 81_CP_TV_NR.dat• 82_S_CP_TV_NR_REV_IMV_ABV_NUEV.dat• solution_valu.dat• 1• 11が出力される．ただし，1 はステップ数が修正された middle_finel.dat，11 はステッ
プ数が修正された 9middle_final.dat である．

5. データの出力
分岐図のデータを gnuplot で出力する (図 6.2.3)．
$ gnuplot $ p "81_CP_TV_NR.dat"␣u␣1:4␣w␣l 

6. データの保存
1 と 11 をそれぞれ middle_final.dat， 9middle_final.dat に上書きする．
$ cp␣1␣middle_final.dat $ cp␣11␣9middle_final.dat 
まとめたデータを SP_F に保存する．
$ ./aout␣SP_F 

Step1.8 定常解の分岐計算
ここでは Step1.7 で求めた定常解の安定性を求める方法について説明する．解の安定性は，⃝
1 解の枝を求めながら安定性を求める方法と
126

u

∞

14

定常解
12

10

8

6

4
0

1

2

3

4

5

6

7

8

d

図 6.2.3: 定常解の分岐図⃝
2 解の枝を求め終わった後、それらの解の安定性を調べる方法の 2 通りある．ここでは，定常解の枝を Step1.7 で既に求めたので 2 の方法を用いる

1. 実行前の設定(a) Initial.dat
Step1.7 で保存したディレクトリ SP_F の中から安定性計算に必要とする 9middle_final.dat
を Initial.dat に変更して計算するディレクトリに持ってくる．
$ cp␣SP_F/9middle_final.dat␣Initial.dat(b) parameter.dat




…………
8
; Initial; (∗1)(∗1) 各 step における解の固有値と固有関数の計算を行う
と変更する．

2. 実行

$ ./aout␣1  で実行する．

127





デモプログラム 1.8
Step1.8 の設定から実行までの操作は
$ cp␣../DEMO/Exothermic-reaction_Model/AS/08_SP07/*␣./ $ ./aout␣1 
で行える．

RD/DEMO/Exothermic-reaction_Model/AS/08_SP07 には Step1.8 の設定で必要なファイル．
• Function.c• Initial.dat• parameter.datが存在する．


3. 実行結果安定性の結果は 82_S_CP_TV_NR_REV_IMV_ABV_NUEV.dat のデータを参照する．gnuplot を開き，
$ gnuplot > plot␣"82_S_CP_TV_NR_REV_IMV_ABV_NUEV.dat"␣u␣1:8と実行することで横軸にステップ数，縦軸に固有値実部を与えたグラフが表示される(図 6.2.4)．
固有値実部

400
"82_S_CP_TV_NR_REV_IMV_ABV_NUEV.dat" u 1:8

350
300
250
200
150
100
1

50

2

4

5

0
3

-50

6

7

8

-100
0

200

400

600

800

1000

1200

Step数

図 6.2.4:

step の変化と共に固有値の実部が 0 を横切るとき，分岐が起きる．今回の場合，8 つの分岐が存在する．
ここでひとつ問題があり，今回扱っているモデルでは，定常解の固有値を安定に求める128

ことが困難である．なので，ここでは分岐の種類についての詳しい説明は省略する．
図 6.2.4 の 1〜8 では次の 3 つの種類の分岐が現れる．

• saddle-node 分岐 (1,5,8)• Hopf 分岐 (4,6)• pitch-fork 分岐 (2,3,7)
4. データの保存
Step1.8 で求めた分岐図のデータを SP というディレクトリに上書きする．
$ ./aout␣SP_F

このとき，ディレクトリを上書きしますかと尋ねられるので，$ yes と入力する．

Step2.1 進行解の数値計算
ここでは Step1.8 で求めた定常解の pitch-fork 分岐点における固有関数の情報を用いて進行解を求める．用いる固有関数の情報は，固有値の実部が正である step のものを用いる．Step1.8でも示したように，pitch-fork 分岐点での固有値の情報は以下のようになっている．


左から 1 列目の step 数，8 列目の固有値実部，9 列目の固有値虚部の値．

274
274
274
274
274
275
275
275
275
275

…………
………………… 2.90284773E+02………………… 8.78498815E-02………………… -4.49039817E-01………………… -1.46363103E+00………………… -7.30440658E+01………………… 2.89238921E+02………………… 1.34657221E-01………………… 9.15488035E-03………………… -1.03611995E+00………………… -7.30439543E+01…………

0.00000000E+00
0.00000000E+00
0.00000000E+00
0.00000000E+00
0.00000000E+00
0.00000000E+00
0.00000000E+00
0.00000000E+00
0.00000000E+00
0.00000000E+00

…………
…………
…………
…………
…………
…………
…………
…………
…………
…………


1. 実行前の設定(a) Initial.dat
pitch-folk 分岐点の固有関数を Initial.dat として使用する．固有ベクトルのデー
タは 6_EigenVector_Data に入っている．
6_EigenVector_Data/ 275_CP . . . 1 /o1 . . . 2 のデータを Initial.dat に変更して持ってくる．
(b) parameter.dat

1
2

ここにはコントロールパラメータの値が入っている
ここには固有値実部の値が入っている

129





; ds
; パラメータの変化量なし; KellerEqn; 擬似弧長法の設定なし…………
; 進行解を選択
2; SolSelect…………
155; MaxStep; 最大 step 回数…………
2; Initial; (∗1)(∗1) 分岐点から伸びる解の枝を求める-0.005000-1と変更する．

(c) 0_BranchChangeParameter.dat
分岐点における零固有値に対する固有関数の情報と 0_BranchChangeParameter.datのパラメータを用いて Newton 法に用いる初期値 (uini , cini , pini ) を構成する.
• 初期値の構成 uini = u0 + εΦR ,• 速度の構成 cini = δc ,• パラメータの構成 pini = δp ,ただし，u0 は定常解の pitch-fork 分岐点における解 u，ΦR は分岐点における零固有値に対する固有関数の実部，ε, δc , δp はそれぞれ 0_BranchChangeParameter.dat
における eps , velocity_delta , par_delta である.
ここでは，






0.05
0.00
0.0
0.00000000
0.000

;
;
;
;
;

eps
period_delta
velocity_delta
par_deltatimeと変更する．

2. 実行

$ ./aout␣1  で実行する

130; uini の構成に必要; ここでは不必要な設定; 速度の微小変量; パラメータの微小変量; ここでは不必要な設定







デモプログラム 2.1
Step2.1 の設定から実行までの操作は
$ cp␣../DEMO/Exothermic-reaction_Model/AS/09_TP1_01/*␣./ $ ./aout␣1 
で行える．

RD/DEMO/Exothermic-reaction_Model/AS/09_TP1_01 には Step2.1 の設定で必要なファイル．
• Function.c• Initial.dat• parameter.dat• 0_BranchChangeParameterが存在する．


3. 実行結果


…………

Step / MaxStep
iter / MaxIter
Control Parameter
Velocity
MachineEps
CC_Norm

=
=
=
=
=
=

155 / 155
2 / 200
0.003779154067403
5.773071012535138
0.000000010000000
0.000000037012342

Step / MaxStep=
155 / 155
iter / MaxIter=
3 / 200
Control Parameter=
0.003779154067403
Velocity=
5.773071012718584
MachineEps=
0.000000010000000
CC_Norm=
0.000000000004661--------------------------------------------------------------Convergence !!
--------------------------------------------------------------



収束した解は速度を持っているため，進行解であるとわかる．
進行解の分岐図を表示すると図 6.2.5 のようになる．
図 6.2.5 を定常解とあわせて表示すると次のようになる．

4. データの保存
Step2.2 で求めた分岐図のデータを TP1_F というディレクトリに保存する．
$ ./aout␣TP1_F




131

u

∞

7.5

"81_CP_TV_NR.dat" u 1:4

7

進行解の分岐図

6.5

6

5.5

5
0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

図 6.2.5: 進行解の分岐図

Step3.1 進行解の分岐計算
7.1 節で求めた Step2 の解を用いて，進行解の分岐図を求める．
1. 実行前の設定(a) Initial.dat
7.1 節の Step で求めた解を Newton 法の初期値に与える．
$ cp␣../TE/TP/final.dat␣Initial.dat (b) parameter.dat

132

0.8

d

u
14

∞

12

10

定常解の分岐図
8

進行解の分岐図
6

4
0

1

2

3

4

5

6

7

8

d

図 6.2.6: 進行解の分岐図


0.000000
0.000101-130
1.0E-03
1
-1
2
2
3
6
-1
-1
1.0E-08
1
1
0
1

;
;
;
;
;
;
;
;
;
;
;
;
;
;
;
;
;
;

ds
PseudoNewton
KellerEqn
MaxIterdt
N_MSM
PeriodFix
ConPar
SolSelect
BC_Select
EE_Select
EigenValue
EigenVector
MachineEpsion
StepCut
MaxStep
Gnuplot
Initial2. 実行
上記の設定が終わったら，$ ./aout␣1 で実行する．

133

;
;
;
;
;
;
;
;
;
;
;
;
;
;
;
;
;
;

OP[1]
OP[2]
OP[3]
OP[4]
OP[5]
OP[6]
OP[7]
OP[7]
OP[8]
OP[9]
OP[10]
OP[11]
OP[12]
OP[13]
OP[14]
OP[15]
OP[16]
OP[17]



デモプログラム 3.1
Step3.1 の設定から実行までの操作は
$ cp␣../DEMO/Exothermic-reaction_Model/AS/10_TP2_INI/*␣./ $ ./aout␣1 
で行える．

RD/DEMO/Exothermic-reaction_Model/AS/10_TP2_INI には Step3.1 の設定で必要なファイル• Function.c• Initial.dat• parameter.datが存在する．


3. 実行結果

134





…………

Step / MaxStep
iter / MaxIter
Control Parameter
Velocity
MachineEps
CC_Norm

=
=
=
=
=
=

1 / 1
3 / 302.000000000000000
12.928757323597416
0.000000010000000
0.032783248608212

Step / MaxStep
iter / MaxIter
Control Parameter
Velocity
MachineEps
CC_Norm

=
=
=
=
=
=

1 / 1
4 / 302.000000000000000
12.928384371030262
0.000000010000000
0.000002199030632

Step / MaxStep=
1 / 1
iter / MaxIter=
5 / 30Control Parameter=
2.000000000000000
Velocity=
12.928384353630330
MachineEps=
0.000000010000000
CC_Norm=
0.000000000231466--------------------------------------------------------------Convergence !!
--------------------------------------------------------------…………


4. データの保存
Step3.1 で求めた定常解のデータを TP2_INI というディレクトリに保存する．
次のコマンドを実行することでデータをディレクトリ TP2_INI に保存できる．

$ ./aout␣TP2_INI

Step3.2 進行解の分岐計算
次に，Step3.1 で求めた解を用いてパラメータを変化させながら進行解の枝を求める．
1. 実行前の設定(a) Initial.dat
初期値は Step3.1 で求めた解を用いる．
(b) parameter.dat

135





0.002

;

1
330

;
;

; パラメータの変化量 0.005

ds…………
StepCut
MaxStep…………

; 10step ごとにデータを出力する; 最大 step 回数と変更する．

2. 実行$ cp␣TP2_INI/final.dat␣Initial.dat





$ ./aout␣1 

デモプログラム 3.2
Step3.2 の設定から実行までの操作は
$ cp␣../DEMO/Exothermic-reaction_Model/AS/11_TP2_01/*␣./ $ ./aout␣1 
で行える．

RD/DEMO/Exothermic-reaction_Model/AS/11_TP2_01 には Step3.2 の設定で必要なファイル
• Function.c• Initial.dat• parameter.datが存在する．


3. 実行結果

136





…………

Step / MaxStep
iter / MaxIter
Control Parameter
Velocity
MachineEps
CC_Norm

=
=
=
=
=
=

330 / 330
6 / 302.659999999999927-0.090468446226469
0.000000010000000
0.000003118381549

Step / MaxStep=
330 / 330
iter / MaxIter=
7 / 30Control Parameter=
2.659999999999927
Velocity=-0.090468666535193
MachineEps=
0.000000010000000=
0.000000001455116
CC_Norm--------------------------------------------------------------Convergence !!
--------------------------------------------------------------…………


進行解の分岐図を表示すると図 6.2.7 のようになる．

u

∞

17
進行解

16.5
16
15.5
15
14.5
14
13.5
13
12.5
12
2

2.1

2.2

2.3

2.4

2.5

2.6

図 6.2.7: 進行解の分岐図
図 6.2.7 を定常解とあわせて表示すると図 6.2.8 のようになる．

4. データの保存

137

2.7

d

u
14

∞

12

進行解の分岐図

108

定常解の分岐図

6

4
0

1

2

3

4

5

6

7

8

d

図 6.2.8: 進行解の分岐図

step3.2 で求めた分岐図のデータを TP2_F のディレクトリに保存する．
$ ./aout␣TP2_F 
Step3.3 進行解の分岐計算
ここでは Step3.2 で求めた進行解の安定性を求める方法について説明する．
1. 実行前の設定(a) Initial.dat
まず，Step3.2 で保存したディレクトリ TP2_01 の中から固有値計算に必要なデー
タ 9middle_final.dat を Initial.dat に変更して計算するディレクトリに持ってくる．
$ cp␣TP2_01/9middle_final.dat␣Initial.dat(b) parameter.dat


8





;

…………
Initial; 安定性の計算を選択と変更する．

2. 実行

$ ./aout␣1  で実行する．

138







デモプログラム 3.3
Step3.3 の設定から実行までの操作は
$ cp␣../DEMO/Exothermic-reaction_Model/AS/12_TP2_02/*␣./ $ ./aout␣1 
で行える．

RD/DEMO/Exothermic-reaction_Model/AS/12_TP2_02 には Step3.3 の設定で必要なファイル
• Function.c• Initial.dat• parameter.datが存在する．


3. 実行結果
82_S_CP_TV_NR_REV_IMV_ABV_NUEV.dat のデータから横軸にステップ数，縦軸に固有値実部をグラフで示すと図 6.2.9 となる．
固有値実部

40
"82_S_CP_TV_NR_REV_IMV_ABV_NUEV.dat" u 1:8
Hopf分岐点

20

0

-20

-40

-60

Step数

-80
0

50

100

150

200

250

300

350

図 6.2.9:
図 6.2.9 をみると，step300 付近で固有値実部が 0 を横切っているのがみてとれる．そこで，82_S_CP_TV_NR_REV_IMV_ABV_NUEV.dat を開いてみてみると，308step で固有値の実部が 2 つ 0 を横切っている．つまりここで Hopf 分岐が起こっていることがわかる．Hopf 分岐点からは脈動進行解の枝が存在する．

139





左から 1 列目の step 数，8 列目の固有値実部，9 列目の固有値虚部の値．

307 …………………
307 …………………
307 …………………
307 …………………
307 …………………
308 …………………
308 …………………
308 …………………
308 …………………
308 …………………

…………
1.61302794E+00 0.00000000E+00-8.27802402E-01 -1.44360682E+02-8.27802402E-01 1.44360682E+02-4.29835612E+00 0.00000000E+00-7.15261082E+01 0.00000000E+00
1.60629017E+00 0.00000000E+00
1.07230184E-02 -1.43964738E+02
1.07230184E-02 1.43964738E+02-4.15896757E+00 0.00000000E+00-7.15701018E+01 0.00000000E+00
…………

…………
…………
…………
…………
…………
…………
…………
…………
…………
…………





4. データの保存
Step3.3 で求めた分岐図のデータを TP2_F というディレクトリに上書きする．
$ ./aout␣TP2_F



Step4 進行解の数値計算
Step4 では Step1.8 で求めた定常解の pitch-fork 分岐点から分岐する進行解の分岐図を求める．その際，分岐点における固有関数の情報を用いるという今までのやり方ではここで求めたい枝を求めるのは困難である．そこで，このでは違う方法を用いて進行解の枝を求めていく．
Step4.1 進行解の数値計算
ここではまず，コントロールパラメーターを ϵ に変更して，ϵ ＝ 0.001025 のときの定常解の分岐図を求める．
1. 実行前の設定(a) Initial.dat
初期値は SP_INI にある middle_final.dat のデータを用いる．
$cp␣SP_INI/middle_final.dat␣Initial.dat (b) parameter.dat

140





-0.000005
0.000101-130
1.0E-04
1
-1
3
1
3
6
-1
-1
1.0E-08
1
35
0
1

;

ds
; OP[1]; PseudoNewton; OP[2]; KellerEqn; OP[3]; MaxIter; OP[4]; dt; OP[5]; N_MSM; OP[6]; PeriodFix; OP[7]; ConPar; OP[7]; SolSelect; OP[8]; BC_Select; OP[9]; EE_Select; OP[10]; EigenValue; OP[11]; EigenVector; OP[12]; MachineEpsion ; OP[13]; StepCut; OP[14]; MaxStep; OP[15]; Gnuplot; OP[16]; Initial; OP[17]と変更する．

2. 実行

$ ./aout␣1 で実行する．


デモプログラム 4.1
Step4.1 の設定から実行までの操作は
$ cp␣../DEMO/Exothermic-reaction_Model/AS/13_TP3_01/*␣./ $ ./aout␣1 
で行える．

RD/DEMO/Exothermic-reaction_Model/AS/13_TP3_01 には Step4.1 の設定で必要なファイル．
• Function.c• Initial.dat• parameter.datが存在する．


Step4.2 進行解の数値計算
Step4.1 で求めた解を用いて定常解の分岐図を求める．
1. 実行前の設定
141

(a) Initial.dat
先ほど求めた解を初期値として用いる．
$cp␣final.dat␣Initial.dat (b) parameter.dat

0.0050000

;

ds
…………
; ConPar…………
; EigenValue; EigenCector…………
; StepCut; MaxStep…………
;
Initial

2
1
1
1
200
1



; 固有値を計算する; 固有ベクトルを計算する

と変更する．Step4.2 では固有値の計算をしながら定常解の枝を追っていく．

2. 実行



$ ./aout␣1 

デモプログラム 4.2
Step4.2 の設定から実行までの操作は
$ cp␣../DEMO/Exothermic-reaction_Model/AS/14_TP3_02/*␣./ $ ./aout␣1 
で行える．

RD/DEMO/Exothermic-reaction_Model/AS/14_TP3_02 には Step4.2 の設定で必要なファイル．
• Function.c• Initial.dat• parameter.datが存在する．


82_S_CP_TV_NR_REV_IMV_ABV_NUEV.dat を開いてみてみると，step166 で固有値実部の値が 0 を横切っている．詳しい説明は省略するが，ここでは pitch-fork 分岐により，進行解が分岐する．

142





165
165
165
165
165
166
166
166
166
166

…………
…………………
8.54593778E-02 0.00000000E+00…………………
-2.39380306E+00 0.00000000E+00…………………
-3.76071247E+01 1.66351451E+02…………………
-3.76071247E+01 -1.66351451E+02…………………
-4.31666931E+01 2.24716640E+01…………………
-3.89203220E-02 0.00000000E+00…………………
-2.25037219E+00 0.00000000E+00…………………
-3.74751869E+01 1.66530530E+02…………………
-3.74751869E+01 -1.66530530E+02…………………
-4.29838236E+01 2.18762624E+01
…………

…………
…………
…………
…………
…………
…………
…………
…………
…………
…………





Step4.3 進行解の数値計算
ここでは，Step4.2 で求めた pitch-fork 分岐点の固有関数の情報を用いて，進行解を求める．
1. 実行前の設定(a) Initial.dat
step165 では正の値だった 1 行目の固有値実部が，step166 では負の値になっている．したがって，初期値として step165 の固有関数の情報を用いる．
$cp␣6_EigenVector_Data/165_CP.../o1...␣Initial.dat (b) parameter.dat


0.0050000
-1
2
-1
-1
200
2





;
;

ds
KellerEqn…………
; SolSelect…………
; EigenValue; EigenCector…………
; MaxStep…………
; Initialと変更する．

2. 実行

$ ./aout␣1 
で実行する

143

; パラメータの変化量なし; 擬似弧長法の設定なし; 進行解を選択

デモプログラム 4.3
Step4.3 の設定から実行までの操作は
$ cp␣../DEMO/Exothermic-reaction_Model/AS/15_TP3_03/*␣./ $ ./aout␣1 
で行える．

RD/DEMO/Exothermic-reaction_Model/AS/15_TP3_03 には Step4.3 の設定で必要なファイル．
• Function.c• Initial.dat• parameter.datが存在する．


Step4.4 進行解の数値計算
ここでは Step4.1 と同様の操作で，ϵ ＝ 0.001025 のときの進行解の分岐図を求める．
1. 実行前の設定(a) Initial.dat
初期値は Step4.3 で求めた final.dat のデータを用いる．
$cp␣final.dat␣Initial.dat (b) parameter.dat

0.000005

;

3

;

35

;

1

;



ds
…………
ConPar…………
MaxStep…………
Initialと変更する．

2. 実行

$ ./aout␣1 で実行する．

144

デモプログラム 4.4
Step4.4 の設定から実行までの操作は
$ cp␣../DEMO/Exothermic-reaction_Model/AS/16_TP3_04/*␣./ $ ./aout␣1 
で行える．

RD/DEMO/Exothermic-reaction_Model/AS/16_TP3_04 には Step4.4 の設定で必要なファイル．
• Function.c• Initial.dat• parameter.datが存在する．


Step4.5 進行解の数値計算
Step4.1〜Step4.4 までの操作により，ϵ=0.0012 のときの進行解の枝を見つけることができた．ここでは Step4.4 で求めた解を用いて進行解の分岐図を求める．
1. 実行前の設定(a) Initial.dat
初期値は Step4.4 で求めた final.dat のデータを用いる．
$cp␣final.dat␣Initial.dat (b) parameter.dat-0.005
-1
2
114
1





;
;

ds
KellerEqn…………
; ConPar…………
; MaxStep…………
; Initialと変更する．

2. 実行

$ ./aout␣1  で実行する

145







デモプログラム 4.5
Step4.5 の設定から実行までの操作は
$ cp␣../DEMO/Exothermic-reaction_Model/AS/17_TP3_05/*␣./ $ ./aout␣1 
で行える．

RD/DEMO/Exothermic-reaction_Model/AS/17_TP3_05 には Step4.5 の設定で必要なファイル．
• Function.c• Initial.dat• parameter.datが存在する．
求めた進行解の分岐図を定常解のものとあわせて表示すると次のようになる (※はデモ
では求めていない部分の進行解の枝)．

u

∞

14

12

※
10

進行解の分岐図
8

定常解の分岐図
6

d

4
0

1

2

3

4

5

6

7

8

図 6.2.10: 定常解と進行解の分岐図

3. データの保存
Step4.5 で求めた分岐図のデータを TP3_F というディレクトリに保存する．
$ ./aout␣TP3_F



Step5.1 脈動進行解の分岐計算
ここでは Step3.3 で求めた進行解の Hopf 分岐点の固有関数の情報を用いて脈動進行解を求める．
1. TP2_F の固有関数の情報
146

TP2_F で求めた固有関数の情報は保存されていないので，求めなおす必要がある．
次のコマンドで保存した TP2_F のデータをもってくる．
$ cp␣TP2_F/*␣./ 
実行する．
$ ./aout␣1 

2. 実行前の設定(a) Initial.dat
Hopf 分岐点の固有関数を Initial.dat として使用する．固有ベクトルのデータ
は 6_EigenVector_Data に入っている．6_EigenVector_Data/ 308_CP . . . 3 /o1. . . 4 のデータを Initial.dat に変更して持ってくる．
$cp␣6_EigenVector_Data/308_CP.../o1...␣Initial.dat (b) parameter.dat
ds; パラメータの変化量なし…………
4; SolSelect; 脈動進行解を選択…………
1; StepCut; 毎 step ごとにデータを出力する
1; MaxStep; 最大 step 回数…………
2; Initial; (∗1)(∗1) 分岐点から伸びる解の枝を求める．
0.000000;と変更する．

(c) 0_BranchChangeParameter.dat
分岐点における零固有値に対する固有関数の情報と 0_BranchChangeParameter.datのパラメータを用いて Newton 法に用いる初期値 (uini , Tini , cini , pini ) を構成する.
• 初期値の構成 uini = u0 + ε(sin(2πt)ΦR + cos(2πt)ΦI ),2π• 周期の構成 Tini =+ δTβ• 速度の構成 cini = δc ,• パラメータの構成 pini = δp ,ただし，u0 は定常解の pitch-fork 分岐点における解 u，ΦR , ΦI はそれぞれ分岐点における零固有値に対する固有関数の実部と虚部，c は解の速度，T は解の周期 p はコン
トロールパラメータ，c0 は進行解の Hopf 分岐点における速度，ε, δT , δc , δp , t はそれ
ぞれ 0_BranchChangeParameter.dat における eps , period_delta , velocity_delta, par_delta , time である.
ここでは，

3
4

ここにはコントロールパラメータの値が入っている
ここには固有値実部の値が入っている

147





0.05
0.00
0.0
0.00000000
0.000

;
;
;
;
;

eps
period_delta
velocity_delta
par_delta
time

;
;
;
;
;

SS_2
SS_2
SS_2

&
&
&
&
&

SS_3
SS_3
SS_3
SS_3

&
&
&
&
&

SS_4
SS_4
SS_4
SS_4
SS_4





と変更する．

3. 実行$ ./aout␣1 
で実行する




デモプログラム 5.1
Step5.1 の設定から実行までの操作は
$ cp␣../DEMO/Exothermic-reaction_Model/AS/18_TB_INI/*␣./ $ ./aout␣1 
で行える．

RD/DEMO/Exothermic-reaction_Model/AS/18_TB_INI には Step5.1 の設定で必要なファイル．
• Function.c• Initial.dat• parameter.dat• 0_BranchChangeParameterが存在する．


4. 実行結果実行すると，Hopf 分岐点から振動進行解の枝を求めるための初期値が得ら
れる．

148





…………
Step / MaxStep
iter / MaxIter
Control Parameter
PeriodVelocity
MachineEps
CC_Norm

=
=
=
=
=
=
=

1 / 1
8 / 302.641842624371586
0.043900163639821
4.887823615525028
0.000000010000000
0.000002305873995

Step / MaxStep=
1 / 1
iter / MaxIter=
9 / 30Control Parameter=
2.641842624544382
Period
=0.043900163647529
Velocity=
4.887823617448827=
0.000000010000000
MachineEps
CC_Norm=
0.000000000632522--------------------------------------------------------------Convergence !!
--------------------------------------------------------------…………


ただし，プログラムのバグで「Convergence !!」と表示されていても，周期 (Period) が負になっていたり，極端に大きな値になっていたりする場合は，求めたい解に収束していないことに注意する．

5. データの保存
step5.1 で求めた分岐図のデータを TB2_INI のディレクトリに保存する．
$ ./aout␣TB2_INI 
Step5.2 進行解の分岐計算
Step5.2 では，Step5.1 で求めた解を用いて，振動進行解の枝を求める．
実行と保存を繰り返すため，詳しい説明は省略してデモプログラムを記しておく．
step5.2 のデモと操作方法
1. デモの読み込み$ cp␣../DEMO/Exothermic-reaction_Model/AS/(*1)/* ␣./ 
2. 実行
$ ./aout␣1 
3. 保存
$ ./aout␣(*2) 
上記の 1\UTF{FF5E}3 の操作を 2 回繰り返す．(*1)，(*2) は，
1 回目 (*1):19_TB_01(*2):TB_01
2 回目 (*1):20_TB_02(*2):TB_02

1493 回目 (*1):21_TB_03
4 回目 (*1):22_TB_04
5 回目 (*1):23_TB_05と変更して実行する．

(*2):TB_03(*2):TB_04(*2):TB_05

分岐図を出力するには，それまでに保存したデータを結合する．
$ cat␣TB_01/81_CP_TV_NR.dat␣TB_02/81_CP_TV_NR.dat␣TB_03/81_CP_TV_NR.dat
␣TB_04/81_CP_TV_NR.dat␣TB_05/81_CP_TV_NR.dat ＞ tmp.dat $ mv␣tmp.dat␣81_CP_TV_NR.dat 
次のコマンドで定常解と Step3 で求めた進行解の分岐図と一緒に分岐図を出力できる．
$ gnuplot $ p "SP_F/81_CP_TV_NR.dat"␣u␣1:4␣w␣l, "TP2_F/81_CP_TV_NR.dat"␣u␣1:4␣w␣l
, "81_CP_TV_NR.dat"␣u␣1:4␣w␣l 

u
14

∞

12

10
定常解
進行解

8

脈動進行解

6

4
0

1

2

3

4

5

6

7

8

d

図 6.2.11:

Step5.2 で求めたデータを TB_F というディレクトリに保存する．
$ ./aout␣TB_F 
Step6 脈動定常解の分岐計算
Step6 では，定常解の Hopf 分岐点より分岐する振動定常解の分岐図を求める．
しかし，脈動定常解の枝を追うのは非常に難しいため，このマニュアルでは取り扱わない．
分岐図は以下のようになる．

150

u

∞

14

12

10

脈動定常解の分岐図
8

定常解の分岐図
6

4
0

1

2

3

4

5

6

図 6.2.12:

Step7 分岐図
今まで求めた分岐図をまとめて出力すると図 6.2.13 のようになる．

151

7

8

d

u

∞

14

進行解2(Step3)

12
進行解3(Step4)
脈動定常解

10

脈動進行解

8

定常解

進行解1(Step2)

6

4
0

1

2

3

4

5

図 6.2.13: 分岐図

152

6

7

8

d

第7章

付録

鳥瞰図

7.1

この節では，第 6 章で求めた解の鳥瞰図を描く工程を説明する．
この節では，GLSC がインストールされていることを前提としている．
ここでは例として，RD/DEMO/Allee-EFFECT_Model/TE/04_bird のデータを用いて鳥瞰図を描く．

Step1 必要なデータを作成する


時間発展方程式を解くディレクトリに移動 ($ cd␣RD/TE  ) する．
RD/DEMO/Allee-EFFECT_Model/TE/04_bird のデータを計算するディレクトリに持ってくる．
$ cp␣../DEMO/Allee-Effect_Model/TE/04_bird/*␣./ 
1. MaxStep と StepCut を決める．
MaxStep と StepCut は好きなように指定してよいが，MaxStep/StepCut の値によって作成される鳥瞰図のメッシュが異なる．
MaxStep/StepCut が大きくなると，線が密になる．図 7.1.1
MaxStep/StepCut が小さくなると，線が疎になる．図??

200

t

0
0

x

2π

図 7.1.1: MaxStep/StepCut = 200
ここでは，MaxStep/StepCut = 200 となるように以下のように設定する．




…………
10 ; StepCut2000 ; MaxStep…………
2; Initial_Select;;OP[12]
OP[13];

OP[15]


153

2. 実行$ ./aout 
3. 実行結果
図 7.1.2 のような振動進行パルスが出る．
1

0.8

0.6

0.4

0.2

0

0

1

2

3

4

5

6

図 7.1.2: 振動進行パルス

実行すると，鳥瞰図を描くために必要なデータ (solution_value.dat) が作られる．

Step2 鳥瞰図を描く
鳥瞰図は RD/TE で描く．

1. glsc の設定
glsc を次のように設定する．




…………
#----------------------------------------------------------------$Initial_Select= 2 ; #(*1)$OutPut_Cut_Data= 1 ; #この数字ごとのデータが出力される．
$U_SolutionValue_Point = 3 ; #solution_value.dat の 3 列目を u とする．
$V_SolutionValue_Point = 4 ; #solution_value.dat の 4 列目を v とする．
#----------------------------------------------------------------$Dx_Zbottom_U= 0.0;#出力するときの u の 0 の位置.
$Dx_Ztop_U= 0.0;#出力するときの u の高さに関係する．
…………
$Graph_X= 30.0;#出力するときの横幅．
$Graph_Y= 20.0;#出力するときの縦幅．
…………




(*1) Initial_Select の値によって出力されるものが違う.
154

• InitialSelect = 1 の時，u と v のそれぞれの鳥瞰図が並んで出力される.
• InitialSelect = 2 の時，u の鳥瞰図が出力される.
• InitialSelect = 3 の時，v の鳥瞰図が出力される.
• InitialSelect = 4 の時，最初に u と v のそれぞれの鳥瞰図が並んで出力され，画像をクリックするたびに u の鳥瞰図，v の鳥瞰図，と出力されていく．
2. 実行
次のコマンドで鳥瞰図が作成される．
$ ./glsc 

3. 実行結果
図 7.1.3 のような鳥瞰図が出力される．また，画像を消すことでディレクトリ TE/DATA が作成され、画像 TE/DATA/U.i00.epsi が作られる．
Initial_Select = 1 の時は，U_V.i00.epsi が作られる．
Initial_Select = 3 の時は，V.i00.epsi が作られる．
Initial_Select = 4 の時は，U.i00.epsi，V.i00.epsi，U_V.i00.epsi が作られる．
次のコマンドで画像が作られていることを確認できる．
$ display␣DATA/U.i00.epsi 
200

t

0
0

x

2π

図 7.1.3: 鳥瞰図 p80

4. 視点を変える
次のコマンドで出力される画像の視点が変わる．
$ ./glsc␣pXX 

ただし，XX には 00 から 90 までの数字が入り，デフォルト ($ ./glsc  ) では XX = 80 として出力される．図 7.1.3
XX = 90 では真上からの視点になる．図 7.1.4

155

XX = 00 では真横からの視点になる．図 7.1.7

200

200

t

t

0

0
0

x

2π

x

0

図 7.1.4: 鳥瞰図 p902π

図 7.1.5: 鳥瞰図 p60

200

t

0

ｘ

0
0

x

2π

図 7.1.7: 鳥瞰図 p00

図 7.1.6: 鳥瞰図 p30

156

2π

7.2

実行時に起きるエラーの対処法$./aout 1 実行時に以下のようなエラーが起きた場合に対する対処法を記載する．

図 7.2.1: エラー時のキャプチャ
図 7.2.1 のようなエラーが発生した場合，以下のコマンドを使うことで解決できる可能性が
ある．







$ yum search glibc-static$ su 
# yum install glibc-static.x86_64


157




