第12回(2009年7月15日)


今回と次回の2回にわたって拡散方程式(熱伝導方程式)の数値解法を行います.拡散方程式(熱伝導方程式)の導出等に関しては,講義にて詳しい説明があったと思います.


演習内容 


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

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

問題 0.1

 

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

問題 0.1.1


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

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


さて,波動方程式に対して行ったのと同様に計算機で扱えるよう離散化し,数値的に近似解を求めます.ここでは,問題0.1.1を例にすすめますが,問題0.1に対しても同様に離散化できます.

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

と書く事にしましょう.

波動方程式と同様にして,

  (1)

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

これを書き換えれば,

      (2)

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

     (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」等とする.)

(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-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-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), 境界条件はディリクレ条件とする


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



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

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

問題 0.2


課題2

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


問題0.2の数値計算結果例


課題3

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



とすればよい。


レポート課題

課題1〜3のどれか, もしくは下の「おまけ」についてのプログラムを rep0715.c というファイル名で report フォルダにおいてください. (プログラムの先頭にコメント文で名前と学籍番号と課題の番号を書く事. またプログラムの説明もコメント文として適宜記入する事.)
期末レポート作成には、拡散方程式の理解が必須になってきます. しっかり習得してください.

自動回収日時:2009年7月21日(火)午後4時



先に出てきた,

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

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

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

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


課題4

課題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倍時間がかかる!)となります.

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


おまけ(フラクタル図形を描く)

(拡散方程式が習得できて余裕のある方、トライしてみてください.)

(1) Ukj が以下の差分方程式に従うとする.

   Ujk+1 = F( Ukj-1 ,  Ukj ,  Ukj+1 )
ただし k = 0, 1, 2, ...M (正整数), j = 1, 2, ..., N (正整数) 


ここで F( Ukj-1 ,  Ukj ,  Ukj+1 ) = Ukj + α( Ukj-1 -2Ukj + Ukj+1 )
と与えられれば、先に扱った拡散方程式の差分化したものになります.


しかし一般的に F( Ukj-1 ,  Ukj ,  Ukj+1 ) は、もっといろいろな形のものが考えられて良いはずです。
そこで今回、以下のような差分方程式系を考えてみます。

 Ukj は 0 もしくは 1 の2値のみをとり、次の差分式に従う.

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

   ただし
       F(0, 0, 0) = 0, F(0, 0, 1) = 1, F(0, 1, 0) = 0, F(0, 1, 1) = 1,
       F(1, 0, 0) = 1, F(1, 0, 1) = 0, F(1, 1, 0) = 1, F(1, 1, 1) = 0
 もしくは
     F( Ukj-1 ,  Ukj ,  Ukj+1 ) = 1  ( Ukj-1 + Ukj+1 = 1 の時 )
     F( Ukj-1 ,  Ukj ,  Ukj+1 ) = 0  ( Ukj-1 + Ukj+1 = 1 以外の時 )
(どちらも同じです.)

境界条件: Uk0 = 0, UkN+1 = 0

初期条件: U0j = 1 ( j = N/2 ), U0j = 0 ( j = N/2 以外 )

とする. その解の様子を図示せよ.

表示方法は, 縦軸 k, 横軸 j, とし, Ukj = 1 のとき, g_circle, g_marker 等で ( j, k ) に丸(点)を描く等とすること.
( M = 200, N =400 程度でこんな絵が描けます.(点は半径 0.5 の円.)

 

(2)「フラクタル」について調べよ

(3) Cellular Automata (セルオートマトン) について調べよ




期末レポート予告

期末レポートを提出してもらいます。提出〆切は7月末の予定です。(詳細は次回)

課題は以下の通り。

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


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

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

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

に従うとする.
ここで初期条件を

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


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

このブラッセレーター系含めた偏微分方程式系の取り扱いは 演習13回〜15回(特に14回〜15回)で行います。まず前回の波動方程式、今回(次回)の拡散方程式を、しっかり習得してください.


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

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

u0 = -u1, uN+1 = -uN

等としていたところを,

u0 = u1, uN+1 = uN

 等と変更する


提出方法:

(1) は数学事務室に提出、(2) はいつも通りプログラムを自動回収します.

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




戻る