第7-1回　陽解法と陰解法[python, julia, fortran]

• English (United States)
• 日本語

第7回では、二次元非定常拡散方程式に対し陰解法の有限体積法で解くことを目標とします。

そこで、まずは陰解法とは何ぞやということで説明します。

1. 陽解法と陰解法の違い

1.1. 陽解法とは何か？

今までの記事では陽解法と呼ばれる方法で離散化していました。これは前の時間の情報から未来の情報を求める方法で下のようなイメージとなります（黒がすでに分かっている情報、赤がこれから求める値）。

1.2. 陰解法とは何か？

一方で陰解法は未来の値を使い解を求める方法になります（黒がすでに分かっている情報、赤がこれから求める値 ）。

陰解法の具体的な解き方は次回の第7-2回で行います。

2. 陽解法と陰解法の安定性比較

今回は陽解法と陰解法の違い、メリット・デメリットについて、フォンノイマンの安定性解析を用いて比較します。

フォンノイマンの安定性解析は、安定的に計算を行うための一つの指標として用いられます。まずは下記の移流方程式に対して行ってみます。

$$\frac{\partial u}{\partial t}+c\frac{\partial u}{\partial x}=0$$

2.1. 陽解法の安定性

今までと同じ形である陽解法で行います。

時間に関して前進差分、空間に対して後退差分を行い整理します。

$$\frac{u^{n+1}_j-u^n_j}{\Delta t}+c\frac{u_j^n-u^n_{j-1}}{\Delta x}=0$$

$$⇔u^{n+1}_j=u^n_j-\nu(u_j^n-u_{j-1}^n)　\ \ \ \ \ \dots ①$$

ただし、\(\nu=c\frac{\Delta t}{\Delta x}\)です。



ここで、フォンノイマンの安定性解析で使われる下記の式を代入します。

$$u^n_j=g^n exp(i(j\theta))$$

\(exp(i(j\theta))\)は複素数平面において単位円を表すので、\(g^n\)が大きさとなります。時間は常に進むため、\(n\)は大きくなります。ここで\(|g|> 1\)のとき、\(g^n\)は無限に大きくなり発散します。

つまり、発散してはいけないので\(|g| \le 1\)が安定条件となります。

それでは代入し、整理していきます。

①より

$$ g^{n+1} exp(i(j\theta)) = g^n exp(i(j\theta))-\nu( g^n exp(i(j\theta)) – g^n exp(i(j-1)\theta) ) $$

$$⇔ g=1-\nu+\nu exp(-i\theta)$$

$$⇔ g=1-\nu+\nu cos\theta-i\nu sin\theta$$

$$⇔ |g|^2=(1-\nu+\nu cos\theta)^2+(\nu sin\theta)^2$$

$$⇔ |g|^2=1+2\nu ^2(1-cos\theta)-2\nu (1-cos\theta)$$

\(|g| \le 1\)より

$$ 1+2\nu ^2(1-cos\theta)-2\nu (1-cos\theta) \le 1$$

$$ ⇔ 2\nu ^2(1-cos\theta)-2\nu (1-cos\theta) \le 0$$

$$ ⇔ \nu \le 1$$

$$ ⇔ c\frac{\Delta t}{\Delta x} \le 1$$

これが 安定的に計算を行える条件になります。今までの記事はこれを満たすように設定していました。

何気なく空間項に対して、現在の時間nを使いましたが、これが陽解法になります。

2.2. 陰解法の安定性

次は空間項に対して、未来の時間n+1を使う陰解法で行ってみます。

$$\frac{u^{n+1}_j-u^n_j}{\Delta t}+c\frac{u_j^{n+1}-u^{n+1}_{j-1}}{\Delta x}=0$$

$$⇔u^{n+1}_j=u^n_j-\nu(u_j^{n+1}-u_{j-1}^{n+1})　\ \ \ \ \ \dots ②$$

$$ ⇔ g^{n+1} exp(i(j\theta)) = g^n exp(i(j\theta))-\nu( g^{n+1} exp(i(j\theta)) – g^{n+1} exp(i(j-1)\theta) ) $$

$$⇔ g=1-\nu g+\nu gcos\theta-i\nu gsin\theta$$

$$⇔g=\frac{1}{1+\nu-\nu cos \theta +i\nu sin\theta}$$

$$⇔g=\frac{1+\nu-\nu cos \theta -i\nu sin\theta}{(1+\nu-\nu cos\theta)^2+(\nu sin\theta)^2}$$

$$⇔|g|^2= \frac{(1+\nu-\nu cos \theta)^2+(\nu sin\theta)^2}{((1+\nu-\nu cos\theta)^2+(\nu sin\theta)^2)^2} $$

$$⇔|g|^2= \frac{1}{(1+\nu-\nu cos\theta)^2+(\nu sin\theta)^2} $$

\(|g| \le 1\)より

$$ \frac{1}{(1+\nu-\nu cos\theta)^2+(\nu sin\theta)^2} \le 1$$

$$ ⇔1 \le (1+\nu-\nu cos\theta)^2+(\nu sin\theta)^2 $$

$$⇔0 \le \nu ^2 (1-cos\theta)+\nu(1-cos\theta)$$

$$⇔0 \le \nu+1$$

\(\nu=c\frac{\Delta t}{\Delta x} \ge 0\)より、常に安定である。

この「常に安定」が陰解法の強みです。

3. まとめ

陽解法では空間の解像度をあげるため、\(\Delta x\)を10倍細かくしたら\(\Delta t\)も10倍細かくしなければいけません。
陰解法では常に安定のため、\(\Delta x\)を10倍細かくしても\(\Delta t\)を変更する必要がありません。

もちろん陰解法にもデメリットがあり「1step当たりの計算時間が長い」「時間刻み幅が荒いと精度がおちる」などありますが、現在は広く使われております。



次回は具体的な陰解法の解き方について行います。下記に拡散方程式の場合も書いておきますので、興味があれば見てみてください。

*コードは7-5回にあります。

by hide

補足.　拡散方程式の安定性解析

$$\frac{\partial u}{\partial t}=D \frac{\partial u}{\partial x^2}$$

陽解法の場合（時間前進：差分、空間：中心差分）

$$\frac{u^{n+1}_j-u^n_j}{\Delta t}=D \frac{u^n_{j+1} -2 u_j^n+u^n_{j-1}}{\Delta x^2}$$

$$⇔u^{n+1}_j=u^n_j+\nu( u_{j+1}^n -2 u_j^n+u_{j-1}^n)　\ \ \ \ \ \dots ①$$

ただし、\(\nu=D\frac{\Delta t}{\Delta x^2}\)です。


①より

$$ g^{n+1} exp(i(j\theta)) = g^n exp(i(j\theta))+\nu( g^n exp(i(j+1)\theta) -2 g^n exp(i(j\theta)) + g^n exp(i(j-1)\theta) ) $$

$$⇔g=1+\nu exp(i\theta)-2\nu +\nu exp(-i\theta)$$

$$⇔g=1-2\nu +2\nu cos\theta =1-2\nu (1-cos\theta )$$

$$⇔|g|^2=1-4\nu(1-cos\theta )+4\nu ^2(1-cos\theta)^2$$

\(|g| \le 1\)より

$$ ⇔-4\nu(1-cos\theta )+4\nu ^2(1-cos\theta)^2 \le 0 $$

$$ ⇔-1+\nu (1-cos\theta) \le 0$$

$$ ⇔nu \le \frac{1}{1-cos\theta } \le \frac{1}{2}$$

$$ ⇔D \frac{\Delta t}{\Delta x^2} \le \frac{1}{2}$$

拡散方程式はより厳しい条件になります。

陰解法の場合（時間前進：差分、空間：中心差分）

$$\frac{u^{n+1}_j-u^n_j}{\Delta t}=D \frac{u^{n+1}_{j+1} -2 u_j^{n+1}+u^{n+1}_{j-1}}{\Delta x^2}$$

$$⇔u^{n+1}_j- \nu (u^{n+1}_{j+1} -2 u_j^{n+1}+u^{n+1}_{j-1})=u^n_j　\ \ \ \ \ \dots ②$$

ただし、\(\nu=D\frac{\Delta t}{\Delta x^2}\)です。

②より

$$ g^{n+1} exp(i(j\theta))- \nu( g^{n+1} exp(i(j+1)\theta) -2 g^{n+1} exp(i(j\theta)) + g^{n+1} exp(i(j-1)\theta) ) = g^n exp(i(j\theta)) $$

$$⇔ g(1- \nu exp(i\theta) +2\nu – \nu exp(-i\theta) ) = 1 $$

$$⇔ g(1+2\nu – 2 \nu cos\theta) = 1 $$

$$⇔ g = \frac{1}{ 1+2\nu( 1- cos\theta ) } $$

$$⇔ |g|^2 = \frac{1}{ 1-4\nu(1-cos\theta )+4\nu ^2(1-cos\theta)^2 } $$

\(|g| \le 1\)より

$$⇔ \frac{1}{ 1-4\nu(1-cos\theta )+4\nu ^2(1-cos\theta)^2 } \le 1 $$

$$⇔ 0 \le -4\nu(1-cos\theta )+4\nu ^2(1-cos\theta)^2 $$

$$⇔ \frac{1}{1-cos\theta} \le \nu $$

$$⇔ 0\le \frac{1}{1-cos\theta} \le \nu $$

常に\(0 \le \nu\) であるので、上記の条件は常に満たされる。

よって陰解法を用いると拡散方程式においても、常に安定である。