無意味なブログ

勉強のこと、趣味のこと、日記など適当に

MENU

微分方程式をプログラムで解こう!(4):差分法[3]「陰的な差分法」

前回(第3回)
微分方程式をプログラムで解こう!(3):差分法[2]「差分法の精度」 - 無意味なブログ

前回は差分法の精度について説明し、分割数nを変化させた時、結果がどう変わるのかみてみました。
今回は"陰的"に解く差分法について説明します。

陽的と陰的

差分法には陽的な差分法と陰的な差分法があります。
これまでに紹介した前進差分、後退差分、中心差分は全て陽的な差分法です。
陽的な差分法と陰的な差分法にはどのような違いがあるのでしょうか。

最近使われている言葉として"陽キャ"、"陰キャ"という言葉がありますね。
陽的な差分法と陰的な差分法もその違いに近いものを持っていると考えてもらって良いです。

つまり、陽的な差分法は”陽キャ”、陰的な差分法は”陰キャ”というわけです。
陽的な差分法は積極的なため行動が早く、シンプルに物事を捉えられますが、結果の正しさに難がある"陽キャ"タイプ。
陰的な差分法は結果の正しさは優れていますが、慎重なため一つの行動にかける時間が長く、考えが複雑化してしまう"陰キャ"タイプと分類できます。

はい。クソ意味不明なことをいいました。僕も疲れているみたいです。

差分法を使う立場からすると、重要となる陽的な差分法と陰的な差分法の違いは、”連立一次方程式を解く必要があるかないか”です。
陽的な差分法は連立一次方程式を解く必要が"なく"、陰的な差分法は連立一次方程式を解く必要が”あり”ます。
陰的な差分法は連立一次方程式を解く必要があるため、陽的な差分法より解法は複雑化しますが、数値安定性はよくなります。

今回は陰的な差分法の一つ、4次精度陰的差分法*1について説明します。

4次精度陰的差分法

今回使う4次精度陰的差分法は以下の式で表されます。

 \displaystyle
 \frac{1}{6} \left( \left( \frac{df}{dx} \right)_{i-1} + 4 \left( \frac{df}{dx} \right)_i + \left( \frac{df}{dx} \right)_{i+1} \right) = \frac{f_{i+1}-f_{i-1}}{2h}

この式は以下の書籍を参考にしました。

上の式を行列形式で表すと以下のようになります(ただし、今回は周期的な関数を微分すると仮定しています)。

 \displaystyle
\frac{1}{6}
\begin{pmatrix}
 4 & 1 & & & 1 \\
 1 & \ddots & \ddots &  \huge{0} & \\
    & \ddots & \ddots & \ddots & \\
    &  \huge{0} & \ddots & \ddots & 1 \\
 1 &    & & 1 & 4 
\end{pmatrix}
\begin{pmatrix}
\left( \frac{df}{dx} \right)_1 \\
\vdots \\
\left( \frac{df}{dx} \right)_i \\
\vdots \\
\left( \frac{df}{dx} \right)_n \\
\end{pmatrix}
=
\frac{1}{2h}
\begin{pmatrix}
f_{2}-f_{n} \\
\vdots \\
f_{i+1}-f_{i-1} \\
\vdots \\
f_{1}-f_{n-1} \\
\end{pmatrix}

この式について左から2個目の行列\begin{pmatrix} \left( \frac{df}{dx} \right)_1 \cdots \left( \frac{df}{dx} \right)_i \cdots \left( \frac{df}{dx} \right)_n \end{pmatrix}が求めたい微分値となります。
この行列を求めるためには、連立一次方程式を解く必要があります。
この連立一次方程式について、今回は試しにLU分解を用いて解いてみます。

sin(x)を微分してみる

さあ、4次精度陰的差分を用いてsin(x)を微分してみましょう!
今回用いたプログラムは以下の通りです。

xの範囲は0から2π、分割数nは30としました。
このプログラムを実行すると、ファイル"diffIm.dat"に結果が出力されます。

結果

4次精度陰的差分による結果と解析解(cos(x))を比較します。
グラフの緑線が解析解(cos(x))、紫の点が差分によって得られた数値解です。

f:id:kouya17:20170212150500p:plain

よく一致していますね。
ちゃんと微分できています。

まとめ

4次精度陰的差分法を使ってsin(x)を微分してみました。
解析解と結果を比較してみたところ、よく一致しました。

シリーズ第4回「差分法[3]陰的な差分法」は以上です。
次回は陰的な差分法と陽的な差分法の結果を比較します。

*1:今回用いた"陰的差分法"は"コンパクト差分"という呼び方もされます。 調べるときは"コンパクト差分"で調べた方が色々出てくるかもしれません。