右上の添字は時間を、右下の添字は空間を表すのが慣例なので、そのようにし、時間項をオイラー陰解法で、2階微分を2次中心差分で離散化します。\label{eq:scale-factor-2}\end{align}$$境界条件と初期条件も前回と同様に考えます。境界条件は\(u^{n+1}\)について考えます。\label{eq:scale-factor-4}\end{align}$$となります。以下のような関係性になっています。陽解法と大きく違うのは、単純には次のステップの値を求められない関係性になっていると言うことです。陽解法で解けるのに、なぜわざわざこんなことをするのかと思うかもしれませんが、陰解法の方が圧倒的に正しい解を求めやすいからです。具体的に言うと、よく数値解析の場合は発散と言う現象が起こります。明らかにありえない解や解が振動してしまう等のことがおきます。t\)が大きすぎたり、はたまた微分方程式の差分化する時の精度が問題だったり、色々です。解く問題によっても全然変わります。熱伝達率が大きかったり熱の差が大きい問題を解く時は、空間も時間も細かくせねばなりませんし、全体的に温度があまり変わらないなら大きくしても構いません。細かくすればいいと言う問題ではなく、計算時間とトレードオフなので、適切な設定をすることが大事です。x\)を大きくしてみてください。その発散と言う現象が陰解法の方が圧倒的に起こりにくいです。数学的に証明もできるのですがややこしくなるので簡単な説明に留めます。陽解法では過去の値から素直に求めているので、誤差が積み重なりやすいですが、陰解法では将来の値を辻褄を合わせるように値を決めていくので、誤差が積み重なりにくいくらいの理解でとりあえずはいいかと思います。行列で表現するとどうなるか、書いてみます。\right]\end{align}$$あとは、この行列の逆行列を求めてやれば、次のステップの値を求めることができます。今回は逆行列を解くという作業が入ります。これは一度やった、3項方程式なのでトーマス法で解けるのですが、またトーマス法でやってもつまらないので、ヤコビの反復法で解きます。トーマス法は特殊な場合にしか使えませんが、ヤコビの反復法は基本どんな行列にも使えます。内容を説明していきます。逆行列を解く方法は直接法と反復法と言う方法に分かれます。トーマス法は直接法に該当します。そのまま数学的に厳密な逆行列を求めてといてやろうという方法です。これは行列が複雑になってくると計算量がバカにならなくなります。一つ、反復法というものがあります。x(0)から始めて,繰り返し計算によって真の解に収束させていく方法です。詳しい説明は東京大学さんの資料が分かりやすいです。http://nkl.cc.u-tokyo.ac.jp/13n/SolverIterative.pdfhttp://nkl.cc.u-tokyo.ac.jp/13n/SolverDirect.pdf簡単な説明はこちらの記事をみてください。反復法の一番基礎となるのが、ヤコビの反復法です。\right]\end{align}$$とします。\label{eq:scale-factor-5}\end{align}$$として解きます。何をやっているのかというと、n個の連立方程式を一つづつ解いていきます。x_i\)のみを未知数、あとは既知の値として解きます。\label{eq:scale-factor-6}\end{align}$$実にシンプルです。これでなんで解けるのでしょうか?イメージでいくと、これは擬似的に逆行列を解いているようなものです。全てのxの調律を満たしているのが完璧な解だとして、お互いへの影響を考えずに、まずは固定して一つのxについて調整して、それぞれやっていって、それを繰り返して行けば、真の解にたどり着くのじゃないかという発送です。繰り返すうちに値が変わらなくなってきたら収束したとみなして計算を終わります。\label{eq:scale-factor-8}\end{align}$$として、これを満たしたとき、計算を終わります。\(\epsilon\)は自分で設定します。計算量は、\(O(ln^2)\)なので、反復回数lをいかに減らすかが重要です。なんとなく分かるでしょうか?完全法で逆行列を解くというのは、一気に全ての調律を満たしているものを見つけて、一気に解を求めますが、反復法はこうやって徐々に求めていきます。反復法は収束するかどうかが保証されていません。収束条件をみてみます。反復式を以下のように書き直します。\label{eq:scale-factor-9}\end{align}$$また、反復回数が十分大きい数Lの時、以下の式が成立します。\label{eq:scale-factor-10}\end{align}$$この2式より、以下の式が成立します。\)はlステップでの誤差である。反復を繰り返すごとに、誤差が\(\mathbf{H}\)倍になっていく。\label{eq:scale-factor-12}\end{align}$$を満たす必要があります。これは、行列1ノルムと考えてもらえばいいです。\label{eq:scale-factor-13}\end{align}$$この式より,係数行列の各行において対角要素の絶対値が非対角要素の絶対値の和より大きい場合に収束します。dominant)と呼びます。解くときはこのままやればいいのですが、今回の問題の場合中身がどうなっているかをみてみます。\label{eq:scale-factor-14}\end{align}$$となってます。初期解は一個前の時間のステップ\(u^{n-1}\)からそのまま持ってこればいいです。{}^{l}u_{i-1}^{n+1}\)のどっちかの項が消えた形になります。このような形で解いてもいいですし、汎用的にも書けるので汎用的にプログラムを書いてもいいです。今回は汎用的に書いてみます。ヤコビの反復法の兄弟分に、ガウス・ザイデル法とSOR法というのもあるのですが、だいたい同じです。プログラムではプログラムを実装しましょう。言語はとりあえずC++でやります。以下のプログラムは、冒頭でも紹介した「流体解析の基礎」からデバッグ・改変して使っています。今回はだいぶC++っぽく書きました。前の陽解法と比べてみてください。本当は、こういう疎行列を0ばっか埋めて作るのはあまりメモリ上よくないのですが、理解しやすいのでこうします。\end{eqnarray}という形をとりあえず使います。以下のファイルをHEAT_I.cppとでもして保存しましょう。c++っぽいプログラムにします。そのそも時間も長さも単位を定義していないので、単位はあまり考えない。0から始める仕様で書きます。本来は境界まで係数行列に入れなくてもいいが、色々めんどいため入れた。そのため、境界はさっき代入したのが、k+1の値でそのままになって欲しいので、対角成分にだけ1を入れる。