1次元熱伝導方程式を差分法で解く（陽解法）
はじめに
今回は1次元熱伝導方程式を差分法で解きます．時間発展は陽解法（explicit method）で行います．陽解法とは，現在の時間ステップの値だけで次の時間ステップの値を決める時間発展の仕方のことです．
支配方程式
まず，1次元の熱伝導方程式と境界条件，初期条件を示します．未知数を温度と考えると熱伝導方程式になりますし，未知数を濃度と考えると拡散
方程式と考えることもできます．また，微分方程式を考える際は，境界条件と初期条件が重要なので常に注意して下さい．境界条件と初期条件のいずれかが不適切であるせいで発散してしまう例もあります．方程式と境界条件，初期条件は
$$
\frac{\partial u}{\partial t} - \alpha \frac{\partial^2 u}{\partial x^2} = 0，0 \leq t \leq T，0 \leq x \leq L，\\
u(t,0) = 0，u(t,L) = 0，u(0,x) = u_0(x)
$$
です．ここで，$${u = u(t,x)}$$は未知数（温度），$${t}$$は時間，$${\alpha}$$は熱拡散率，$${x}$$は空間座標，$${T}$$は最終時刻，$${L}$$は考えている領域の大きさを表す定数，$${u_0(x)}$$は初期条件をあらわす関数です．今回はDirichlet条件（未知数の値を境界で直接指定すること）を考えます．次回以降は，Neumann条件（未知数の勾配を境界で直接指定すること）やRobin条件（未知数の値と勾配の和を境界で直接指定すること）も考えてみようと思います．
離散化
空間を等分割して節点を配置します．そのとき，境界の節点が計算領域と一致するようにします．また，時間方向にも等分割に分割し，各々の離散的な時間をステップと呼びます．これらの離散点を用いて離散化していきます．$${i}$$番目の節点における，$${n}$$ステップ目における未知数の値を$${u_i^n}$$とあらわします．まず，時間項は陽解法なので
$$
\left( \frac{\partial u}{\partial x} \right)_i^n = \frac{u_{i}^{n+1} - u_i^n}{\Delta t}
$$
と離散化します．一方，拡散項は2階微分なので
$$
\left( \frac{\partial ^2 u}{\partial x^2} \right)_i^n = \frac{u_{i+1}^n -2 u_i^n + u_{i-1}^n}{(\Delta x)^2}
$$
と離散化します．差分の導出はこちらを参考にして下さい．
ここで，拡散項の離散化に時刻$${n}$$ステップ目の値を用いていることに気を付けて下さい．こうすることによって，陽的に次の時刻の値を計算できるようになります．拡散項の離散化に時刻$${n+1}$$ステップ目の値を用いることも可能で，これは陰解法（implicit method）と呼ばれています．
さて，結局離散化された方程式は
$$
\frac{u_{i}^{n+1} - u_i^n}{\Delta t} = \alpha \frac{u_{i+1}^n -2 u_i^n + u_{i-1}^n}{(\Delta x)^2}
$$
となります．$${n+1}$$ステップの値について整理すると
$$
u_{i}^{n+1} = u_i^n + \frac{\alpha \Delta t}{(\Delta x)^2} (u_{i+1}^n -2 u_i^n + u_{i-1}^n)
$$
を得ます．これが熱伝導方程式を陽解法で解く場合の離散化式です．左辺に未知数である$${n+1}$$ステップ目の値があり，右辺は既知である$${n}$$ステップ目の値から構成されているので，逐次次ステップ目の値を計算していくことができます．
さて，上式が成り立つのは実は領域の内部だけです．今，節点番号を$${1}$$から$${Nx}$$としましょう．すなわち，$${Nx}$$個に分割したとします．このとき，上式に$${i = 1}$$を代入すると，拡散項のところで$${u_0^n}$$という定義されていない量が出てきて困ってしまいます．実は，境界の値はDirichlet条件より
$$
u_1^n = 0
$$
と直接決まります．これは，$${i = Nx}$$も同様です．以上から，結局離散化式を解く必要があるのは$${i = 2 \sim Nx-1}$$であることがわかります．ちなみに，このようにDirichlet条件を用いて境界点の値を直接指定できるのは，境界の節点が計算領域と一致するように設定したからです．次回以降，境界の格子が半分だけ計算領域とずれている場合も扱います．こちらは流体計算の分野でスタガード格子（staggered grid）と呼ばれる非常に重要な格子だからです．
さて，これで$${n}$$ステップ目の値から$${n+1}$$ステップ目の値を計算する方法がわかりました．計算の最初の$${1}$$ステップ目の値は，初期条件$${u_0(x)}$$からわかります．
数値計算例
ここでは，$${L = 1}$$，$${T = 10}$$，初期条件を
$$
u_0 (x) = \sin \left( \frac{\pi x}{L} \right)
$$
とします．以下が計算結果です．
sin波がだんだんとけていく様子が観察できます．
おわりに
今回は1次元熱伝導方程式を差分法，陽解法で解きました．次は，陰解法で解いてみようと思います．
参考文献
・登坂宣好，大西和榮『偏微分方程式の数値シミュレーション 第2版』
差分法だけでなく，有限要素法と境界要素法も解説されている．1次元の場合の楕円型方程式，放物型方程式，双曲型方程式をそれぞれ差分法，有限要素法，境界要素法による離散化を説明し，次に同様のことを2次元で行うという非常にシステマティックな構成です．第5章では粘性流れのシミュレーションも扱われており，定常移流拡散方程式，Burgers方程式，Navier-Stokes方程式の差分法，有限要素法，境界要素法による離散化が説明されています．非常に素晴らしい名著なのですが，現在中古が高騰しています．安価で手に入るのは内容が少ない1版なのでご注意ください．表紙が違います．
・梶島岳夫『乱流の数値シミュレーション 改訂版』
差分法でNavier-Stokes方程式を以下に離散化すべきであるかを非常に丁寧に解説しています．差分の基礎からはじめて乱流の数値計算まで到達するすごい本です．おすすめです．
この記事が気に入ったらサポートをしてみませんか？