技術計算：有限差分法
現在工事中です。
- 1. 非線型方程式とは
- 4.1 Newton法の基礎
【注】数式表示にはMathJaxを利用しています。IE8以下では表示が遅くなる可能性があります。FireFox などIE8以外のブラウザを利用下さい。
メッシュ分割
連立の常・偏微分方程式を解くことを考える。ここでは微分方程式を解きやすい形に変形するまでの全体を展望し、個々の項目は次章以下に詳述する。
時・空間をメッシュに分割
人手による分割方法には、等分割と非等分割（注１）とがある。解析空間を分割し、離散化(discretization)を行う。メッシュ点（分割点）上で、従属変数（濃度、温度など）を定義する。メッシュ点以外の点は補間することになる。
【注１】メッシュ間隔が小さい程、解析の精度が増加する。一方、メッシュ点が増大すると計算負荷も増大する。 解析空間で、従属変数が大きく変化する領域を細かく刻むことを行う。 等比級数的な刻み方とは、j+1番目の刻み幅 \( h_{j+1} \) が、その前のメッシュ幅 \( h_j \) と次の比例関係で表され、 \[ h_{j+1} = \alpha h_j \tag{1} \]
比例定数（等比係数）αが１より大きいとき、メッシュ幅は拡大方向、αが１より小さいとき、メッシュ幅は縮小方向。α=1のとき等間隔となる。 メッシュを刻む際に、隣り合うメッシュ幅が急変しないよう、αの値を適当に決める必要がある。 \( 0.5 \lt \alpha \lt 2 \) 程度に収まるよう刻むのがよい。
高階微分を差分化
微分方程式と境界・初期条件をともに、差分（前進差分、後退差分、片側差分、風上差分、中心差分など）で近似する。このとき近似精度・収束安定性（注２）に注意を払う必要がある。
差分化し、解くにあたり、陽的(Explicit)解法(Forward Euler法、Runge-Kutta法、Gear法など初期値問題を前進的に解く方法)と、陰的(Implicit)解法（Backward Euler、Implicit RKなど）の２つの解き方がある。差分化するときの考え方が少し違う。詳細は次章。
陽的解法では収束安定性に問題がある。また初期値問題には良いが、境界値特に出口境界指定のとき収束計算となり、メリット少ない。基本的にRK法などの陽解法は解が急変する領域近傍で物質・熱収支などを満足しないことが多い（例えば、モノマーの暴走反応をSimulateするとき、暴走後の天井温度が無意味な値となる）。
一方、陰解法を用いることにより、陽解法の欠陥・欠点を防ぐことが可能となる。また収束安定性（刻み幅による制限なし）、物質・熱収支の同時満足する解が得られる。ただし、計算機リソース（特にメモリ容量、CPU性能）が必要であるが、現状PC(Pentiumマシン、Widows XP以上、GBクラスのメモリ)では特にネックにならない。
【注２】一般に陽的解法ではメッシュの刻み方（メッシュ幅）により収束不安定となることがある。一方、陰的解法では刻み方に依存せず絶対安定である。陽解法のクーランの安定条件（後出）を満たす刻み方をする。
連立線形方程式化
通常化学反応を伴うとき、反応項があるゆえ、差分化しても非線形連立方程式となる。
非線形連立方程式を、Newton法などで線形化し、収束計算により解く。詳細は次回「非線形方程式の解法」を参照。
正方行列化、バンド行列化、三角行列(tridiagonal)など計算環境に依存するが、連立線形方程式系に帰着させる。
連立の線形方程式は、マルチCPUに特化した高速解法（ベクトル化、OpenMP)が利用できる。
線形化された連立方程式の解法
ガウス消去法、LU分解法、…など既存のライブラリ（高級言語に付属するライブラリIMSL、MKS、インターネット上に公開されているライブラリLAPACK、BLAS、ODEPACKなど、図書に記載されているコード）を利用可能で、コードの入手もネックにはならない。
デバッグや改良を行うためには、これらのソースコードを入手しておく方がよい。
差分化
ここでは微分方程式に含まれる微分項を差分化（離散化、discretization）する方法を解説する。
関数 \( f(x) \) を \( x=x_1 \) のまわりでTaylor展開するとは、次式となる。 \[ \begin{align*} f(x) & =f(x_1)+f'(x_1)(x-x_1)+ \frac1{2!}f''(x_1)(x-x_1)^2+ \cdots \cr & +\frac 1{n!}f^{(n)}(x_1)(x-x_1)^n+ \cdots \tag{2} \end{align*} \] Taylor展開を基礎とし、以下に示す差分化を行う。
前進、後退、中心差分
濃度の流れ方向（ｚ方向）の一階微分を考える。等分割とし、分割幅 \( \Delta z \)とする。入口側からメッシュ番号を付け、下流に向かうにつれメッシュ番号が増加する。
前進（Forward）差分は、Tayler展開の1次項までを下流側（風下側）とったもので、次式となる。 \[ \begin{align*} \frac {dC}{dz} = \frac 1{\Delta z}(C_{j+1}-C_j)+O(\Delta z) \tag{3} \end{align*} \]
後退(Backward)差分は、Taylor展開の1次項までを上流側（風上側）にとり、同じく1次精度となる。 \[ \begin{align*} \frac {dC}{dz} = \frac 1{\Delta z}(C_{j}-C_{j-1})+O(\Delta z) \tag{4} \end{align*} \]
中心(Center)差分は、前進差分と後退差分の算術平均をとったもの。Tayler展開後、二階微分項を消去。２次精度となる。 \[ \begin{align*} \frac {dC}{dz} & = \frac 12 \biggl\{ \frac 1{\Delta z}(C_{j}-C_{j-1}) + \frac 1{\Delta z}(C_{j+1}-C_j) \biggr\} + O(\Delta z^2) \cr & = \frac 1{\Delta 2z}(C_{j+1}-C_{j-1}) + O(\Delta z^2) \tag{5} \end{align*} \]
前進、後退、中心差分を一般化した一階差分式は、重みパラメータωとし、ωを0から１の間の値を適当に設定することで、前進、後退、中心差分となる。なお精度は刻み幅の1次から2次の間になる。 \[ \begin{align*} \frac {dC}{dz} & = \biggl\{ \frac {1-\omega}{\Delta z}(C_{j}-C_{j-1}) + \frac {\omega}{\Delta z}(C_{j+1}-C_j) \biggr\} + O(\Delta z \sim \Delta z^2) \tag{6} \end{align*} \]
二階微分も、Taylor展開した、表２の式から、一階微分項を消去し、次式が得られる。 \[ \begin{align*} \frac {d^2C}{dz^2} & = \frac 1{\Delta z^2}(C_{j+1}-2C_j+C_{j-1}) + O(\Delta z^2) \tag{7} \end{align*} \]
| \[ \begin{align*} u_{\pm 1} =u_0 \pm h u'_0+\frac 1{2!}h^2u''_0 \pm \frac 1{3!}h^3u'''_0 + \frac 1{4!}h^4u_0^{(4)}+\cdots \tag{8} \end{align*} \] |
| \[ \begin{align*} u_{\pm 2}=u_0 \pm 2h u'_0+\frac 1{2!}(2h)^2u''_0 \pm \frac 1{3!}(2h)^3u'''_0 + \frac 1{4!}(2h)^4u_0^{(4)}+\cdots \tag{9} \end{align*} \] |
参考として、不等分割でTaylor展開し、１階微分および２階微分の差分式を表３に示す。
| \[ \begin{align*} u'_0 = \frac 1{h_{-1}+h_0} \biggl\{ - \frac {h_0}{h_{-1}}(u_{-1}-u_0) + \frac {h_{-1}}{h_0}(u_1-u_0) \biggr\} + O(h^2) \tag{10} \end{align*} \] |
| \[ \begin{align*} u''_0 = \frac 2{h_{-1}+h_0} \biggl\{ - \frac 1{h_{-1}}(u_{-1}-u_0) + \frac 1{h_0}(u_1-u_0) \biggr\} \cr - \frac 13 (h_0-h_{-1})u'''_0+O(h^2) \tag{11} \end{align*} \] |
片側差分
境界条件を差分化するときに利用される。また、計算された離散点の値から境界での微係数を求めるときに利用される。たとえば多孔性触媒の有効係数を求めるとき、表面での拡散フラックス（濃度勾配）を必要とし、この算出に高次片側差分を利用する。
表４に左側（入口側）の片側差分近似を、表５に右側（出口側）の片側差分近似を示す。
| \[ \begin{align*} u'_0 & = \frac {-u_0+u_1}h + O(h) \tag{12} \\ u'_0 & = \frac {-3u_0+4u_1-u_2}{2h} + O(h^2) \tag{13} \\ u''_0 & = \frac {u_0-2u_1+u_2}{h^2} + O(h) \tag{14} \\ u''_0 & = \frac {2u_0-5u_1+4u_2-u_3}{h^2} + O(h^2) \tag{15} \\ \end{align*} \] |
| \[ \begin{align*} u'_0 & = \frac {-u_{-1}+u_0}h + O(h) \tag{16} \\ u'_0 & = \frac {u_{-2}-4u_{-1}+3u_0}{2h} + O(h^2) \tag{17} \\ u''_0 & = \frac {u_{-2}-2u_{-1}+u_0}{h^2} + O(h) \tag{18} \\ u''_0 & = \frac {-2u_{-3}+4u_{-2}-5u_{-1}+2u_0}{h^2} + O(h^2) \tag{19} \\ \end{align*} \] |
風上差分
次のBurgers方程式を考える。Burgers方程式とは、典型的な非定常・一次元で反応のない拡散方程式（エネルギー収支式、運動方程式）に相当する、二階偏微分方程式である。左辺は非定常項と移流項（通常係数ｃは流れの速度である）、左辺は拡散項である。 \[ \begin{align*} \frac {\partial u}{\partial t} + c \frac{\partial u}{\partial x} = \nu \frac {\partial ^2u}{\partial x^2} \tag{20} \end{align*} \]
移流項が拡散項に比べ大きくなると、解が振動する傾向にあり、陽的解法では収束不安定となることがある。これを防ぐため、風上差分が用いられる。(28)式の移流項を風上差分を用い差分化した式を表６に示す。 \[ \frac{u_0^{n+1}-u_0}{\Delta t} + \frac{c+|c|}{2} \frac{u_0-u_{-1}}{\Delta x} + \frac{c-|c|}{2} \frac{u_1-u_0}{\Delta x} = \nu \frac {u_{-1}-2u_0+u_1}{\Delta x^2} \tag{21} \]
流れ速度ｃが正のとき、左辺第２項が残り、第３項はゼロとなり消去される。すなわち 項の差分が有効になり、これは流れの上流側の差分をとったことになる。逆に、速度ｃが負のとき、第２項はゼロ、第３項がのこり 項が有効になり、やはり流れの上流側（upstream)の差分をとったことになる。
移流項をこのように風上差分とすることで、収束安定性を増大させることが可能となる。CFDなどの流体解析ではよく利用される手法であるが、運動方程式を解かない反応器解析の計算では使うことは少ない。計算の安定性を増大させるために、連立して解くことを行う（陰的解法を使う）ことが多い。
なお、(29)式で左辺第１項の時間微分は、Forward Euler法（オイラー前進差分法）で差分化しており、次の時間ステップを表す上付添字n+1が記述しているが、ほかの変数はすべて時間ステップｎのときの値を採用する（陽解法）。
なお、Courant数Cと格子Reynolds数Rとを次のように定義する。 \[ \begin{align*} C \equiv \frac {c \Delta t}{\Delta x},~~ R \equiv \frac {c \Delta x}{\nu} \tag{22} \end{align*} \]
このとき、陽解法の収束安定条件として、表７のCourantの安定条件（必要条件）と表８のハートの条件も満足する必要がある。 \[ \begin{align*} |C| \lt 1,~ & {\rm or} ~ \Delta t \lt \frac{\Delta x}{|c|} \tag{23} \\ \frac {2C}R \lt 1, ~ & {\rm or} ~ \Delta t \lt \frac{\Delta x^2}{2\nu} \tag{24} \end{align*} \] \[ \begin{align*} CR \lt 2, ~ {rm or} ~ \Delta t \lt \frac {2\nu}{c^2} \tag{25} \end{align*} \]
Crank-Nicholson法
上のBurgers方程式(28)式の右辺の拡散項を差分化するとき、Crank-Nicholson法では、時刻ｎのときの値と、n+1のときの値との算術平均を用い、次式で表す二階微分近似をする。 \[ \begin{align*} \frac{\partial^2 u}{\partial x^2} = \frac 12 \biggl\{ \frac {u_{j+1}^{n+1}-2u_{j}^{n+1}+u_{j-1}^{n+1}}{\Delta x^2} + \frac {u_{j+1}^{n}-2u_{j}^{n}+u_{j-1}^{n}}{\Delta x^2} \biggr\} + O(\Delta x^2) \tag{26} \end{align*} \]
従って、Crank-Nicholson法を用いたBurgers方程式の差分化式では、新しい時刻n+1の変数が式の両辺に現れ、連立の線形方程式となる。求めたい変数の時刻n+1の値を使い差分化するとき、解法を陰的解法、陰解法(Implicit)という。通常、連立方程式となる。逆に時刻ｎで既に求まっている値だけを用い、差分化するときを陽的解法、陽解法（Explicit)といい、見かけ上の計算は楽になる。上式右辺括弧内第１項がないことになる。ストレートな計算が可能。
通常のRunge-Kutta法、Forward Euler法などは陽解法の典型例。
反応を伴うBurgers方程式では、反応項をn+1での値で代表させるのか、ｎでの値で代表させるかにより、線形となったり、非線形となったりする。n+1での値を使い差分化するときやはり陰解法となる。
陰解法
Burgers方程式をBackward Euler法で差分化した式を、(35)式に示す。時間微分項は前進差分、移流項は時刻n+1の値を用いた中心差分、拡散項はn+1のときの中心差分で近似している。 \[ \begin{align*} \frac{u_0^{n+1}-u_0^n}{\Delta t} +c \frac {u_{1}^{n+1}-u_{-1}^{n+1}}{2 \Delta x} = \nu \frac {u_{1}^{n+1}-2u_{0}^{n+1}+u_{-1}^{n+1}}{\Delta x^2} + O(\Delta t,\Delta x^2) \tag{27} \end{align*} \]
これを変形して、 \[ \begin{align*} \biggl( - \frac{c}{2 \Delta x} -\frac{\nu}{\Delta x^2} \biggr)u_{-1}^{n+1} & + \biggl( \frac{1}{\Delta t} +\frac{2\nu}{\Delta x^2} \biggr)u_{0}^{n+1} \\ & + \biggl( \frac{c}{2 \Delta x} - \frac{\nu}{\Delta x^2} \biggr)u_{1}^{n+1} = \frac{u_0^n}{\Delta t} \tag{28} \end{align*} \]
と表され、前後各１点の値を使い表現できる。すべてのメッシュ点の値を得るには、(36)式を連立して解けばよい。その時の係数行列はTridiagonal Matrixとなり、Gaussの消去法で簡単に解くことができる。
Backward Euler法は陰解法の一つであるが、収束安定性のための必要十分条件はなく、無条件安定である。刻み幅などの制約がないため、極めて収束性に対しては有効である。ただし、式（アルゴリズム）が複雑になったり、メモリを一時に消費する（行列のサイズが大きくなるとサイズの2乗に効く）。
初期条件・境界条件の差分化
初期条件（注１）や、入出口境界の条件も差分化する必要がある。これら条件が微分形式で与えられるときには、片側差分を使い、メッシュ点で変数を満足するよう差分化する。これら条件が微分でなく、直接値のときはその変数のメッシュ点の値をダイレクトに指定すればよい。
【注１】初期条件とは、非定常解析における時刻t=0での変数値をいう。境界条件はより一般に解析空間座標系の境界での条件をいう。初期条件も広義には境界条件といえる。
は、ここで取り扱っています。
Literature Cited
- 引用文献
- ) 矢川：「流れと熱伝導の有限要素法入門」、培風館(1983).