第6回

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

次の１次元拡散方程式（１次元熱伝導方程式とも呼ばれる）の初期値境界値問題を考える。

プログラミングの題材としては、アルゴリズム中の2重ループの使い方が学びどころだろう。

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

問題 0.1.1 について，実際の物理現象としてどのようなものを考えているのかを簡単に説明しよう．問題 0.1.1 の方程式は，拡散方程式とも熱伝導方程式とも呼ばる．拡散とは物が散らばっていく様子，例えばコップ中にインクをたらすとそれは時間とともに散らばり最終的には一様な状態になる．熱伝導とは熱が伝わる様子，例えばスプーンの端を炎で熱すると，熱は徐々に伝わり持っている部分も暑くなり火傷することになる．そういった異なる現象が同じ方程式で表現されるという面白さがある（方程式を導出する際の物理法則が違うが，それらが式として表現した場合に同じである事による）．ここでは，問題 0.1.1 を熱伝導方程式，つまり，熱の伝わる様子を表現する方程式としてみることにする．（ここでの境界条件などが，熱伝導方程式としてみたほうが自然である為．）

ここでは，大変細い棒のようなものの温度分布が時間とともにどのように変化するかといった問題をイメージする．いま，棒の長さは 1 であって，時刻 t における棒の温度分布が関数 u(x,t) で与えられるとする．問題 0.1.1 の第2行目は境界条件と呼ばれ，棒の端の状況を示している．ここでは，すべての時刻において，端の温度は 0 であるとしてある．状況としては，棒の両端を氷で冷やしている，つまり，両端からどんどん熱が奪われていっている状況である．ここで注意が必要なのは，実際の棒では，棒の側面からも熱が奪われるが，この問題 0.1.1 では，それは考えられていない．つまり，棒の側面からは熱が逃げないような状況で，端からのみ熱が逃げる状況を考えている．もちろん，側面からの熱の収支（棒の側面を炎などで暖める状況など）も考えることができ，その場合の方程式は問題 0.1 のようになる（関数　f(x,t)　が熱の収支部分を表現する）．問題 0.1.1 の第4式は，初期条件、初期関数と呼ばれ，観察をはじめる場合のはじめの時間の棒の温度分布を与える．はじめの状態が決まらないと，その後の状態は決まらないであろうことは直感的にわかるだろう．はじめの状態をひとつ決めた場合にその後の状態がただひとつ決まるかどうかは数学として重要な問題である．熱方程式について，そのような性質があるかどうか調べてみるとよだろう（解の一意性）．

さて，この問題を常微分方程式でも行った様に計算機で扱えるよう離散化し，数値的に近似解を求める．ここでは，問題 0.1 1 を例にすすめるが，問題 0.1 に対しても同様である．

さて，数値計算する為にまず(問題 0.1.1)を離散化する必要がある。これはこれまでの常微分方程式の場合と共通で，要は微分を差分で置き換えるのである．分割方法は以下の２種類が考えられる。

ここでは時間方向の離散化は (1) を，空間方向の離散化は (2) を採用する．（これまでの内容で常微分方程式の数値解法を扱ったが，そこでは時間方向の離散化として (1) を採用してきた。空間方向の離散化に (2) を採用する理由は幾つかあるが、境界条件の扱いが比較的自然で、境界条件の切り替えが楽であるという事が大きい。）

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

と書くことにすると、

境界条件は（問題 0.1.1）の第二式だが，まずは境界条件について触れる必要がある．ここでは，境界の処理の為に架空のセルを用意することにする．

様々な境界条件があるが，今回の問題 0.1.1 における境界条件は Dirichlet 条件（この場合各時間 0 という値に固定されている．一般の場合，例えば a という値に固定されている場合には，次の式の右辺 0 が a となる．）と呼ばれるものである．空間離散化の方法として (1) を採用した事を思い出すと，

となればよい．すなわち，

となる。

準備はできたので、問題 0.1.1 の離散化を行う（微分を差分で置き換える）。

は次の差分方程式で近似される（２階微分の差分化については，このページ最後のスライドを参照）．

これを書き換えると、

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

次に拡散方程式を計算機にて解くためのアルゴリズムを考えよう．解を求める時刻の上限を T とし，適当な自然数 M を定めて，τ=T/M とする．問題 0.1.1 を解くアルゴリズムは次のようになる．ここでは、適当な時間感覚でグラフを重ね描きすることとする。

アルゴリズムにおいて、 k の繰り返し文が時間に関する繰り返し、j に関する繰り返し文（赤色）が空間に関する繰り返しである。各時間を止める毎に空間の繰り返しという2重の繰り返し文になっている。

アルゴリズム・1次元拡散（熱伝導）方程式（陽解法）

（１）初期パラメータ設定

N, M, T, D を設定する

h := 1/N, τ:=T/M, α:= Dτ/h2

x[N+2], u[N+2], new_u[N+2] を配列（リスト）として用意

（２） 初期値設定

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

xj = (j - 0.5)*h

uj := 2xj(1-xj)

を繰り返す

t=0におけるグラフを描画 

（３）

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

u0 := -u1, uN+1 = -uN 　　（境界条件）

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

new_uj := αuj-1 + (1 - 2α)uj + αuj+1

を繰り返す

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

uj := new_uj

を繰り返す

一定間隔毎にグラフ描画

を繰り返す

課題1（演習）

上のアルゴリズムを参考にプログラムを完成させよ．（ただし T=0.2, M=1000, N=50, D=1.0 とする。）

上の動画では、最後の方のグラフの仕上げでもたもたしてしまったが、こういうことはよくあるのでそのままにしておく。最終的に ug というリストを用意したが、今考えるとグラフを描く直前のみ u[0] と u[N+1] に 0 を入れるだけでよかった。境界条件は毎回設定するので、これで問題ない。ただし、掲載コードはそのままにしておく。次のおまけ演習動画で修正する。

境界条件について

また課題1において境界条件はディリクレ(Dirichlet)境界条件を扱っているが、ノイマン(Neumann)境界条件の例としては、下の方のおまけ演習2がある。上のコードで言えば、

u[0] = -u[1]

u[N+1] = -u[N]

の部分を、

u[0] = u[1]

u[N+1] = u[N]

と変更するだけである(上の方の境界の処理のスライド参照)。やってみよう(グラフで両端を表示する部分で問題が生じるが、自分で修正して欲しい)。この場合、全体の熱量（物質量）が保存される。ノイマン境界条件は、断熱壁とか no flux 条件とも呼ばれる。

また、次のように変更すると、周期的境界条件となる。

u[0] = u[N]

u[N+1] = u[1]

但し、上記のコードをそのまま変更すると、初期条件が境界条件を満たさない（両端で微分が0ではない）。しかしながら、計算自体はうまく行く。境界の扱いは問題によっては結構繊細なので注意が必要だが、何を見るための数値計算であるかを意識することが大事であろうと思う。初期条件が境界条件を満たす条件を適合条件と呼ぶこともある。

おまけ演習1

課題１ではグラフの重ね描きをしたが、matplotlib の機能でアニメーション表示を行う。

余談：Animated GIF 形式での保存には、以前は ImageMagick というアプリが必要であったが、今は pillow というものが Anaconda に標準で入っており、別にアプリを入れなくてもできるようになったようだ。mp4 形式で保存する場合には、マシンに ffmpeg というアプリが入っている必要がある。 macOS の下準備 にある MacPorts で導入できる。他のOSは、それぞれの方法があるだろう。

上の動画の中で、時刻 t を画面内に入れることに失敗したが、次のコード Code6-1-b.py で実現できた。赤色部分が変更点である。もっとスマートな方法があるのかもしれないが、matplotlib の使い方の問題であって本論では無いのでこれで良しとする。

課題２

問題 0.2 の解を数値的に求めるプログラムを作成せよ。（ただし T=0.5, M=1000, N=80 とする。）

安定性条件

先に出てきた，

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

天下り的になるが，拡散方程式の数値解法として陽公式を用いる場合， τ と h と D は次の関係を満たす必要があります（安定性条件）．（何故そうなるかは参考書を見て欲しい。）

α = Dτ/h2 < 1/2　　（S）

ここでは，この安定性条件を破った場合に計算がどのようになるかを見てみる．つまり、次の課題では、やってはいけないといわれていることをわざとやってみる。

課題3（演習）

課題2 で作成したプログラムについて、安定性条件 （S）を満たしていることを確認し、続いて (S) 破ったパラメータで実行するとどのような結果になるか確かめよ。

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

注意：　上の計算結果では対称製の崩れ（初期値は π/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倍時間がかかる！）となる．

つまり，陽解法では，空間解像度を細かくしたい場合に，時間刻み幅もそれにあわせて変更する必要がある．できることならそういう事を気にせず（全く気にしなくても良いわけではないが・・）自由に空間解像度を上げ下げしたい．その要求を満たす解法として陰解法と呼ばれる方法があり，陰解法は無条件安定であることがわかっている（安定な計算に対するαに関する条件が無い）．ただし、ここでは陰解法については扱わない予定である。

但し、陽解法は単純であることと、近年の並列コンピュータの発展とともに見直されているという事実がある。問題によっては、陽解法で十分であることも多く、特に空間3次元の問題等を並列コンピュータで解く場合には、単純な陽解法を用いて並列度を上げて高速化をはかる場合がある。

おまけ演習2

1次元拡散（熱伝導）方程式の数値計算ができるようになったので、ここでシステム（連立系）に拡張してみる。ここでは、 Gray-Scottモデルと呼ばれる反応拡散系を扱おう。次の文献（大昔に私が書いた）を参考にしつつコードを書く。

• AUTOによる大域ダイナミクスの解析の試み, 上山大信, 数学第61巻第3号(2009), pp.316-326.

空間2次元においては、次の様な自己複製スポットとよばれるパターンが観察される方程式だ。これの、空間1次元版である、自己複製パルスが観察されるパラメーターでのシミュレーションを行う。パラメーター等は、上の文献にあるものを使う。

拡散（熱伝導）方程式の解は、本ページの前半の話からもわかるように、より平坦な形になるのであるが、拡散方程式の一種である Gray-Scottモデルでは模様ができあがる。不思議に思ってもらえると幸いである。

Gray-Scottモデルの解の動きを見ると、やはりゾクゾクする。全く単純な方程式であるが、このような複雑でかつ意味ありげな解が見られるのは面白い。そこにはまったのであるが。ちなみに、文献中の図や上の2次元シミュレーションなどは全てC言語で書いている。2000年頃の事なので、コンピュータの速度が遅く、当時Pythonで数値計算ということは考えもしなかった。

１次元熱方程式の数値解法（その３）

　次の問題を考えます．

課題4

u(x,t) = sin(3πx)exp(-(3π)2t)　　　　(5)

が，問題 0.3の解であることを確かめよ．

ここでの目標

問題 0.3を数値計算で解いた差分解と微分解 (5) を比べることで，差分法による計算が，どの程度解析的に求まった解 (5) と近いのかを調べる．直感的には h や τ を小さくすると，近似の程度はよくなると思われるが，実際どの程度よくなるかを確かめる．

問題 0.3 を数値計算で解いた解を  u kj とし，解 (5) の格子点状での値を，u(xj, tk) と書くとき，それらの相対誤差を次のように定義する．

もしも，数値計算で求まった解と，解析的に求まった解が完全に一致していれば err(N) は 0 となるが，実際にはならない．N （空間刻み数）を変化させて，相対誤差がどのように変化するかを見る．

err(N) は k （時刻）を止めるごとに求まるが，ここでは k=M とすることにする．ただし M とは T を時間の上限とした場合の　M=T/τ．

N を色々かえてみて err(N) を求めればよいが，上で紹介した安定性条件（Ｓ）があるため，N を独立に自由に動かすことはできない．

α = Dτ/h2 < 1/2　　（Ｓ）

h は N と L (N:空間分割数 ,L:領域長)から決まるが，いま L=1 とするので，h は N からのみ決まる．（Ｓ）を満たしていなければならないので，ひとまず，

α=1/3

とする(ここに 1/3 は 1/2 より小さければ何でも良い訳で特別意味は無い)．すると，

τ=h2/3

となる．t の上限である T を与えれば次の式で時間刻み数 M を得ることができる．

M=T/τ

先に示したプログラムでは，N, M, T から τ, αが決まっていたが，今回は α, n, T からM, τが決まるようにする．

解を求める時刻の上限を T ，空間刻み数を N とし，αは1/3とすることにすると，アルゴリズムは次のようになる．前回のプログラムを変更すればよいが，新たに U という配列が必要になる（微分解の値をいれる配列．差分解と比較するために必要）．

アルゴリズム

（１）初期パラメータ設定

N, T, α を設定する

α:=1/3, h := 1/N, τ:=αh2, M:=T/τ

（２） 初期値設定

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

uj := sin(3πxj)

を繰り返す

（３）

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

tk := kτ

u0 := -u1, uN+1 = -uN 　　（境界条件）

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

new_uj := αuj-1 + (1 - 2α)uj + αuj-1

を繰り返す

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

uj := new_uj

を繰り返す

を繰り返す

（４）解析解と数値解の比較

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

xj := (j - 0.5)*h

Uj := sin(3πxj)*exp(-(3πxj)T)

を繰り返す

時刻 T におけるerr(N) をもとめ， h と err を画面に表示する

課題5

N = 50, 100, 200, 400 についてそれぞれ err(N) を求め，h と err(N)　の間にどのような関係があるか調べて見よ．

この問題のポイントは max を求める部分。恐らく関数として用意されているが、繰り返し文を使って求めてみると良い。例えば、配列として、でたらめな数値が入っている a[0]〜a[10] があったときに、その配列の中の数値の絶対値が最大の値を求めるには、次の様なコードが考えられる。

max_a = abs(a[0])

for i in range(1, 11):

if abs(a[i]) > max_a :

max_a = abs(a[i])

print(max_a)

基本中の基本だが、結構これができない人も多い。確認すべきポイントだろう。

課題5 から得られる結果を書いておく。

結果として，h を 1/2 にすると err(N) は 1/4 となることがわかる（つまり，err(N) が h*h に比例しているという事実）．両対数グラフを書けば良いが、例えば次のようになる．本ページの最後の部分に、2階微分係数の離散化についてのスライドがあるが、そこで無視した部分が O(h*h) であるから、真の解との誤差は h*h に比例すると考えられ、実際そのようになっていることを確かめることができる．もしも，x*x の直線と違う傾きの直線なり曲線になった場合は，プログラムのミスである可能性が高い．この様に，プログラムのミスを発見する一つの有用な手段になりうる．（これは， matplotlib のグラフではないが， matplotlib でも同様のグラフが作成できる。赤線の 0530.data なるものが、課題5 の結果に対応する直線で、緑色の直線が x*x に対応する。）

差分近似について