無意味なブログ

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

MENU

【Python3で数値計算】4次精度陰的差分コードの解説

前回の記事で紹介した、4次精度陰的差分コードの解説をしていきます。

インポート

1,2行目

import numpy as np
import matplotlib.pyplot as plt

numpyは行列演算、matplotlib.pyplotはグラフの描画に使います。

定数の宣言

4-6行目

N = 10                       # 分割数
diff_range = 2.0 * np.pi     # 微分範囲
h = diff_range / N           # 刻み幅

差分の計算条件を設定しています。

参考サイト

stackoverflow.com
numpy.piと他のライブラリに含まれているpiの関係について書かれています。

解析解の設定

9行目

exact = np.array([np.cos(h * i) for i in range(0, N)])

あとで数値解と比較するため、解析解を設定しています。今回はsin(x)を微分するので、解析解はcos(x)です。

参考サイト

NumPy で数学系の関数を使ってみよう – Python でデータサイエンス
np.cos(x)の使用例が書かれています。

右辺の行列の初期化

12-17行目

b = np.zeros(N)
b[0] = np.sin(h) - np.sin(h * (N - 1))
b[N-1] = np.sin(0) - np.sin(h * (N - 2))
for i in range(1, N-1):
    b[i] = np.sin(h * (i + 1)) - np.sin(h * (i - 1))
b = b * 3.0 / h

今回解くべき方程式の右辺
 \displaystyle
\frac{1}{2h}
\begin{pmatrix}
f_{2}-f_{n} \\
\vdots \\
f_{i+1}-f_{i-1} \\
\vdots \\
f_{1}-f_{n-1} \\
\end{pmatrix}
を配列bに代入しています。

参考サイト

NumPy配列ndarrayの生成 | hydroculのメモ
要素をすべて0で埋めた配列を生成するnp.zerosの使用例が書かれています。

係数行列の初期化

20-22行目

A = np.diag([4] * N, k = 0) + np.diag([1] * (N - 1), k = 1) \
    + np.diag([1] * (N - 1), k = -1) \
    + np.diag([1], k = N-1) + np.diag([1], k = -N+1)

今回の係数行列
 \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}
を配列Aに代入しています。

参考サイト

numpy.diag — NumPy v1.13 Manual
numpy.diagのリファレンスです。numpy.diagは対角成分を作成するのに使えます。

d.hatena.ne.jp
numpy.diagの使用例が日本語で書かれています。

d.hatena.ne.jp
複数行の文を書くときの書き方について書かれています。

線形方程式を解く

25行目

x = np.linalg.solve(A, b)

線形方程式
 \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つ目の行列
 \displaystyle
\begin{pmatrix}
\left( \frac{df}{dx} \right)_1 \\
\vdots \\
\left( \frac{df}{dx} \right)_i \\
\vdots \\
\left( \frac{df}{dx} \right)_n \\
\end{pmatrix}
を求めて、配列xに代入しています。

参考サイト

org-technology.com
numpyを使った線形連立方程式の解き方が解説されています。linalg.solveの使用例が書かれています。

グラフの描画

28-32行目

plot_x = np.array([h * i for i in range(0, N)])
plt.plot(plot_x, x, label = "compact")
plt.plot(plot_x, exact, linestyle = "", marker = ".", label = "cos(x)")
plt.legend()
plt.show()

matplotlib.pyplotを使って、求めた数値解と解析解をグラフに描画しています。

参考サイト

lines — Matplotlib 2.0.2 documentation
グラフ描写のオプションについてはここを参考にしました。

ハマったところ

途中で

TypeError: 'float' object is not callable

というエラーが出るようになりました。

"float"って書いてあるので型関係で何かミスしたかなと思いましたが、原因は微分範囲を表す変数に"range"という名前をつけていたことでした。pythonにはrange()という組み込み関数があるのでそこと衝突していたようです。

参考サイト

stackoverflow.com