import csv
import numpy as np
from math import exp
from matplotlib import pyplot as plt
"""
数値微分における異なる近似法の誤差を評価するスクリプトです。
このスクリプトは、指数関数 `exp(x)` の導関数を、様々なステップサイズ `h` と異なる数値微分近似法（2点前方、2点後方、3点中心、5点中心、7点中心差分）を用いて計算します。
計算された数値微分結果とその解析解 `exp(x)` との誤差を比較し、その結果をコンソール、CSVファイル、およびMatplotlibグラフに出力します。
関連リンク: :doc:`diff_order_usage`
"""
#===================
# parameters
#===================
outfile = 'diff_order.csv'
h = 0.5
kh = 0.8
nh = 101
x0 = 1.0
[ドキュメント]
def func(x):
"""
微分対象となる関数を定義します。
この関数は、指数関数 `exp(x)` を返します。
:param x: float - 関数の評価点。
:returns: float - `exp(x)` の計算結果。
"""
return exp(x)
[ドキュメント]
def dfdx(x):
"""
`func(x)` の解析的な導関数を定義します。
この関数は、`exp(x)` の導関数である `exp(x)` を返します。
:param x: float - 導関数の評価点。
:returns: float - `exp(x)` の導関数の計算結果。
"""
return exp(x)
[ドキュメント]
def diff2forward(func, x, h):
"""
2点前方差分法を用いて関数の数値を微分します。
この方法は、`f'(x) = (f(x+h) - f(x)) / h` という公式に基づいています。
:param func: callable - 微分対象の関数。
:param x: float - 微分を評価する点。
:param h: float - 微小なステップサイズ。
:returns: float - `x` における関数の数値微分値。
"""
return (func(x+h) - func(x)) / h;
[ドキュメント]
def diff2backward(func, x, h):
"""
2点後方差分法を用いて関数の数値を微分します。
この方法は、`f'(x) = (f(x) - f(x-h)) / h` という公式に基づいています。
:param func: callable - 微分対象の関数。
:param x: float - 微分を評価する点。
:param h: float - 微小なステップサイズ。
:returns: float - `x` における関数の数値微分値。
"""
return (func(x) - func(x-h)) / h;
[ドキュメント]
def diff3(func, x, h):
"""
3点中心差分法を用いて関数の数値を微分します。
この方法は、`f'(x) = (f(x+h) - f(x-h)) / (2h)` という公式に基づいています。
2点差分法よりも高精度な近似を提供します。
:param func: callable - 微分対象の関数。
:param x: float - 微分を評価する点。
:param h: float - 微小なステップサイズ。
:returns: float - `x` における関数の数値微分値。
"""
return (func(x+h) - func(x-h)) / 2.0 / h;
[ドキュメント]
def diff5(func, x, h):
"""
5点中心差分法を用いて関数の数値を微分します。
この方法は、より広範囲の点を使用することで、3点中心差分法よりもさらに高精度な近似を提供します。
:param func: callable - 微分対象の関数。
:param x: float - 微分を評価する点。
:param h: float - 微小なステップサイズ。
:returns: float - `x` における関数の数値微分値。
"""
return (-func(x+2.0*h) + 8.0*func(x+h) - 8.0*func(x-h) + func(x-2.0*h)) / 12.0 / h;
[ドキュメント]
def diff7(func, x, h):
"""
7点中心差分法を用いて関数の数値を微分します。
この方法は、さらに多くの点を使用することで、非常に高精度な近似を提供します。
:param func: callable - 微分対象の関数。
:param x: float - 微分を評価する点。
:param h: float - 微小なステップサイズ。
:returns: float - `x` における関数の数値微分値。
"""
return (func(x+3.0*h) - 9.0*func(x+2.0*h) + 45.0*func(x+h)
- 45.0*func(x-h) + 9.0*func(x-2.0*h) - func(x-3.0*h)) / 60.0 / h;
#===================
# main routine
#===================
[ドキュメント]
def main(x0, h, nh):
"""
異なる数値微分近似の誤差を計算し、結果をCSVに書き出し、グラフをプロットします。
指定された関数 `func` の `x0` における導関数を、様々なステップサイズ `h` と異なる差分近似法で計算し、
解析解との誤差を比較します。結果はコンソール、CSVファイル、およびMatplotlibグラフに出力されます。
:param x0: float - 微分を評価する中心点。
:param h: float - 初期ステップサイズ。この値から `kh` で減少させながら計算を繰り返します。
:param nh: int - ステップサイズを変化させる回数。
:returns: None
"""
print("Numerical differentiation using differnet approximations")
print("Write to [{}]".format(outfile))
# open outfile to write a csv file
f = open(outfile, 'w')
fout = csv.writer(f, lineterminator='\n')
fout.writerow([
'Ndiv', 'h',
'Df/Dx(2-point)',
'Df/Dx(3-point)',
'Df/Dx(5-point)',
'Df/Dx(7-point)'
])
fx = func(x0)
f1x = dfdx(x0)
print("")
print("Analytical values:")
print(" f({})={}".format(x0, fx))
print(" df/dx({})={}".format(x0, f1x))
print("{:^3}:\t{:^10}\t"
"{:^16}\t{:^16}\t{:^16}\t{:^16}"
.format(
'Ndiv', 'h',
'Df/Dx(2-point)',
'Df/Dx(3-point)',
'Df/Dx(5-point)',
'Df/Dx(7-point)'
)
)
xh = [h * kh**ih for ih in range(nh)]
ydiff2 = []
ydiff3 = []
ydiff5 = []
ydiff7 = []
errdiff2 = []
errdiff3 = []
errdiff5 = []
errdiff7 = []
for ih in range(nh):
print("{:^3}:\t{:^10.4e}\t".format(ih+1, xh[ih]), end = '')
ydiff2.append(diff2forward(func, x0, xh[ih]))
errdiff2.append(abs(ydiff2[ih] - f1x))
ydiff3.append(diff3(func, x0, xh[ih]))
errdiff3.append(abs(ydiff3[ih] - f1x))
ydiff5.append(diff5(func, x0, xh[ih]))
errdiff5.append(abs(ydiff5[ih] - f1x))
ydiff7.append(diff7(func, x0, xh[ih]))
errdiff7.append(abs(ydiff7[ih] - f1x))
print("{:^16.12f}\t{:^16.12f}\t{:^16.12f}\t{:^16.12f}\t".format(
ydiff2[ih], ydiff3[ih], ydiff5[ih], ydiff7[ih]), end = '')
fout.writerow([ih+1, xh[ih], ydiff2[ih], ydiff3[ih], ydiff5[ih], ydiff7[ih]])
print("")
f.close()
#=============================
# Plot graphs
#=============================
fig = plt.figure()
ax1 = fig.add_subplot(2, 1, 1)
ax2 = fig.add_subplot(2, 1, 2)
ax1.plot(xh, ydiff2, label = '$\Delta$f/$\Delta$x (2-points)', linewidth = 0.5)
ax1.plot(xh, ydiff3, label = '$\Delta$f/$\Delta$x (3-points)', linewidth = 0.5)
ax1.plot(xh, ydiff5, label = '$\Delta$f/$\Delta$x (5-points)', linewidth = 0.5)
ax1.plot(xh, ydiff7, label = '$\Delta$f/$\Delta$x (7-points)', linewidth = 0.5)
ax1.plot(ax1.get_xlim(), [f1x, f1x], label = 'df/dx (exact)', linestyle = 'dashed', linewidth = 0.5)
ax1.set_xlabel("h")
ax1.set_ylabel("$\Delta$f/$\Delta$x")
ax1.legend(bbox_to_anchor=(1.05, 1.0), loc='upper left', borderaxespad = 0)
ax2.plot(xh, errdiff2, label = 'error (2-points)', linewidth = 0.5)
ax2.plot(xh, errdiff3, label = 'error (3-points)', linewidth = 0.5)
ax2.plot(xh, errdiff5, label = 'error (5-points)', linewidth = 0.5)
ax2.plot(xh, errdiff7, label = 'error (7-points)', linewidth = 0.5)
ax2.set_xscale('log')
ax2.set_yscale('log')
ax2.set_xlabel("h")
ax2.set_ylabel("|$\Delta$f/$\Delta$x - df/dx|")
ax2.legend(bbox_to_anchor=(1.05, 1.0), loc='upper left', borderaxespad = 0)
plt.tight_layout()
plt.pause(0.1)
print("Press ENTER to exit>>", end = '')
input()
if __name__ == '__main__':
main(x0, h, nh)