差分法で2階常微分方程式を解く まずはトーマス法

2019年9月28日非圧縮性流体

まずは簡単なやつから

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

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

コメント書いたり修正したプログラムも、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はあると思います。

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

 

gnuplotで可視化するために、output.gpというファイルを作ります。

これでgnuplotで可視化してみます。

グラフが書かれるかと思います。

この微分方程式が本質的にどのようなものかというと、まず時間項がありません。つまりは、定常解析と呼ばれるものです。実際に、流体でも定常解析をする時は時間項を無くしてときます。すると、ある解に落ち着くはずです。

対象が熱であれば、熱が広がり切って安定した状態の空間的な分布を求めている、といった状況になります。

---

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

---