今回のターゲット：１次元熱拡散方程式
前回は数値計算に関して、基本に立ち返った勉強をしました。
今回は物理現象に戻って、１次元の熱拡散方程式の数理計算を試みます。
まずは熱拡散方程式と差分法について勉強します。
勉強の際に参考にさせていただいた書籍もご紹介させていただきます。具体的なPythonのコードも載っており、大変勉強になる１冊です。
１次元熱拡散方程式
１次元熱拡散方程式は、以下の式で表すことができます。
$$
\cfrac{\partial{u\left(x, t\right)}}{\partial{t}} = \alpha\cfrac{\partial{^2}{u\left(x, t\right)}}{\partial{x^2}} \tag{1}
$$
ここで、$\alpha$ [m$^2$/s]は熱拡散率であり、熱伝導率 $k$ [W/m•K]、物体の密度 $\rho$ [kg/m$^3$]、比熱 $c_{p}$ [J/kg•K]を使って、
$$
\alpha = \cfrac{k}{\rho{c_{p}}} \tag{2}
$$
と表すことができます。
式(1)は２つ以上の変数を持つ微分方程式、すなわち偏微分方程式になっています。もっと言えば、時間 $t$ と位置 $x$ の偏微分方程式となります。
このような偏微分方程式を数値計算するために、今回は中心差分を使います。
差分法
差分法
まずは差分法について勉強しました。ここではわかりやすいように以下の式を比べてみます。ただし、$x_{i+1} = x_{i} + h$です。
$$
\cfrac{\partial{f}}{\partial{x}} = \lim_{h\to 0}\cfrac{f\left(x + h\right) – f\left(x\right)}{h} \tag{3}
$$
$$
\left(\cfrac{\partial{f}}{\partial{x}}\right)_{i} = \cfrac{f\left(x_{i+1}\right) – f\left(x_{i}\right)}{h} \tag{4}
$$
上の式は微分の定義式です。下の式が差分法（前進差分法）の式になります。
イメージとしては、”連続である上の式の $x$ を離散化したものが下の式になる” という感じですか。
また、前進差分法は $f(x_{i} + h)$ を$h$についてTaylor展開したものから得られます。つまり、
$$
f\left(x _{i} + h\right) = f\left(x_{i}\right) + h\left(\cfrac{\partial{f}}{\partial{x}}\right)_{i} + O\left(h^{3}\right) \tag{5}
$$
これを $\left(\cfrac{\partial{f}}{\partial{x}}\right)_{i}$ について整理すれば、式(4)が得られます。
中心差分
中心差分は $f\left(x_{i} \pm h\right)$ を$h$について $\pm h$ 両方の場合についてTaylor展開したものから得られます。つまり、
$$
f\left(x _{i} \pm h\right) = f(x_{i}) \pm h\left(\cfrac{\partial{f}}{\partial{x}}\right)_{i} + \cfrac{h^2}{2}\left(\cfrac{\partial{^{2}f}}{\partial{x^2}}\right)_{i} + O\left(h^{3}\right) \tag{6}
$$
について、$+$ の場合と $-$ の場合の差を取ると、以下の式が得られます。
$$
\left(\cfrac{\partial{f}}{\partial{x}}\right)_{i} = \cfrac{f\left(x_{i+1}\right) – f\left(x_{i – 1}\right)}{2h} \tag{7}
$$
これが中心差分の式となります。式を見ればわかるように、精度は2次となります。
さらに式(6)の$+$ の場合と $-$ の場合の和を取って整理すると、
$$
\left(\cfrac{\partial{^{2}f}}{\partial{x^2}}\right)_{i} = \cfrac{f\left(x_{i+1}\right) -2f\left(x_{i+1}\right) + f\left(x_{i – 1}\right)}{h^2} + O\left(h^{2}\right) \tag{7}
$$
となります。これが熱拡散方程式を解く際に使う2階微分に関する中心差分の式となります。
最後に
いかがでしたでしょうか。今回は1次元熱拡散方程式を数値計算するために、差分法について勉強してみました。
もうお腹いっぱいなので、実際にPythonで計算するのは次回にしたいと思います。
それでは今回はこのあたりで終わりにします。
他にも色々な記事があるのでぜひご覧ください。
コメント