第13回(2011年7月13日)


今回は拡散方程式(熱伝導方程式)と、簡単な反応拡散方程式を取り扱います。次回、波動方程式の取り扱いへと移っていく予定です.

お知らせ:次週(7/20) の演習では、期末テストを予定しています。その他に期末レポートも出題の予定です。

期末テストのお知らせ


演習内容 


1変数反応拡散系(前回紹介した内容です)

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

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

 

という形で与えられます(プラス, 適当な境界条件). 例えば ui(x, t) を位置 x 時刻 t における物質 i の濃度と考えると、Fi( { ui(x, t) } ) がその場所、時刻での i の(化学反応等に依る)濃度の変化規則を表しており、後の項が、i の拡散を意味しています( Di は 物質 i の拡散係数.)。

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

 

この方程式は、先の拡散方程式に -au(x,t) という項が付け加わったものとなっています。これは、例えば熱伝導をイメージするならば、系の両端以外に熱が(寒い)外に逃げている効果を表し、また物質の拡散をイメージするなら、物が途中で漏れている(分解されている)効果を表しています。



では、このような「漏れ」のある場合の様子を, L = 1, a = 1, D = 0.1, g(x) = 1 として見てみましょう.


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

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」, βは「beta」等とする.)

xj , uj (j = 0, 1, 2, ... N, N+1) は,配列変数 x[j] , u[j] 等として与える。今回は j = 0 ~ N+1 に対する合計 N+2 個の x[j] , u[j] が必要なので, この変数を定義(宣言)する部分では, x[N+2] , u[N+2] 等とする.

(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 の順に


(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 ) + beta( uj-1 - 2uj + uj+1 )

を繰り返す

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

uj = new_uj

を繰り返す

u0 = 2A - u1 , uN+1 = 2B - uN   (境界条件)


を繰り返す

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


課題6

L = 1, a = 1, D = 0.1, A = 2, B = 0.2 , g(x) = 1 として 上記の方程式を数値的に解け。また a の値をいろいろ変え、その様子を考察せよ ( a を0から大きくして行きつつ、課題4の結果との比較をせよ.)。

実行サンプル (mpg)


課題7

A = 1 + sin(wt), B = 1 とした場合の様子を、様々な w, a, D について数値的に解き、その様子を考察せよ。( a を変えつつ課題5の結果との比較をせよ.)。

実行サンプル (mpg)




2変数反応拡散系

 拡散方程式は熱や物質が拡散していく様子を表現しており、それが示す通り、物質や熱は一様にならされる方向に時間変化していきます。しかし例えば2種類以上の物質が反応しながら拡散する場合、そのような一様化が阻害され、むしろ不均一性が増加する場合がありえます。このとき現れる化学成分の不均一な分布パターンは、チューリングパターン(Turing pattern)と呼ばれます. 例えば次のブラッセレーターモデル は、チューリングパターンを生成する典型的な方程式系となっています. (このような現象の解析は、後期の「現象数理学」等で扱われる予定です。)

このモデルでは、時刻 t, 位置 x (0 < x < π) における2種類の化学物質 u, v の密度 u(x,t), v(x,t) が、以下の方程式

に従って反応と拡散をすると考えます。

ここで初期条件を

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


この条件で a, b, Du, Dv を適度に選べば、u, v は時間、空間的に均一には成らず、大きな不均一性を持った状態が実現します。
今回は例えば a = 2.0, b = 4.0, Dv = 0.1 等として, Du を適度に与えてみることで、この様子を考察してみてください.

アルゴリズムは以下のようになっています。

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

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

h, τはそれぞれ x, t に関する格子点の間隔で,h = L/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 の順に


(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   (境界条件)


を繰り返す


課題8

a = 2.0, b = 4.0, Dv = 0.1 (もしくはもっと小さく) 等として, Du を適度に与えてみることで、空間的に不均一な u , v の分布kが形成される様子を考察してみてください, また a, b, Du, Dv を変化させた場合も考察してください

実行サンプル(Dv = 0.01 としてます) (mpg)



期末レポート

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

課題は以下の通り。

(1) 反応拡散系、(課題8で扱っている)チューリングパターン (またはアラン・チューリング本人) について調べ, A4 用紙 2〜4枚程度に、手書きでまとめよ. (数学事務室のレポートボックスに提出すること)
 (自分で資料を探し、頭の中で何か考えてまとめる事. 参考文献を記す事. 必ず手書きの事.)

(2) 以下で紹介する波動方程式をに関し、以下の課題(9, 10, 11, 12) の何れかを解け。作成したプログラムはファイル名をrep0727.c として 下の提出用ページより提出する事. プログラムの先頭にコメント文で名前、学籍番号と選んだ課題番号. またプログラムの説明もコメント文として適宜記入する事.

提出用ページ

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




波動方程式の数値解法(7/13)

今回は1次元波動方程式(wave equation)を扱います。ここでは,1次元波動方程式をコンピュータを用いて解くために差分化を行い,コンピュータを用いて波動方程式を解くアルゴリズムを示します.今日の内容に関する詳しい情報は講義で説明があります(ました).演習では細かい事は気にせずに一気にアルゴリズムの紹介まで行きたいと思います.

さて,天下り的ですが1次元波動方程式は、下図のような弦や(深さのある水の)水面の時空間的な振動を表現する方程式系です(導出等は講義で説明されます. ).下記方程式 (W) における u(x,t) は弦(もしくは水面)の変位をあらわします.演習では1次元波動方程式を数値計算を用いて解き,その解 u(x,t) を画面上にアニメーションとして表示する事を目標とします.つまり,コンピュータ内部に弦の振動を再現し,それをアニメーションとして表示するわけです.




(W) 

時間ー空間2階偏微分方程式の差分化

さて,次の1次元波動方程式の初期値境界値問題、

(W) 

但し、c >0, f(0)=0, f(L )=0, g(0)=0, g(L )=0.

を考えます。境界条件,初期条件が適切に与えられれば(W)の解が弦の振動を表現していると言えるでしょう.これまで扱ってきた拡散方程式との違いは、
i) u(x,t) の時間 ( t ) に関する2階微分が含まれている事,
ii) それに伴って, 初期条件に u(x,t) の時間の1階微分に関する条件がつく事,
です

さて,数値計算する為にまず(W)を離散化する必要があります。拡散方程式のときにも紹介しましたが、分割方法は以下の2種類が考えられますが,

ここでは時間方向の離散化は「京都方式」を,空間方向の離散化は「札幌方式」を採用します.(拡散方程式と同様です.)

[0,L ] 区間を N 等分するとし,上記離散化方法に従い, (δx = ) h = L /N, τ >0 とし,それを用いて xj = ( j - 1/2)h, tk = として,

と書くことにします。(W)の第一式は、左辺・右辺それぞれ u(x,t) の t に関する2階微分、x に関する2階微分を含んでいます。 これらは拡散方程式における, 空間に関する2階微分の差分化(以下のスライドを参照)と同様、次のように離散化されます。

よって λ=/h として、まとめると

という式が得られます。


境界条件

拡散方程式と同様に、空間に関しても端を適度に扱う必要があり, そのために境界条件が与えられます. 境界の処理の仕方も拡散方程式と同様です,

上記の Dirichlet 条件は、両端を u(x,t) = 0 に固定した場合の物と成っています、拡散方程式の時と同様に
  すなわち 

  すなわち 

で与えると、両端の u の値を好きな値 A , B で固定した境界条件を与えられます。


初期条件

拡散方程式同様, 初期条件を空間全域に対して与える必要があります. ここで拡散方程式と波動方程式の大きな違いが現れます。今回の波動方程式を差分化した上記の式では、 時刻 tk-1 , tk での U ( つまり Uik-1 , Uik ) から次の時刻 tk+1 での U ( つまり Uik+1 ) が計算できるという形をしています。したがって計算を開始するためには、初期値として、 が必要となります。

まず は(W)の第三式より、

となります。 を決める方法はいくつか考えられますが、例えば次のようにします。まず、テイラーの公式を用いると、

となります。(W)より

であり,さらに t=0 でも(W)の第一式が成り立つとすると、

となります。この式の右辺は2階差分商を用いて

と近似できます。つまり,

と求まります.以下ではこのように求まった を用います。


解くべき(差分)方程式

以上の考察から解くべき式をまとめると、以下のようになります.

(2)

但し、 ( λ=/h ) が安定性のための必要十分条件(安定性条件とは,数値計算がうまく行く為の条件です.なぜにこういった条件が必要であるかは,講義で説明があります).

めでたく差分方程式が得られました。(2)をコンピューターを用いて解くのが今週の目標です。


アルゴリズム

1次元波動方程式を解く為のアルゴリズムは次のようになります。安定性条件を満たすため λ ≦ 1 となるように h, τ, を与えます.

アルゴリズム中の 「 xj 」, 「 uj 」, 「 old_uj 」, 「 new_uj 」 はそれぞれ実数型の配列として定義すると良い.境界条件処理用に配列要素はN個よりも余分に2つ必要となり,計 N+2 個必要となる事に注意(new_uj については,実際には境界処理用の要素は用いない).例えば uj について,C 言語風に配列を定義すると,

float u[N+2];

となる.このようにすれば,u[0], u[1], ... , u[N+1] が用意される.実際そのように用意すれば良い.

アルゴリズム

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

T(時間上限), M(時間刻み数), L(空間幅), N(空間刻み数), c を設定する

h = L/N, τ=T/M, λ = τ*c/h (λ ≦ 1 となるように. M, N を調整, プログラム中では, λは「ramda」, τは「tau」などとするとよい.)

(1.5) GLSC の準備

(2) 初期値設定

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

xj = j*h - 0.5*h

を繰り返す ( f(x) は課題9以降を参照 )

(境界条件)


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

を繰り返す ( g(x) は課題9以降を参照 )

(境界条件)


(3)計算, 動画描画

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

(x[1], u[1]) 〜 (x[N], u[N]) を画面に表示 (g_move + g_plot, x[0], u[0], x[N+1], u[N+1]

は境界処理用なので表示しない)

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

を繰り返す


( uj , old_uj 両方の値を更新します。以下の順序を守らないと計算がうまく行きません。(なぜだか分かりますか?))

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

を繰り返す

(境界条件)


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

を繰り返す

(境界条件)


を繰り返す


課題9

L=1,  c=1,  f(x)=2*x*(1-x),  g(x)=0,  u(0,t)=u(1,t)=0  として (2) を解くプログラムを作成してください。(N は100程度, T ~ 100 で λ << 1 となるように M を適度に選び, 描写頻度も適度に選ぶ事.)

実行サンプル (mpg)


課題10

L=1, c=1, u(0,t)=u(1,t)=0 (N = 100 程度) はそのままで、f(x) および g(x) を他のもにかえてみよ。例えば関数 g(x)=0、関数 f(x) を

f(x)=exp(-100*(x-0.4)*(x-0.4)) ... (i)

としたり

f(x)=exp(-100*(x-0.4)*(x-0.4)) - exp(-50*(x-0.7)*(x-0.7)) + exp(-10*(x-0.5)*(x-0.5)) ... (ii)

とする等。
(局在化した初期値からはじめると、波の伝播が観察される。)

実行サンプル (i)  実行サンプル (ii) (mpg)


課題11

他の境界条件,(i) Neumann 境界条件もしくは (ii) 周期境界条件、を適応したプログラムを作成し,実際に数値計算してみよ.

実行サンプル (i) 実行サンプル (ii) (mpg)


課題12

境界条件として (端に適当な値を与える) ディレクレ条件で, A = sin(wt), B = 0 とした場合について、様々な c, w について数値的に解き ( 初期条件 : f(x)=0 , g(x)=0 )、、c, w を適度に選んで、(バネ振動のときに見たような)「共鳴(共振)」を観察せよ. (以前の「バネ振動子」に比べ、より「吊り橋の共振」に近い状況がシミュレーションできます.)

実行サンプル(非共鳴的) 実行サンプル(共鳴的) (mpg)


おまけ

境界条件として Neumann 条件を用いて,初期条件として
g(x)=0、
f(x)
  = 1 ( 0.4 ≦ x ≦ 0.6 ) ,
  = 0 ( Otherwise )
等として計算しその様子を観察せよ.

実行サンプル (N = 5000 としています. ) (mpg)


先日大きな被害を出した津波は、地震によって生じた海底の隆起(沈降)によって生じる急激な海面の変形(今回の初期条件のような形状)が波として伝搬したものです。その様子は波動方程式で観察されます。



浅い所では波と海底が擦れて海底付近の海水の流れが遅くなるため、波の上部が多いかぶさるようになります.このような現象は、今回扱っている波動方程式の適用外になっています。


戻る