第12回(2011年7月6日)


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

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


演習内容 


前回の内容(復習)


1次元拡散方程式の数値解法

1次元拡散方程式(1次元熱伝導方程式とも呼ばれる)の初期値境界値問題を数値的に解きます.それは,次のようなものです.(記号など多少異なるかもしれません.)

問題 0.1

 

この問題の近似解を数値計算にて求めるのが今回の目標です.まずは,上問題の特別な場合として次の問題を考えましょう.なお,式中の D は拡散係数または熱伝導率とよばれます.

問題 0.1.1


この方程式は、拡散方程式もしくは熱伝導方程式と呼ばれます。拡散とは物が散らばっていく様子(例:コップにインクを垂らすと、時間とともに広がって最後には一様になる)であり、熱伝導とは熱が伝わる様子(例:スプーンを熱湯に浸けていると徐々に熱が伝わって持っている部分まで熱くなる)であります。これらの現象を表現しているのが上式なのです。物の散らばりと熱の伝わりという異なる現象が同じ方程式で表されるという面白さがあります。熱伝導方程式だとみるとu(x,t)は時刻tでの温度分布、拡散方程式だと思うとu(x,t)は粒子の密度分布に対応しています。

ここでは, x = 0 と x = Lでいつも u = 0 となる条件を境界(空間の端)に課しています。温度の分布の言葉でいうと、両端を氷か何かでつねに冷やしている状況を考えていることになります。つまり、熱が両端からどんどんと奪われているわけです。


偏微分方程式の差分化

さてこの拡散方程式(問題0.1、問題0.0.1)の数値的に近似解を求めます.数値計算をするためには、微分を差分で置き換える必要があります。これまで演習で扱ってきた常微分方程式系では、t (時間) に関する微分のみを扱ってきましたが、拡散方程式には x (空間)に関する微分も含まれているため、こちらに対しても適当な差分化を行う必要があります。 差分化の方法については前回紹介しましたのでそちらを参考にしてください。 そのおさらいも兼ねて、問題0.1.1を離散化して解いてみましょう.

[0,L ] 区間を N 等分するとし,h = L /N とします. τ >0 ,xj = ( j - 1/2)h ,  tk = として,

と書く事にしましょう.すると拡散方程式

  (1)

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

これを書き換えれば,

      (2)

が得られます.ただし,α = Dτ/h2 としました. 問題0.1.1 における境界条件は Dirichlet 条件ですので.

すなわち,

となります.

よって解くべき式が

   

      (2)
   

     (3)
   

    (4)

と得られ、これで準備ができました.

次に拡散方程式を計算機にて解くためのアルゴリズムを考えましょう.解を求める時刻の上限を T とし,適当な自然数 M を定めて,τ=T/M とします.問題0.1.1を解くアルゴリズムは次のようになります. (D=1.0 とします).


アルゴリズム

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

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

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

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 = 2xj(1-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-1 + (1 - 2α)uj + αuj+1

を繰り返す

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

uj = new_uj

を繰り返す

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


を繰り返す

このアルゴリズムを実現するC のプログラムを作成せよというのが当面の課題となります.


課題1

上のアルゴリズムに従ってプログラムを完成させよ.ただし T = 0.5, M = 5000, L = 1.0, N = 50, D = 1.0 程度,
初期条件は u(x,0) = 2*x*(1 - x), 境界条件はディリクレ条件とする


実際の計算部分を正しくプログラムにすることができれば,次のようなアニメーションを見ることができます.両端から熱(物質)が奪われ,全体の温度(密度)が下がることがわかります.





今日(7/6)の内容


1次元拡散方程式の数値解法(その2)

次の問題を考えます.先の問題0.1.1と少し違いますので注意してください(初期値と x の範囲が異なる).

問題 0.2


課題2

問題0.2の解を数値的に求めるプログラムを作成してください.π は pi = 3.1415926 等と近似値を与えてください.


問題0.2の数値計算結果例


このように拡散方程式では、初期にできている空間の不均一性のうち、細かい構造(細かい波)から速く減衰して均一化していく事が分かります。


課題3

境界条件をノイマン条件に変更すると、どのような振る舞いをするか調べよ。
ノイマン条件にするには、式(4)の代わりに



とすればよい。


課題4

境界条件を先に述べたディレクレ条件をやや改良した形、
  すなわち 

  すなわち 

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

これを用いて、例えば A = 2.0 , B = 0.2 等とし、拡散方程式(課題2)を数値的に解き、青の様子を動画で観察せよ。


これは例えば、両端の温度を固定したときの熱伝導、両端に濃度の異なる物質の容器をつなげたパイプ中で、粒子の拡散する状況等をシミュレートしてる事になります。





課題5

A = 1 + sin(wt), B = 1 とした場合の様子を、様々な w , D について数値的に解き、その様子を考察せよ。

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の結果との比較をせよ.)。

課題7

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



2変数反応拡散系

 拡散方程式は熱や物質が拡散していく様子を表現しており、それが示す通り、物質や熱は一様にならされる方向に時間変化していきます。しかし例えば2種類以上の物質が反応しながら拡散する場合、そのような一様化が阻害され、むしろ不均一性が増加する場合がありえます。例えば次のブラッセレーターモデル がその一例となっています.

このモデルでは、時刻 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 を変化させた場合も考察してください




差分化の誤差について

先に出てきた,

を繰り返し計算し,近似解を求める方法を陽解法と呼びます.(前の時刻 tn での値から次の時刻 tn+1 の値を直ちに計算できる形の差分方程式を陽公式(explicit scheme)と呼びます.陽公式を使った解法なので陽解法と呼びます.)

天下り的になりますが,拡散方程式の数値解法として陽公式を用いる場合, τと h と D は次の関係を満たす必要があります(安定性条件)

α = Dτ/h2 < 1/2  (S)

ここでは,この安定性条件を破った場合に計算がどのようになるかを見ます.


課題9

課題2のプログラムにおいてのパラメータ n (空間刻み数) やパラメータ M (時間刻み数) 等を調節し,安定性条件(S)をわずかに破って見よ. 

課題4に従い,安定性条件(S)をわずかに破った場合,次のような結果が得られます.(S)の左辺の値が 0.5066 と,1/2 をわずかに越えているにすぎませんが,最終的には真の解からはほど遠いギザギザの状態になってしまいます.1/2 を大きく越えると,計算がおかしくなるのも早まります.(ここでは n および M をα>1/2 となるように調整した.)


上の,ギザギザの解は,直感的に考えても熱方程式の解としてはおかしいことがわかると思います.(スプーンの熱分布がこの様になるとは思えません.)陽解法においては安定性条件を満たさないパラメータで計算すると,計算が不安定となり,上記のような現象が見られます.(ただし,これは,差分方程式の解としては正しいものです.しかし,近似したい熱方程式の近似解とはなっていところが問題なのです.)

注意: 上の計算結果では対称性の崩れ(初期値は π/2 で対称)がみられます.これは,πを近似値として与えたためであると思われます.


この様に,陽解法においては,計算が安定に進むためには時間刻み幅および空間刻み幅が安定性条件(S)を満たさねばなりません.ところが,この条件は計算量の点で問題があります.例えば,時刻0から1まで計算を行うとします.時間方向の格子点の総数は 1/τとなります.計算の手間(計算時間)を減らす為には τ をなるべく大きく取りたい(τ が大きければ少ない繰り返し計算で目的の時刻まで計算できるので).そこで,安定性条件(S)より,

τ ≒ h2/(2D)

とすれば良いでしょう((S)を余裕をもって満たす τ を用いる).一方,D は一定とした場合,精度の良い計算をするために空間に関する分解能を高めたいとします.そこで h を 1/10 倍にしたとします.すると,τ  は上の関係から 1/100 倍せざるを得ません.ということは,時間方向の格子点数が 100倍 に増えることになります.これでは, h を小さくしていったときに,あまりに格子点数が多くなりすぎて計算の効率が悪くなります.例えば今の例では,h を 1/10 倍したため,計算量は 10 倍となり(これは仕方が無い),更に安定性条件(S)を満たすために τ を 1/100 倍したので計算量は更に100倍となりました.つまり,もとの計算量に比べて 1000倍の計算量が必要(1000倍時間がかかる!)となります.

つまり,陽解法では,空間解像度を細かくしたい場合に,時間刻み幅もそれにあわせて変更する必要があります.できることならそういう事を気にせずに(全く気にしなくても良いわけではありませんが・・)自由に空間解像度を上げ下げしたいものです.その要求を満たす解法として陰解法と呼ばれる方法があり,陰解法は無条件安定であることがわかっています(安定な計算に対するαに関する条件が無い).(演習では時間の関係から陰解法についてはできませんが,各人陰解法のプログラムの作成をしてみることをおすすめします.陰解法では連立一次方程式を解く必要がありますが,それらのプログラムは良いプログラミングの練習にもなります.)


戻る