第14, 15回(2009年7月22, 29日)


前回と今回の2回にわたって反応拡散方程式を扱います。反応拡散系については、後期の「現象数理学」で扱われるかと思われます。(昨年までは。。。)


演習内容 


1変数反応拡散系

空間のある局所的な部分における何らかの量(例えば物質の濃度の変化、温度の変化、磁力の強さ、魚や人の密度、、、等等、定量化できるものなら何でも)の変化と、それが空間全体に広がっていく様子を記述する場合、よく「反応拡散系」と呼ばれるクラスの微分方程式系が用いられます。もともとは近くにいる分子同士が化学反応しつつ拡散する様を記述していたので、反応拡散系と呼ばれています。

1次元反応拡散方程式は、一般的に変数の組 { ui(x, t) } に対して

 

という形で与えられます(プラス, 適当な境界条件). 例えば ui(x, t) を位置 x 時刻 t における物質 i の濃度と考えると、Fi( { ui(x, t) } ) がその場所、時刻での i の(化学反応等に依る)濃度の変化規則を表しており、後の項が、i の拡散を意味しています( Di は 物質 i の拡散係数.)。期末レポートで扱うブラッセレーター系(このページの下の方参照)もこの方程式系の1種であることが分かると思います。

今回はまず例題として、次の1変数反応拡散系について、差分化とアルゴリズムを示します。

 

(ここで u = u(x,t) です )この方程式は「ギンツブルグ-ランダウ方程式」と呼ばれるもので、歴史的には超伝導の理論等で活躍したり、磁石のモデル方程式として考えられてきています。多数決のような事をする方程式だと思えば良いかもしれません。(実際我々の周りにある磁石の磁性(向き)は、それを構成する分子の持つ小さな磁気が、多数決で向きを合わせる感じに決まります。この場合 u(x,t) は、プラス or マイナスの値をとりますが、これが 局所的な磁気の向き(プラス方向、マイナス方向)を表します。)

L = 1, a = 1, b = 1, D = 0.1 として, 境界条件をノイマン条件, 初期条件を P を 0 〜 L の定数として

  g(x) = -1 ... (x > P)
  g(x) = 1 ...  (x ≦ P)

とした場合の様子を見てみましょう.

拡散方程式に対して行ったのと同様に計算機で扱えるよう離散化し,数値的に近似解を求めます.

h, τはそれぞれ x, t に関する格子点の間隔で,h = 1/Nとします.また,xj = (j - 1/2)h, tk = とします(分割の仕方は前回を参照).さらに,

と書く事にしましょう.

拡散方程式と同様にして,

  

は次の差分方程式で近似されます.

これを書き換えれば,

     

が得られます.ただし,α = Dτ/h2 としました.この式は見てのとおり,時刻 tk での u から次の時刻 tk+1 での u の値が計算できる形をしています.ようは,これをプログラムにすればよいわけです.後は,初期条件と境界条件(今回はノイマン条件)ですが,それらはそれぞれ次のようになります.

    

これで,準備はできました.

次にアルゴリズムを考えましょう.とはいえもう分かると思いますが、拡散方程式の場合に少し手を加えるだけです。解を求める時刻の上限を T とし,適当な自然数 M を定めて,τ=T/M とします.

アルゴリズム

(1)初期パラメータ設定

T(時間上限), M(時間刻み数), L(空間幅), N(空間刻み数), D(拡散係数), a, b を設定する

h = L/N, τ = T/M, α = Dτ/h2 (プログラム中では, τは「tau」, αは「alpha」等とする.)

(1.5) GLSC の準備

(2) 初期値設定

j := 1,....,N の順に

xj = (j - 0.5)h

uj = g( xj )

を繰り返す

u0 = u1, uN+1 = uN   (境界条件)

(3)

k := 0,1,.....,M-1 の順に


(3.1)(x1, u1) 〜 (xN, uN) を画面に表示

g_move 関数で曲線の始点 (x1, u1)与える.

j := 1,2,....,N の順に

g_plot 関数で (xj, uj) を線でつなぐ.

(ただし x0, u0, xN+1, uN+1 は境界処理用なので表示しない)


(3.2)計算

j := 1,2,....,N の順に

new_uj = uj + tau( auj - buj3 ) + α( uj-1 - 2uj + uj+1 )

を繰り返す

j := 1,2,....,N の順に

uj = new_uj

を繰り返す

u0 = u1, uN+1 = uN   (境界条件)


を繰り返す

このアルゴリズムを実現するC のプログラムを作成してみてください.


課題1

L = 1, a = 1, b = 1, D = 0.1 として, 境界条件をノイマン条件, 初期条件を P を 0 〜 L の定数として

  g(x) = -1 ... (x > P)
  g(x) = 1 ...  (x ≦ P)

とした場合で、p=0.45, 0.5, 0.55 等とした場合の様子を見てください。またD = 0 の場合(拡散が無く、各場所が孤立している場合)と見比べてください。(T = 10, M = 50000, N = 50 程度で.)

D = 0 つまり拡散が無い場合, u(x,t) は初期条件に依存して +1 or -1 に近づきます。しかし D が有限で拡散がある場合, 互いに値を「揃えよう」とするので、綱引きのような事が起こります.

課題2

D = 0.01, T = 50, g(x) = sin(20x)+cos(10x)+0.1  等とした場合の様子を見てください。

また g(x) = sin(20x)+cos(10x)-0.1 とした場合はどうでしょう。


このプログラムを2変数( U(x,t), V(x,t) )にそのまま拡張すれば, 期末レポートの(2)は出来ると思います。




期末レポート

期末レポートを提出してもらいます。提出〆切は8/4の予定です。

課題は以下の通り。

(1) 反応拡散系, チューリングパターンに関することを, A4 2枚程度にまとめよ.(自分の頭の中で何か考えてまとめる事. 必ず手書きの事.)


(2) チューリングパターンを生成する偏微分方程式系を数値的に解き, チューリングパターンを生成するプログラムを作成せよ

 例えば次のブラッセレーターモデル を解くプログラムを作成し、チューリングパターンが生成される様子をみてください.

時刻 t, 位置 x (0 < x < π) における化学物質 u, v の密度を u(x,t), v(x,t) とすると それらが

に従うとする.

ここで初期条件を

とし、境界条件はノイマン条件とする.


この条件で a, b, Du, Dv が適度に選ばれていれば、チューリングパターンが生成されるはずです。 (この適度なパラメーターのセットを見つけてください. 課題(1) をやっているうちに、パラメータの条件が分かると思われます. )

参考までに ... 差分化の際 T(時間上限)= 100, M(時間刻み数)= 100000, L(空間幅)= π, N(空間刻み数)= 100 程度と設定する

ノイマン条件を課すので、ディレクリ条件で

u0 = -u1, uN+1 = -uN

等としていたところを,

u0 = u1, uN+1 = uN

 等と変更する


提出方法:

(1) は数学事務室のレポートボックスに提出、(2) は作成したプログラムをrep0728.c として report フォルダに置く事. (8/5に自動回収します.)プログラムの先頭にコメント文で名前と学籍番号. またプログラムの説明もコメント文として適宜記入する事.

〆切はどちらも8/4中とします。

   成績は [これまでのレポートの点数]×[期末レポートの点数] でつけます. (つまり未提出では単位はありませんのでご注意ください.)




ブラッセレーターモデルの差分化とアルゴリズム

1変数系に対して行ったのと同様に計算機で扱えるよう離散化し,数値的に近似解を求めます.

h, τはそれぞれ x, t に関する格子点の間隔で,h = 1/Nとします.また,xj = (j - 1/2)h, tk = とします(分割の仕方はこの回を参照).さらに,

と書く事にしましょう.

1変数系と同様にして,

は次の差分方程式で近似されます.

これを書き換えれば,

     

が得られます.ただし,αU = Du×τ/h2, αV = Dv×τ/h2 としました.後は,初期条件と境界条件(今回はノイマン条件)ですが,それらはそれぞれ次のようになります.

    

これで,準備はできました.

次にアルゴリズムを考えましょう.とはいえもう分かると思いますが、1変数系の場合に少し手を加えるだけです。解を求める時刻の上限を T とし,適当な自然数 M を定めて,τ=T/M とします.

アルゴリズム

(1)初期パラメータ設定

T(時間上限), M(時間刻み数), L(空間幅), N(空間刻み数), Du(U の拡散係数), Dv(V の拡散係数), a, b を設定する

h = L/N, τ = T/M, αU = Du×τ/h2, αV = Dv×τ/h2 (プログラム中では, τは「tau」, αUは「alpha_u」等とする.)

(1.5) GLSC の準備

(2) 初期値設定

j := 1,....,N の順に

xj = (j - 0.5)h

uj = sin( xj ) + sin( 7xj )

vj = sin( xj ) + sin( 7xj )

を繰り返す

u0 = u1, uN+1 = uN

v0 = v1, vN+1 = vN   (境界条件)

(3)

k := 0,1,.....,M-1 の順に


(3.1)(x1, u1) 〜 (xN, uN) を画面に表示

g_move 関数で曲線の始点 (x1, u1)与える.

j := 1,2,....,N の順に

g_plot 関数で (xj, uj) を線でつなぐ.

(ただし x0, u0, xN+1, uN+1 は境界処理用なので表示しない)

を繰り返す

(3.2)(x1, v1) 〜 (xN, vN) を (色を変えてから) 画面に表示

g_move 関数で曲線の始点 (x1, v1)与える.

j := 1,2,....,N の順に

g_plot 関数で (xj, vj) を線でつなぐ.

(ただし x0, v0, xN+1, vN+1 は境界処理用なので表示しない)

を繰り返す


(3.3)計算

j := 1,2,....,N の順に

new_uj = uj + tau( a + (uj)2 vj - (b+1)uj ) + αU( uj-1 - 2uj + uj+1 )

new_vj = vj + tau( -(uj)2 vj + buj ) + αV( vj-1 - 2vj + vj+1 )

を繰り返す

j := 1,2,....,N の順に

uj = new_uj

vj = new_vj

を繰り返す

u0 = u1, uN+1 = uN

v0 = v1, vN+1 = vN   (境界条件)


を繰り返す

このアルゴリズムを実現するC のプログラムを作成してみてください.


課題3

a, b, Du, Dv を変化させ、チューリングパターン以外の解の様子を調べよ.

課題4

(1) 次の方程式系

     

において, a = 0.025, b = 0.0542, Du = 0.00001, Dv = 0.00002 等とし、解の様子を調べよ. またそれらの値 (特に a: 0.02 〜 0.03 辺りで) や初期条件を変化させたときの解の様子を調べよ. (T = 20000, L = 2.0, M = 1000000, N = 200 程度、境界はノイマン条件。)

この方程式は Gray-Scott(グレイ-スコット)モデルと言われるもので、自己増殖パルス(スポット)解という特徴的な解を持つ事で知られています。(上のパラメータで現れます。)またパラメータ次第で、非常にリッチに解が変化します。

(2) Gray-Scott モデルについても調べてみてください。

課題5

ローレンツ方程式で各変数が拡散するとしたときの反応拡散方程式を書き下し、その解の様子を調べよ。(ローレンツ方程式は単体でカオスになりますが、拡散を考える事で乱流を記述するモデルになります.)

おまけ(交通渋滞のモデル...*先週と同じです...)

(余裕のある方、トライしてみてください.)

(1) 一車線道路上の車の動きがシンプルなルールに従うとして、交通渋滞の様子をモデル化してみます。簡単のため、空間は j = 1, 2, ... N の格子点に区切られていて、各格子点上には車が1台だけ居られるとします。時刻も k = 0, 1, 2, ... と離散的な量として扱います。車は前方( 格子点番号 j が大きくなる方向 )に進むとします.

ある時刻 m で n 番目の格子に車が居る場合、その車の動きは以下の規則に従うとします。

 n+1 番目の格子(前方の格子点)に別の車が居なければ、次の時刻(m+1)に n+1 番目の格子に移動する。 (前が空いていたら進む。)

 n+1 番目の格子(前方の格子点)に別の車が居る場合、次の時刻(m+1)も n 番目の格子にとどまる。 (前が空い無ければ進まない。)

 


ここで時刻 k, 場所 j に車が存在したら Ukj = 1, いなければ Ukj = 0 とすると、Ukj の変化は以下のような差分方程式系で記述できることが分かります。

   Ujk+1 = F( Ukj-1 ,  Ukj ,  Ukj+1 )

       F(0, 0, 0) = 0, F(0, 0, 1) = 0, F(0, 1, 0) = 0, F(0, 1, 1) = 1,
       F(1, 0, 0) = 1, F(1, 0, 1) = 1, F(1, 1, 0) = 0, F(1, 1, 1) = 1
 

この方程式系について、境界条件を周期境界条件: Uk0 = UkN, UkN+1 = Uk1 とし、初期条件を

U0j = 1 ( sin( 0.5*j ) + P > 0 の場合 )

U0j = 0 ( それ以外 )

として. その様子を図示してみてください. (P : -1 〜 1 で大きいほど車の台数が増えます。)

表示方法は, 縦軸 k, 横軸 j, とし, Ukj = 1 のとき, g_circle, g_marker 等で ( j, k ) に丸(点)を描く等とすること.
( M = 50, N =50 程度でこんな絵が描けます.(点は半径 0.25 の円と長さ 0.5 の線. 時間お方向は下→上、車は左→右に進む。)

 

 P = -0.95(車少ない)各車はスムーズに前方に進む

 

 P = -0.5(車少し多い)各車は(まだ)スムーズに前方に進む (1格子点分の車間距離は空いているので、止まらない。)

 

 P = 0.5(車多い)渋滞ができて進めない車がいる。渋滞自身は後方に移動していく(車は前方にすすむ。)。

このモデルは、非常に単純化されていますが現実の交通流のもつ特徴を良く捉えています。この方法をもっと複雑な道路網やネットワークに拡張して、都市交通やインターネットのパケットの渋滞等の解析も行われています。



戻る