差分法で2階常微分方程式を解く まずはトーマス法
まずは簡単なやつから
プログラムの計算では、以下の書籍の内容やプログラムを参考にして作ってます。プログラムは、一部改修・デバッグして使用しています。超いい本なので、ぜひ読んでください。
プログラムは朝倉書店のホームベーじからダウンロードできます。書籍名で検索してみてください。
コメント書いたり修正したプログラムも、githubにアップしておきます。
Cも修正してないプログラムは動かないこともあり、fortranも77なので、暇な時にいずれ修正する予定です。
まずは、シンプルな問題を解いていきましょう。
$$
\begin{align}
&\frac{ d^2 u }{ dx^2} + u + x = 0 \qquad (0 < x < 1) \\\\
&u(0) = 0, u(1) = 0
\tag{1} \label{eq:scale-factor-1}
\end{align}
$$
を考えます。この問題は厳密解を持ちます。
$$
\begin{align}
&u(x) = \frac{\sin x}{\sin 1} – x
\tag{2} \label{eq:scale-factor-2}
\end{align}
$$
差分法で解く場合、まずはメッシュを切る、つまり格子点に分割します。
念の為説明すると、格子点上だけが値を持つようにして、空間を離散的にします。値は滑らかな値を持たなくなり、飛び飛びの場所の値しかわかりません。もし、間の値が知りたければ補完します。
こうすると、テイラー展開等で解け、近似解を得ることができます。
こんな感じで、N個の格子、N+1個格子点に離散化します。格子点の番号を左から0,1,…,Nとし、\(x\)と\(u\)の値にも同じように\(x_0,x_1,..,\)、\(u_0,u_1,..,\)、と名前をつけていきます。代表点を、\(x_j, u_j\)とします。
今回は簡単のため、格子幅を等間隔で\(\Delta x\)とします。
$$
\begin{align}
&\Delta x = \frac{1}{N}, \quad x_j = j \Delta x \qquad (j=0,1,…,N)
\tag{3} \label{eq:scale-factor-3}
\end{align}
$$
次に、微分方程式\eqref{eq:scale-factor-1}を差分近似します。
$$
\begin{align}
&\frac{u_{j+1} -2u_j + u_{j-1}}{\Delta x^2} + u_j + x_j = 0 \\\\
& u_{j-1} + (\Delta x^2 – 2)u_j + u_{j+1} = \ -j \Delta x^3
\tag{4} \label{eq:scale-factor-4}
\end{align}
$$
となります。ここまではいいでしょう。境界条件の、\(u_0=u_N=0\)を考慮して、行列として書くと、
$$\begin{align}
& \left[ \begin{array}{cccccc} \Delta x^2 -2 & 1 & & & & 0 \\ 1 & \Delta x^2 -2 & 1 & & & \\ & \ddots & \ddots & \ddots & & \\ & & 1 & \Delta x^2-2 & 1 & \\ 0 & & & & 1 & \Delta x^2 -2 \end{array} \right] \hspace{ 3pt } \left[ \begin{array}{c} u_{1} \\ u_{2} \\ \vdots \\ u_{N-2} \\ u_{N-1} \end{array} \right] \hspace{ 3pt }=\hspace{ 3pt } \left[ \begin{array}{c} -\Delta x^3 \\ – 2 \Delta x^3 \\ \vdots \\ -(N-2)\Delta x^3 \\ -(N-1)\Delta x^3 \end{array} \right]
\end{align}$$
できました!あとは、この連立方程式を解けば終わりです。
行列を書くとき、境界条件も入れてもいいですが、プログラム上では境界条件は別に書くことが多いので、基本的に行列でも書かないことが多いです。境界条件を満たすように行列の中に入れて別に書いても大丈夫です。
トーマス法
今回の連立方程式は3項方程式と呼ばれるもので、トーマス法と呼ばれる方法で解くことができます。
3項方程式を一般的に書いてみます。
$$
\begin{align}
b_1x_1 + &c_1x_2&=d_1 \\\\
a_2x_1 + &b_2x_2 + c_2x_3 &=d_2 \\\\
&a_3x_2 + b_3x_3 + c_3x_4 &=d_3 \\\\
&\hspace{60pt}\vdots & \quad \\\\
&\hspace{50pt}a_{M-1}x_{M-2} + b_{M-1}x_{M-1} + c_{M-1}x_M &=d_{M-1} \\\\
&\hspace{100pt}a_{M}x_{M-1} + b_{M}x_{M} \hspace{10pt} &=d_{M} \\\\
\tag{5} \label{eq:scale-factor-5}
\end{align}
$$
\(x_1\)と\(u_1\)対応していて、あとは同じものだとわかります。
それを、こんな感じに順番に消去していって最後まで行きついてから、また逆に辿っていくという解法で解けます。
一番目の式より、
$$
\begin{align}
&x_1 = \frac{d_1-c_1x_2}{b_1} = \frac{s_1-c_1x_2}{g_1} \tag{6} \label{eq:scale-factor-6}
\end{align}
$$
ただし、\(g_1=b_1, s_1=d_1 \)
とおく。こんな感じで\(x_1\)を消去します。
2番目の式の\(x_1\)のところにこれを代入して、
$$
\begin{align}
&x_2 = \frac{s_2-c_2 x_3}{g_2} \tag{7} \label{eq:scale-factor-7}
\end{align}
$$
ただし、$$\begin{align}g_2=b_2 – \ \frac{a_2 c_1}{g_1} , s_2=d_2 – \ \frac{a_2 s_1}{g_1} \end{align}$$
3番目の式の\(x_2\)のところにこれを代入して、
$$
\begin{align}
x_3 = \frac{s_3-c_3x_4}{g_3} \tag{8} \label{eq:scale-factor-8}
\end{align}
$$
$$\begin{align}g_3=b_3 – \ \frac{a_3 c_2}{g_2} , s_2=d_3 – \ \frac{a_3 s_2}{g_2}\end{align}$$
置き換えで常に同じ形にするのがポイントです。i番目は、
$$
\begin{align}
x_i = \frac{s_i – c_i x_{i+1}}{g_i} \tag{9} \label{eq:scale-factor-9}
\end{align}
$$
$$\begin{align}g_i=b_i – \ \frac{a_i c_{i-1}}{g_{i-1}} , s_i=d_i – \ \frac{a_i s_{i-1}}{g_{i-1}}\end{align}$$
となります。この式はi=2,…Mで成り立ち、i=Mの時は、\(C_M\)がない為、\(C_M\)の部分を0とし、ここで\(x_M\)が求まります。
$$
\begin{align}
x_M = \frac{s_M}{g_M} \tag{10} \label{eq:scale-factor-10}
\end{align}
$$
これで\(x_M\)求まり、\eqref{eq:scale-factor-9}式から、i=M-1の時を考え、\(x_{M-1}\)が求まります。
このようにして順々に求めていくことができます。
プログラム
ではプログラムを実装しましょう。僕の好きな数値解析ソフトOpenFOAMがC++なのと、地味に好きな数値解析専用言語freefem++という言語もC++ベースなので、C++で書きますが、githubの方にはCとfortranとpythonで書いたものを上げようかと思います。元々のCからC++仕様に書き換えているので、確実に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 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 |
#include <iostream> #include <string> #include <vector> #include <fstream> #include <cmath> using namespace std; // ベクトルの参照渡し void thomas(int il,int iu,vector<double>& a,vector<double>& b,vector<double>& c,vector<double>& d); int main(void){ float h,x,u,err; int n,i; string filename = "output.dat"; ofstream writing_file; writing_file.open(filename, ios::out); n = 20; // mesh number // 要素数nで、0.0に初期化した配列を作成するが、なぜか2から始めるため一個多く確保する。 vector<double> a(n+1,0.0),b(n+1,0.0),c(n+1,0.0),d(n+1,0.0); h=(double)1./n; for(i=2;i<=n;i++){ a[i]=1.0; b[i]=h*h-2.0; c[i]=1.0; d[i]=(double)-(i-1)*h*h*h; } thomas(2,n, a,b,c,d); for(i=2;i<=n;i++){ x=(double)h*(i-1); u=sin(x)/sin(1.0)-x; err=(d[i]-u)/u*100.0; // 出力 writing_file <<"\t"<<scientific<<x <<"\t"<<scientific<<d[i] <<"\t"<<scientific<<u <<"\t"<<scientific<<err<<endl; } writing_file.close(); return 0; } void thomas(int il,int iu, vector<double>& a,vector<double>& b,vector<double>& c,vector<double>& d) { int ip,i,j; float r; ip=il+1; // 最初はg=b, s=dなのでそのまま使用。 for(i=ip;i<=iu;i++){ r=c[i]/b[i-1]; b[i]=b[i]-r*a[i-1]; // gを算出。最初がg=bなのでbをそのままg d[i]=d[i]-r*d[i-1]; // sを算出 } d[iu]=d[iu]/b[iu]; for(i=ip;i<=iu;i++){ j=iu-i+il; d[j]=(d[j]-a[j]*d[j+1])/b[j]; } return; } |
解説です。まず、解説ファイルが付いているので、それを読みます。
このプログラムは常微分方程式の境界値問題を解くプログラムである。nは格子数である。最初のループにおいて三項対角方程式の係数を計算した後、thomas でトーマス法の関数を呼んでいる。次のループは結果の表示のためのものであり、厳密解と比較している。
プログラムの仕様
配列 a —- 差分方程式のUi-1 (i=2,n)の係数を記憶する配列
配列 b —- 差分方程式のUi (i=2,n)の係数を記憶する配列
配列 c —- 差分方程式のUi+1 (i=2,n)の係数を記憶する配列
配列 d —- 差分方程式の右辺の値を記憶する配列、サブルーチン
THOMASを呼んだあとでは近似解が記憶される。
変数 n —- 格子数(格子点数はn+1になる)
変数 h —- 格子幅
変数 x —- 格子点のX座標
変数 u —- Xでの厳密解
変数 err —- 厳密解との相対誤差(%)
といった感じです。
一応C++をあまり書いたことのない人向けに少しだけ解説を入れます。
vectorは配列のクラスです。動的に配列を確保できますが、最初に要素数を定義してしまっています。これは初期化のためで、最初から配列の要素にアクセスできるようにするためです。mx+1個にしているのは、C言語の配列は0から始まるため、u[0]が一個目の値で、u[mx+1]がmx個目の値なのですが、分かりやすくするためこうしてます。慣れたら、0から始めます。
C++書いたことない人は、まずはC言語の苦しんで覚えるC言語でも軽くみて、一週間で身につくC++言語の基本でも流し読みすれば十分です。
BVP.cppというファイルに保存して、実行してみます。
1 2 |
g++ BVP.cpp ./a.out |
これでoutput.datというファイルが書き出されると思います。
gnuplotで可視化するために、output.gpというファイルを作ります。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 |
set pm3d map set contour set ticslevel 0 unset colorbox set cntrparam levels incremental -1,0.01,0 set xr[1:64] set yr[1:64] set palette defined (-1 "white", 0 "white") set cbrange[-1:0] set term gif animate set output "2D-cavity_unsteady-anime.gif" n0=100 n1=29900 dn=100 load "anime.plt" |
これでgnuplotで可視化してみます。
1 |
gnuplot output.gp |
グラフが書かれるかと思います。
この微分方程式が本質的にどのようなものかというと、まず時間項がありません。つまりは、定常解析と呼ばれるものです。実際に、流体でも定常解析をする時は時間項を無くしてときます。すると、ある解に落ち着くはずです。
対象が熱であれば、熱が広がり切って安定した状態の空間的な分布を求めている、といった状況になります。
ディスカッション
コメント一覧
まだ、コメントがありません