無意味なブログ

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

MENU

【Python3で数値計算】4次精度陰的差分法の実装

Numpyの勉強がてら、昔の記事で紹介した、4次精度陰的差分法をpythonで実装したいと思います。

4次精度陰的差分法

昔の記事のコピペになりますが、今回使う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}が求めたい微分値となります。
↑コピペ終わり↑

今回もsin(x)を0から2πの範囲で微分したいと思います。

実装

さて、実装です。
完成したソースコードを以下に載せます。

行列Aを初期化する部分(20行目)はすごい読みづらい式になってますね。ただ、係数行列の初期化を一つの式でできちゃうっていうのはなかなか驚きです。
線形方程式を解くメソッドがある(25行目)というのも、便利ですね。

計算結果

matplotlibのpyplotというモジュールを使ってグラフ化します。

分割数N=10の時の結果です。

f:id:kouya17:20170813031119p:plain

青線が4次精度陰的差分法により得られた数値解、オレンジの点が解析解(cos(x))です。よく一致しています。


続いて分割数N=100の時の結果です。

f:id:kouya17:20170813031411p:plain

うん、なめらかですね。

まとめ

4次精度陰的差分法をpythonで実装しました。
fortranしか触ってこなかった自分にとって、Numpyはもの凄く便利に感じます。
もし余裕があったら、今回作成したプログラムの細かい説明なども記事にしたいと思います。