第11回(2011年6月29日)


今回から3回にわたって偏微分方程式の数値解法を紹介します。まず拡散方程式(熱伝導方程式)をとり扱い、その後、反応拡散方程式、波動方程式の取り扱いへと移っていく予定です.拡散方程式(熱伝導方程式)の導出等に関しては,講義にて詳しい説明があったと思います.

お知らせ: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 (空間)に関する微分も含まれているため、こちらに対しても適当な差分化を行う必要があります。時間(t)、空間(x) の離散化(分割)は、以下のように進めていきます。



テーラー展開の結果を利用して, U(x) の1階、2階微分を以下のように差分近似します。







常微分方程式の時間 (t) 同様、積分区間 (x と t) を離散化する必要があります。離散化の方法は以下のように2種類あります。



後々便利なので、空間、時間共々以下の離散化方法を用います。



上記の事柄を用いると、拡散方程式は以下のように離散化されます。





この式は見てのとおり,時刻 tk での u から次の時刻 tk+1 での u の値が計算できる形をしています.ですので繰り返し文を適切に使ってプログラムを構成すれば、与えられた初期条件から順々に計算する事で、拡散方程式が数値的に解かれる事に成ります。。。と言いたい所ですが、実はもう少し条件が必要です。


境界条件

空間に関しても端を適度に扱う必要があり, そのために境界条件が与えられます. 境界条件とは例えば問題0.1, 問題0.0.1の第二式で与えられるような条件で、空間(x)の端の満たす条件です. ここでは,境界の処理の為に, という架空のセルを用意することにします. 良く利用される境界条件として、以下の3つの条件があります。





それぞれ3つの境界条件は稲野用なイメージです。
 . ディレクレ条件: 両端の状態(温度や物質の濃度)が常に一定. (固定端)
 . ノイマン条件: 両端に鏡がついていて、鏡の向こう側とやり取りしているかのような条件. (自由端)
 . 周期境界条件: システムがループ状に接続


以上の境界条件と更に初期条件が整う事で、ようやく以下のように偏微分方程式を解く事ができるようになります.







ではおさらいも兼ねて、問題0.1.1を離散化して解いてみましょう.(dt → τ , dx → h と成っています.)

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


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



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

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

問題 0.2


課題2

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


問題0.2の数値計算結果例


課題3

境界条件をノイマン条件に変更すると、どのような振る舞いをするか調べよ。
ノイマン条件にするには、式(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倍時間がかかる!)となります.

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


戻る