第7-4回　二次元拡散方程式（実装までの式変形）[python, julia, fortran]

• English (United States)
• 日本語

$$\rho C \frac{\partial T}{\partial t}=\frac{\partial}{\partial x} \left( \lambda \frac{\partial T}{\partial x} \right) + \frac{\partial}{\partial y} \left( \lambda \frac{\partial T}{\partial y} \right)$$

時間一次精度、陰解法で離散化します。第7-2回で行った方法を用います。

左辺は今までと同じです。

$$\rho C \frac{\partial T}{\partial t} \simeq \rho C \frac{\Delta T^n_{i,j}}{\Delta t}$$

$$\Delta T^n_{i,j}= T^{n+1}_{i,j}- T^{n}_{i,j}$$

右辺第一項は次のように書きます。詳しくは第7-2回をご覧ください。

$$\frac{\partial}{\partial x} \left( \lambda \frac{\partial T}{\partial x} \right) \simeq \frac{\partial}{\partial x} \left( \lambda \frac{\partial T^{n+1}}{\partial x} \right)$$

$$ \simeq \frac{ \lambda_{i+1/2,j} }{\Delta x^2} (\Delta T_{i+1,j}^n +T_{i+1,j}^n-\Delta T_{i,j}^n -T_{i,j}^n )$$

$$-\frac{ \lambda_{i-1/2,j} }{\Delta x^2} (\Delta T_{i,j}^n +T_{i,j}^n -\Delta T_{i-1,j}^n -T_{i-1,j}^n) $$

右辺第二項も同様に行います。

$$\frac{\partial}{\partial y} \left( \lambda \frac{\partial T}{\partial y} \right) $$

$$\simeq \frac{ \lambda_{i,j+1/2} }{\Delta y^2} (\Delta T_{i,j+1}^n +T_{i,j+1}^n-\Delta T_{i,j}^n -T_{i,j}^n )$$

$$- \frac{ \lambda_{i,j-1/2} }{\Delta y^2} (\Delta T_{i,j}^n +T_{i,j}^n -\Delta T_{i,j-1}^n -T_{i,j-1}^n) $$

まとめると次のようになります。

$$\rho C \frac{\Delta T^n_{i,j}}{\Delta t}=\frac{1}{\Delta x^2}(\lambda_{i+1/2,j} (\Delta T_{i+1,j}^n -\Delta T_{i,j}^n)-\lambda_{i-1/2,j}(\Delta T_{i,j}^n -\Delta T_{i-1,j}^n))$$

$$+ \frac{1}{\Delta y^2}(\lambda_{i,j+1/2} (\Delta T_{i,j+1}^n -\Delta T_{i,j}^n )-\lambda_{i,j-1/2}(\Delta T_{i,j}^n -\Delta T_{i,j-1}^n))+R_{i,j}^n$$

$$ R_{i,j}^n =\frac{1}{\Delta x^2}(\lambda_{i+1/2,j} (T_{i+1,j}^n – T_{i,j}^n)-\lambda_{i-1/2,j}( T_{i,j}^n – T_{i-1,j}^n))$$

$$+ \frac{1}{\Delta y^2}(\lambda_{i,j+1/2} (T_{i,j+1}^n – T_{i,j}^n )-\lambda_{i,j-1/2}( T_{i,j}^n – T_{i,j-1}^n))$$

さらに、反復数kを用いてガウスザイデル法の形に変形します。詳しくは第7-3回を参照してください。

$$ \left ( \frac{\rho C }{\Delta t}+ \frac{ \lambda_{i+1/2,j} }{\Delta x^2} + \frac{ \lambda_{i-1/2,j} }{\Delta x^2}+ \frac{ \lambda_{i,j+1/2} }{\Delta y^2} + \frac{ \lambda_{i,j-1/2} }{\Delta x^2} \right ) \Delta T^{k+1}_{i,j} =\frac{1}{\Delta x^2} \left[\lambda_{i+1/2,j} \Delta T_{i+1,j}^k -\lambda_{i-1/2,j} \Delta T_{i-1,j}^k \right]$$ $$+ \frac{1}{\Delta y^2}\left[\lambda_{i,j+1/2} \Delta T_{i,j+1}^k -\lambda_{i,j-1/2} \Delta T_{i,j-1}^k \right]+R_{i,j}^n$$

$$ \Delta T^{k+1}_{i,j} =\left ( \frac{1}{\Delta x^2} \left[\lambda_{i+1/2,j} \Delta T_{i+1,j}^k -\lambda_{i-1/2,j} \Delta T_{i-1,j}^k \right]+ \frac{1}{\Delta y^2} \left [ \lambda_{i,j+1/2} \Delta T_{i,j+1}^k -\lambda_{i,j-1/2} \Delta T_{i,j-1}^k \right]+R_{i,j}^n \right ) $$ $$ / \left ( \frac{\rho C }{\Delta t}+ \frac{ \lambda_{i+1/2,j} }{\Delta x^2} + \frac{ \lambda_{i-1/2,j} }{\Delta x^2}+ \frac{ \lambda_{i,j+1/2} }{\Delta y^2} + \frac{ \lambda_{i,j-1/2} }{\Delta x^2} \right ) $$

これで式としては実装できる形となりました。（ガウスザイデル法の場合，格子の並び順によっては一部kではなくk+1と書きますが，理解しやすい形のためこのまま行きます。）

しかし、反復法の収束性を見極める基準をまだ設けていません。

今回は2-ノルムが\(10^{-4}\)以下となったら計算を終了させたいと思います。

定常反復法は一定割合で収束していきます。そこで、真の解との差がほとんどなくなったら収束したとみなします。その基準を2-ノルムが\(10^{-4}\)以下としました。

下記のように行列があったとします。\(\bf A\)は\(i × i\)行列, \(\bf x\)と\(\bf b\)は\(i\)の列行列です。

$$\bf A \bf x =\bf b$$

反復数\(k\)のときを考えると誤差があるため、以下の式の右辺は0となりません。

$$\bf A \bf x^{k} -\bf b=\bf e$$

そこで、2-ノルムをとります。

$$2-norm=\sqrt{\sum_{k=0}^i (a_{k,j}x_{k}-b_{k})^2}$$

いつでも使えるように\(\bf\)を用いて規格化します。

$$2-norm=\sqrt{\frac{\sum_{k=0}^i (a_{k,j}x_{k}-b_{k})^2}{\sum_{k=0}^i (b_{k})^2}}$$

以上の値を用いて

$$2-norm \le 10^{-4}$$

となったとき計算終了とします。

実装までの式変形がすべて終わったので、次回は実際にコードに落とし込みたいと思います。

*コードは第7-5回にあります。

By hide