1次元熱伝導を陰解法で解く ヤコビの反復法
次は陰解法で
プログラムは、以下の書籍の内容やプログラムを参考にして作ってます。プログラムは、一部改修・デバッグして使用しています。超いい本なので、ぜひ読んでください。
プログラムは朝倉書店のホームベーじからダウンロードできます。書籍名で検索してみてください。
コメント書いたり修正したプログラムも、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項方程式なのでトーマス法で解けるのですが、またトーマス法でやってもつまらないので、ヤコビの反復法で解きます。
トーマス法は特殊な場合にしか使えませんが、ヤコビの反復法は基本どんな行列にも使えます。
内容を説明していきます。まず、逆行列を解く方法は直接法と反復法と言う方法に分かれます。
トーマス法は直接法に該当します。つまり、そのまま数学的に厳密な逆行列を求めてといてやろうという方法です。これは行列が複雑になってくると計算量がバカにならなくなります。
もう一つ、反復法というものがあります。
適当な初期解 x(0)から始めて,繰り返し計算によって真の解に収束させていく方法です。
詳しい説明は東京大学さんの資料が分かりやすいです。
http://nkl.cc.u-tokyo.ac.jp/13n/SolverIterative.pdf
http://nkl.cc.u-tokyo.ac.jp/13n/SolverDirect.pdf
簡単な説明はこちらの記事をみてください。そして、反復法の一番基礎となるのが、ヤコビの反復法です。やってみましょう。
一般的な形にすると、
$$\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}
$$
この式より,係数行列
解くときはこのままやればいいのですが、今回の問題の場合中身がどうなっているかをみてみます。
$$
\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とでもして保存しましょう。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 |
#include <iostream> #include <string> #include <vector> #include <cmath> #include <sys/stat.h> using namespace std; /*********************************************************************** * DIFFUSION EQUATION EULER IMPLICIT METHOD * *********************************************************************** */ // よりc++っぽいプログラムにします。 void output_anime(const vector<double>& uw, int step); int JacobiIteration(vector< vector<double> > &A, int n, int &max_iter, double &eps); int main(void) { int i, ih, k, mx, km, ii, kop, counter, end_of_cal, max_iter; double dt, dx, r, x, alpha, eps; /***** INPUT & CALCULATE PARAMETERS */ mx=20; km=250; // mx:メッシュ数. km:ステップ数 alpha = 1.; // 熱伝達率。そのそも時間も長さも単位を定義していないので、単位はあまり考えない。 dt=0.001; // Δt dx = 1./(double)( mx - 1 ); // メッシュ幅 r = alpha * dt/(dx*dx); ih = (mx + 1)/2; // 初期条件に使うパラメータ kop = 5; // 書き出すステップ間隔 counter = 1; // ステップ間隔で書き出すのたるいので使う変数 // 0から始める仕様で書きます。 vector<double> u(mx,0.0), uu(mx,0.0); // 値と次の値を入れる。 // A n×nの係数行列とn×1の定数項(b)を併せたn×(n+1)の拡大行列. vector< vector<double> > A(mx, vector<double>(mx+1, 0.0)); /***** INITIAL CONDITION */ for( i = 0; i < mx; i++ ){ x = (double)( i )/(double)( mx - 1 ); if( i < ih ){ u[i] = x; } else{ u[i] = 1. - x; } } /***** MAIN LOOP *****/ for( k = 1; k <= km; k++ ){ u[0] = 0.; u[mx-1] = 0.; if( (k%kop) == 1 ){ output_anime(u, counter); counter++; } // 係数の設定 // 本来は境界まで係数行列に入れなくてもいいが、色々めんどいため入れた。 // そのため、境界はさっき代入したのが、k+1の値でそのままになって欲しいので、対角成分にだけ1を入れる。 A[0][0] = 1; A[mx-1][mx-1] = 1; A[0][mx] = u[0]; // 右辺 A[mx-1][mx] = u[mx-1]; for(i=1; i<mx-1; i++){ A[i][i-1]= -r; A[i][i]= 1.+ 2.*r; A[i][i+1]= -r; A[i][mx]= u[i]; // 右辺 } max_iter = 1000; eps = 0.0001; // ここで陰解法でとく。 JacobiIteration(A, mx, max_iter, eps); // 解を入れる。Aの一番右が解になっている。 for(i=0; i<mx; i++){ u[i] = A[i][mx]; } } } /*! * ヤコビ反復法(Jacobi iterative method) * - 解が収束するのは * ・対角有利(diagonal dominant, 対角要素の絶対値>その行の他の要素の絶対値の和) * ・係数行列が対称(symmetric)かつ正定(positive definite) * のどちらかの場合 * @param[inout] A n×nの係数行列とn×1の定数項(b)を併せたn×(n+1)の拡大行列.n+1列目に解が入る. * @param[in] n n元連立一次方程式 * @param[inout] max_iter 最大反復数(反復終了後,実際の反復数を返す) * @param[inout] eps 許容誤差(反復終了後,実際の誤差を返す) * @return 1:成功,0:失敗 */ int JacobiIteration(vector< vector<double> > &A, int n, int &max_iter, double &eps) { vector<double> x(n, 0.0); // 初期値はすべて0とする vector<double> y(n, 0.0); double e = 0.0; int k; for(k = 0; k < max_iter; ++k){ // 現在の値を代入して,次の解候補を計算 for(int i = 0; i < n; ++i){ y[i] = A[i][n]; for(int j = 0; j < n; ++j){ //y[i] -= (j != i ? A[i][j]*x[j] : 0.0); こうやっても書ける if (i != j){ y[i] -= A[i][j]*x[j]; } } y[i] /= A[i][i]; } // 収束判定 e = 0.0; for(int i = 0; i < n; ++i){ e += fabs(y[i]-x[i]); // 絶対誤差の場合 //e += fabs((y[i]-x[i])/y[i]); // 相対誤差の場合 } if(e <= eps){ break; } swap(x, y); } max_iter = k; eps = e; for(int i = 0; i < n; ++i){ A[i][n] = y[i]; } return 1; } void output_anime(const vector<double>& uw ,int step) { int i; char out_dir[256] = "./anime"; char out_path[256]; sprintf(out_path,"./anime/temp_%05d.dat",step); mkdir(out_dir, 0777); // sys.stat.hがないと機能しない。その場合は手動でanimeのフォルダを作ってコメントアウトする。 FILE *file; file = fopen(out_path, "w"); for(i=0; i<(int)uw.size(); i++){ fprintf(file, "\t%d\t%f\n",i,uw[i]); } fclose(file); } |
まあ、だいたいコメントに書いておきました。
今回は2次元配列を使っています。2次元配列に、行列の係数と右辺を格納しています。A[i][j]は(0から数えて)i行j列目です。
そして、その2次元配列をヤコビの反復法で解きます。式と色々見比べてみてください。
Cでトーマス法で書いたのもgibhubに上げておいたので、生成物を見比べて見ると面白いと思います。max_iter=100とかにすると、かなり雑で誤差がだいぶ出てきます。
まあ前回と全く同じプロセスですが、計算を回してみます。
1 2 |
g++ HEAT_I.C -g ./a.out |
で実行できるかと思います。-g オプションはデバッグ用の情報を追加で記載する設定です。Macだと、
1 |
lldb ./a.out |
とかでデバッカーを使えます。普通gdbなのですが、Macがサポートしてないので、lldbです。
実行すると、animeというフォルダが作成されて、そこにファイルが書き出されると思います。
前回作った、anime.gpを使ってアニメーションも書いてみてください。
全く前回と同じですが、一応アニメーションのプログラムを。今回はC++じゃなくて、Cverを貼っときます。
以下のプログラムをheat_animetion_make.Cとでもして保存してください。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 |
#include <stdio.h> #include <math.h> #include <stdlib.h> #include <string.h> int main(void) { int mx = 20; // 描画する幅。メッシュ数に合わせる。 int step, step_end = 50; // 書き出したファイルのステップ分 char base[256]; char out_path[256]; sprintf(base,"./anime/temp_"); sprintf(out_path, "anime.gp"); FILE *file; file = fopen(out_path, "w"); fprintf(file,"#\n"); fprintf(file,"# gnuplot script generated by heat_animation_make.C\n"); fprintf(file,"#\n"); fprintf(file,"set xlabel 'j' # x-axis\n"); fprintf(file,"set ylabel 'temperature' # y-axis\n"); fprintf(file,"set xrange[0:%d] # j-grid min & max\n", mx+1); fprintf(file,"set yrange[0.0:0.5] # temperature min & max\n"); for(step=1; step<step_end; step++){ fprintf(file,"plot './anime/temp_%05d.dat' w lp\n", step); if(step==1){ fprintf(file,"pause 5\n"); } else{ fprintf(file,"pause 1\n"); } } fclose(file); } |
こいつで、anime.gpといった奴を作り、それを使ってgnuplotで簡単にアニメーションを書きます。
1 2 3 |
gcc heat_animation_make.C ./a.out gnuplot anime.gp |
陰解法の方が精度が高いので、メッシュ数とかdtとか色々いじってみると面白いです。
ディスカッション
コメント一覧
まだ、コメントがありません