1次元熱伝導を陰解法で解く ヤコビの反復法

2019年10月5日非圧縮性流体

次は陰解法で

プログラムは、以下の書籍の内容やプログラムを参考にして作ってます。プログラムは、一部改修・デバッグして使用しています。超いい本なので、ぜひ読んでください。

プログラムは朝倉書店のホームベーじからダウンロードできます。書籍名で検索してみてください。

コメント書いたり修正したプログラムも、githubにアップしておきます。

前回、1次元熱伝導方程式を陽解法解きましたが、今回は陰解法で解きます。同様の問題、

$$
\begin{align}
&\frac{ \partial u }{ \partial t}  = k\frac{\partial^2 u}{\partial x^2} \qquad (k > 0, t > 0,0<x<1) \\\\
& u(0,t) = u(1,t) = 0 \\\\
& u(x,0) = f(x)
\tag{1} \label{eq:scale-factor-1}
\end{align}
$$

を考えます。

前回と同じように、以下のようにメッシュを切ります。

J-1個の格子、J個の格子点に離散化します。格子点の番号を左から1,…,Jとし、n時点での、\(u\)の値を、\(u_{1}^{n},u_{2}^{n},..,\)、と名前をつけていきます。代表点を、\(u_{j}^{n}\)とします。

なので、番号はx=0の点で格子番号は1, x=1の点で格子番号はJ, t=0の点で時間の方の格子番号は1になってます。ややこしいです。

格子幅を等間隔で\(\Delta x, \Delta t \)とします。右上の添字は時間を、右下の添字は空間を表すのが慣例なので、そのようにし、時間項をオイラー陰解法で、2階微分を2次中心差分で離散化します。

$$
\begin{align}
&\frac{u_{j}^{n+1} – u_{j}^{n}}{\Delta t} = k \frac{u_{j-1}^{n+1} – 2u_{j}^{n+1} + u_{j+1}^{n+1}}{(\Delta x)^2} \\\\
& – r u_{j-1}^{n+1} + (1+2r) u_{j}^{n+1} – r u_{j+1}^{n+1} = u_{j}^{n} \qquad (r = k \Delta t/ (\Delta x)^2) \\\\
&(j=2,,…,J-1, n=1,2,…)
\tag{2} \label{eq:scale-factor-2}
\end{align}
$$

境界条件と初期条件も前回と同様に考えます。境界条件は\(u^{n+1}\)について考えます。

$$
\begin{align}
& u(0,t) = u(1,t) = 0 \\\\
& u(x,0) = f(x) \\\\\
& より、\\\\
& u_{1}^{n+1} = u_{J}^{n+1} = 0 \qquad (n=1,2,…)\\\\
& u_{j}^{1} = f(x_j) \qquad (j=1,2,…,J)
\tag{4} \label{eq:scale-factor-4}
\end{align}
$$

となります。以下のような関係性になっています。陽解法と大きく違うのは、単純には次のステップの値を求められない関係性になっていると言うことです。

陽解法で解けるのに、なぜわざわざこんなことをするのかと思うかもしれませんが、陰解法の方が圧倒的に正しい解を求めやすいからです。具体的に言うと、よく数値解析の場合は発散と言う現象が起こります。明らかにありえない解や解が振動してしまう等のことがおきます。それは空間の離散化の\(\Delta x\)が大きすぎたり、時間の離散化の\(\Delta t\)が大きすぎたり、はたまた微分方程式の差分化する時の精度が問題だったり、色々です。

また、解く問題によっても全然変わります。熱伝達率が大きかったり熱の差が大きい問題を解く時は、空間も時間も細かくせねばなりませんし、全体的に温度があまり変わらないなら大きくしても構いません。ただ細かくすればいいと言う問題ではなく、計算時間とトレードオフなので、適切な設定をすることが大事です。

試しに前回の陽解法のプログラムで\(\Delta t\)を大きくしたり、メッシュ数を減らして、\(\Delta x\)を大きくしてみてください。

その発散と言う現象が陰解法の方が圧倒的に起こりにくいです。数学的に証明もできるのですがややこしくなるので簡単な説明に留めます。
陽解法では過去の値から素直に求めているので、誤差が積み重なりやすいですが、陰解法では将来の値を辻褄を合わせるように値を決めていくので、誤差が積み重なりにくいくらいの理解でとりあえずはいいかと思います。

 

 

行列で表現するとどうなるか、書いてみます。

$$\begin{align}
& \left[ \begin{array}{cccccc} 1+2r & -r &  &  &  & 0 \\ -r & 1+2r & -r &  &  & \\  & \ddots & \ddots & \ddots &  &  \\  &  & -r & 1+2r & -r & \\ 0 &  &  &  & -r & 1+2r \end{array} \right] \hspace{ 3pt } \left[ \begin{array}{c} u_{2}^{n+1} \\ u_{3}^{n+1} \\ \vdots \\ u_{J-2}^{n+1} \\ u_{J-1}^{n+1} \end{array} \right] \hspace{ 3pt }=\hspace{ 3pt } \left[ \begin{array}{c} u_{2}^{n} \\ u_{3}^{n} \\ \vdots \\ u_{J-2}^{n} \\ u_{J-1}^{n} \end{array} \right]
\end{align}$$

あとは、この行列の逆行列を求めてやれば、次のステップの値を求めることができます。

---

解法

今回は逆行列を解くという作業が入ります。これは一度やった、3項方程式なのでトーマス法で解けるのですが、またトーマス法でやってもつまらないので、ヤコビの反復法で解きます。

トーマス法は特殊な場合にしか使えませんが、ヤコビの反復法は基本どんな行列にも使えます。

内容を説明していきます。まず、逆行列を解く方法は直接法と反復法と言う方法に分かれます。

トーマス法は直接法に該当します。つまり、そのまま数学的に厳密な逆行列を求めてといてやろうという方法です。これは行列が複雑になってくると計算量がバカにならなくなります。

もう一つ、反復法というものがあります。

簡単な説明はこちらの記事をみてください。そして、反復法の一番基礎となるのが、ヤコビの反復法です。やってみましょう。

一般的な形にすると、

$$\begin{align}
& \left[ \begin{array}{cccc} a_{11} & a_{12} & \cdots & a_{1n} \\ a_{21} & a_{22} & \cdots & a_{2n} \\  \vdots & \vdots & \ddots & \vdots \\ a_{n1} & a_{n2} & \cdots & a_{nn} \end{array} \right] \hspace{ 3pt } \left[ \begin{array}{c} x_{1} \\ x_{2} \\ \vdots \\ x_{n} \end{array} \right] \hspace{ 3pt }=\hspace{ 3pt } \left[ \begin{array}{c} b_{1} \\ b_{2} \\ \vdots \\ b_{n} \end{array} \right]
\end{align}$$

で、l回目の反復に置ける回の推定値を、

$$\begin{align}
& ^{l} \mathbf{x} \hspace{ 3pt }=\hspace{ 3pt } \left[ \begin{array}{c} b_{1} \\ b_{2} \\ \vdots \\ b_{n} \end{array} \right]
\end{align}$$

とします。この時、次のステップの推定値\(^{l+1} \mathbf{x}\)を、

$$
\begin{align}
& ^{l+1} x_1 = \frac{1}{a_{11}} \Big( b_1 – a_{12} {}^{l}x_{2} – a_{13} {}^{l} x_3 \cdots – a_{1n} {}^{l} x_{n} \Big)\\\\
& ^{l+1} x_2 = \frac{1}{a_{22}} \Big( b_2 – a_{21} {}^{l}x_{1} – a_{23} {}^{l} x_3 \cdots – a_{2n} {}^{l} x_{n} \Big)\\\\
& \hspace{10px} \vdots \\\\
& ^{l+1} x_1 = \frac{1}{a_{nn}} \Big( b_1 – a_{n1} {}^{l}x_{1} – a_{n2} {}^{l} x_2 \cdots – a_{n,n-1} {}^{l} x_{n-1} \Big)\\\\
\tag{5} \label{eq:scale-factor-5}
\end{align}
$$

として解きます。何をやっているのかというと、n個の連立方程式を一つづつ解いていきます。i番目の方程式を解くときは、i番目の\({}^{l+1} x_i\)のみを未知数、あとは既知の値として解きます。

$$
\begin{align}
& a_{i1} {}^{l} x_1 +a_{i2} {}^{l} x_2 + a_{ii-1} {}^{l} x_{i-1} + a_{ii} {}^{l+1} x_i + a_{ii+1} {}^{l} x_{i+1} = b_i \\\\
& ^{l+1} x_i = \frac{1}{a_{ii}} \Big( b_1 – a_{n1} {}^{l}x_{1} – a_{n2} {}^{l} x_2 \cdots – a_{in} {}^{l} x_{n} \Big) \\\\
\tag{6} \label{eq:scale-factor-6}
\end{align}
$$

実にシンプルです。これでなんで解けるのでしょうか？イメージでいくと、これは擬似的に逆行列を解いているようなものです。全てのxの調律を満たしているのが完璧な解だとして、お互いへの影響を考えずに、まずは固定して一つのxについて調整して、それぞれやっていって、それを繰り返して行けば、真の解にたどり着くのじゃないかという発送です。

そして、繰り返すうちに値が変わらなくなってきたら収束したとみなして計算を終わります。絶対誤差を用いて収束を判定するには、

$$
\begin{align}
& \sum_{ i = 1 }^{ n } | {}^{l+1}x_{i} –  {}^{l}x_{i} | \leq \epsilon \\\\
\tag{7} \label{eq:scale-factor-7}
\end{align}
$$

相対誤差を用いた場合は、

$$
\begin{align}
& \sum_{ i = 1 }^{ n } | \frac{{}^{l+1}x_{i} –  {}^{l}x_{i}}{{}^{l+1}x_{i}} | \leq \epsilon \\\\
\tag{8} \label{eq:scale-factor-8}
\end{align}
$$

として、これを満たしたとき、計算を終わります。\(\epsilon\)は自分で設定します。計算量は、\(O(ln^2)\)なので、反復回数lをいかに減らすかが重要です。

なんとなく分かるでしょうか？完全法で逆行列を解くというのは、一気に全ての調律を満たしているものを見つけて、一気に解を求めますが、反復法はこうやって徐々に求めていきます。

また、反復法は収束するかどうかが保証されていません。収束条件をみてみます。
まず、反復式を以下のように書き直します。

$$
\begin{align}
& {}^{l+1} \mathbf{x} = \mathbf{H} \ {}^{l} \mathbf{x} + \mathbf{z}
\tag{9} \label{eq:scale-factor-9}
\end{align}
$$

また、反復回数が十分大きい数Lの時、以下の式が成立します。

$$
\begin{align}
& {}^{L} \mathbf{x} = \mathbf{H} \ {}^{L} \mathbf{x} + \mathbf{z}
\tag{10} \label{eq:scale-factor-10}
\end{align}
$$

この2式より、以下の式が成立します。

$$
\begin{align}
& {}^{l+1} \mathbf{x} – \mathbf{H} \ {}^{l} \mathbf{x} = {}^{L} \mathbf{x} – \mathbf{H} \ {}^{L} \mathbf{x}  \\\\
& \Leftrightarrow {}^{L} \mathbf{x} –  {}^{l+1} \mathbf{x} = \mathbf{H} ( {}^{L} \mathbf{x} – {}^{l} \mathbf{x}) \\\\
\tag{11} \label{eq:scale-factor-11}
\end{align}
$$

ここで、\({}^{L} \mathbf{x} –  {}^{l+1} \mathbf{x} \)は、l+1ステップでの真の解からの誤差、\({}^{L} \mathbf{x} –  {}^{l} \mathbf{x} \)はlステップでの誤差である。つまり、反復を繰り返すごとに、誤差が\(\mathbf{H}\)倍になっていく。そのため、解が収束するためには、

$$
\begin{align}
&\| \mathbf{H} \| < 1  \\\\
\tag{12} \label{eq:scale-factor-12}
\end{align}
$$

を満たす必要があります。これは、行列1ノルムと考えてもらえばいいです。\eqref{eq:scale-factor-5}式の反復式をみて、

$$
\begin{align}
\frac{\sum_{j \neq i} | a_{ij} |  }{| a_{ii} |} < 1 \\\\
\tag{13} \label{eq:scale-factor-13}
\end{align}
$$

この式より，係数行列の各行において対角要素の絶対値が非対角要素の絶対値の和より大きい場合に収束します。このような条件を対角優位(diagonal dominant)と呼びます。

解くときはこのままやればいいのですが、今回の問題の場合中身がどうなっているかをみてみます。

$$
\begin{align}
& a_{ii-1} = -r, \ a_{ii} = 1+2r, \ a_{ii+1} = -r , 他0より、\\\\
& ^{l+1} u_{i}^{n+1} = \frac{1}{1+2r} \Big( u_{i}^{n} + r \ {}^{l}u_{i-1}^{n+1} + r \ {}^{l} u_{i+1}^{n+1} \Big) \quad (i=3,…,n-2) \\\\
\tag{14} \label{eq:scale-factor-14}
\end{align}
$$

となってます。そして、初期解は一個前の時間のステップ\(u^{n-1}\)からそのまま持ってこればいいです。i=1,nは境界条件で決まり、i=2,n-1は一個\(r {}^{l}u_{i-1}^{n+1}\)のどっちかの項が消えた形になります。

このような形で解いてもいいですし、汎用的にも書けるので汎用的にプログラムを書いてもいいです。今回は汎用的に書いてみます。

ヤコビの反復法の兄弟分に、ガウス・ザイデル法とSOR法というのもあるのですが、だいたい同じです。どっかでやります。

---

プログラム

ではプログラムを実装しましょう。言語はとりあえずC++でやります。

以下のプログラムは、冒頭でも紹介した「流体解析の基礎」からデバッグ・改変して使っています。

今回はだいぶC++っぽく書きました。前の陽解法と比べてみてください。本当は、こういう疎行列を0ばっか埋めて作るのはあまりメモリ上よくないのですが、理解しやすいのでこうします。

そして、初期条件に、

\begin{eqnarray} f(x)  = \begin{cases} x & (0 \leq x \leq 0.5 ) \\ 1-x & (0.5 \leq x \leq 1 ) \end{cases} \end{eqnarray}

という形をとりあえず使います。

以下のファイルをHEAT_I.cppとでもして保存しましょう。

 

まあ前回と全く同じプロセスですが、計算を回してみます。

で実行できるかと思います。-g オプションはデバッグ用の情報を追加で記載する設定です。Macだと、

とかでデバッカーを使えます。普通gdbなのですが、Macがサポートしてないので、lldbです。

実行すると、animeというフォルダが作成されて、そこにファイルが書き出されると思います。

前回作った、anime.gpを使ってアニメーションも書いてみてください。

全く前回と同じですが、一応アニメーションのプログラムを。今回はC++じゃなくて、Cverを貼っときます。
以下のプログラムをheat_animetion_make.Cとでもして保存してください。

こいつで、anime.gpといった奴を作り、それを使ってgnuplotで簡単にアニメーションを書きます。

 

陰解法の方が精度が高いので、メッシュ数とかdtとか色々いじってみると面白いです。

---

[前の記事へ]　 [非圧縮性流体の目次へ]　 [次の記事へ]

---