import csv
import numpy as np
from math import exp, sin, cos
from random import random, uniform
from matplotlib import pyplot as plt
"""
ノイズを含む数値データに対する異なる近似次数による数値微分の効果を評価するスクリプトです。
このスクリプトは、ノイズを含むサンプルデータに対して、様々な次数の数値微分スキーム（2点前方差分、
2点後方差分、3点中心差分、5点中心差分、7点中心差分）を適用し、その結果を真の導関数と比較します。
計算されたデータはコンソールに出力されるとともに、CSVファイルに保存されます。
また、元の関数とノイズ付きデータ、および各微分スキームによる導関数の推定値をグラフで可視化し、
ノイズが数値微分に与える影響と、近似次数の違いによる精度の変化を示します。
:doc:`diff_order_noise_usage`
"""
# define function to be differentiated
[ドキュメント]
def func(x):
"""
微分対象の関数f(x)を定義します。
この関数は、与えられた独立変数xに対してsin(x)の値を返します。
:param x: float 独立変数x。
:returns: float sin(x)の値。
"""
return sin(x)
# define analytical deviation of f(x)
[ドキュメント]
def dfdx(x):
"""
関数f(x)の解析的な導関数f'(x)を定義します。
この関数は、`func(x) = sin(x)` の解析的な導関数である `cos(x)` の値を返します。
:param x: float 独立変数x。
:returns: float cos(x)の値（f(x)=sin(x)の導関数）。
"""
return cos(x)
[ドキュメント]
def make_data(x0, xstep, nx, noise):
"""
指定された範囲とステップで、ノイズを含むデータセットを生成します。
関数 `func(x)` の値に一様乱数によるノイズを加算します。
ノイズは `noise * (np.random.rand() - 0.5) * 2.0` の形式で加えられ、
`[-noise, +noise]` の範囲で変動します。
データは中心の `nx` 点の前後3点ずつ（計 `nx+6` 点）生成されます。
:param x0: float データの開始x座標。
:param xstep: float 各点のx間隔。
:param nx: int データ点の数（中心部分の点数）。
:param noise: float データに加えるノイズの振幅。
:returns: tuple of list (xa, ya) - x座標のリストと、ノイズを含むy座標のリスト。
"""
xa = []
ya = []
for i in range(-3, nx+3):
x = x0 + i * xstep
xa.append(x)
ya.append(func(x) + noise * (np.random.rand() - 0.5) * 2.0)
# ya.append(func(x) + noise * (random() - 0.5) * 2.0)
# ya.append(func(x) + noise * uniform())
return xa, ya
# numerical differentiation by two-point forward difference method
[ドキュメント]
def diff2forward(x, y, i):
"""
2点前方差分法を用いて数値微分を計算します。
対象インデックス `i` と `i+1` のデータ点 `(x[i], y[i])` と `(x[i+1], y[i+1])` を使用します。
配列の範囲外 (`i` が `n-1` の場合) では計算できず、空文字列を返します。
:param x: list[float] x座標のリスト。
:param y: list[float] y座標のリスト。
:param i: int 微分を計算する点のインデックス。
:returns: float または str 指定された点における数値微分の値、または計算不可の場合は空文字列。
"""
n = len(x)
if(0 <= i <= n-2):
return (y[i+1] - y[i]) / (x[i+1] - x[i]);
return ''
[ドキュメント]
def diff2backward(x, y, i):
"""
2点後方差分法を用いて数値微分を計算します。
対象インデックス `i` と `i-1` のデータ点 `(x[i], y[i])` と `(x[i-1], y[i-1])` を使用します。
配列の範囲外 (`i` が `0` の場合) では計算できず、空文字列を返します。
:param x: list[float] x座標のリスト。
:param y: list[float] y座標のリスト。
:param i: int 微分を計算する点のインデックス。
:returns: float または str 指定された点における数値微分の値、または計算不可の場合は空文字列。
"""
n = len(x)
if(1 <= i <= n-1):
return (y[i] - y[i-1]) / (x[i] - x[i-1]);
return ''
[ドキュメント]
def diff3(x, y, i):
"""
3点中心差分法を用いて数値微分を計算します。
対象インデックス `i-1` と `i+1` のデータ点 `(x[i-1], y[i-1])` と `(x[i+1], y[i+1])` を使用します。
配列の範囲外 (`i` が `0` または `n-1` の場合) では計算できず、空文字列を返します。
:param x: list[float] x座標のリスト。
:param y: list[float] y座標のリスト。
:param i: int 微分を計算する点のインデックス。
:returns: float または str 指定された点における数値微分の値、または計算不可の場合は空文字列。
"""
n = len(x)
if(1 <= i <= n-2):
return (y[i+1] - y[i-1]) / 2.0 / (x[i] - x[i-1]);
return ''
[ドキュメント]
def diff5(x, y, i):
"""
5点中心差分法を用いて数値微分を計算します。
対象インデックス `i-2`, `i-1`, `i+1`, `i+2` のデータ点を使用します。
配列の範囲外 (`i` が `0, 1` または `n-1, n-2` の場合) では計算できず、空文字列を返します。
:param x: list[float] x座標のリスト。
:param y: list[float] y座標のリスト。
:param i: int 微分を計算する点のインデックス。
:returns: float または str 指定された点における数値微分の値、または計算不可の場合は空文字列。
"""
n = len(x)
if(2 <= i <= n-3):
return (-y[i+2] + 8.0*y[i+1] - 8.0*y[i-1] + y[i-2]) / 12.0 / (x[i] - x[i-1]);
return ''
[ドキュメント]
def diff7(x, y, i):
"""
7点中心差分法を用いて数値微分を計算します。
対象インデックス `i-3`, `i-2`, `i-1`, `i+1`, `i+2`, `i+3` のデータ点を使用します。
配列の範囲外 (`i` が `0, 1, 2` または `n-1, n-2, n-3` の場合) では計算できず、空文字列を返します。
:param x: list[float] x座標のリスト。
:param y: list[float] y座標のリスト。
:param i: int 微分を計算する点のインデックス。
:returns: float または str 指定された点における数値微分の値、または計算不可の場合は空文字列。
"""
n = len(x)
if(3 <= i <= n-4):
return (y[i+3] - 9.0*y[i+2] + 45.0*y[i+1]
- 45.0*y[i-1] + 9.0*y[i-2] - y[i-3]) / 60.0 / (x[i] - x[i-1]);
return ''
#===================
# parameters
#===================
outfile = 'diff_order_noise.csv'
x0 = 0.0
xstep = 0.2
nx = 51
noise = 0.2
#===================
# main routine
#===================
[ドキュメント]
def main():
"""
プログラムのメインエントリポイントです。
ノイズを含むデータセットを生成し、異なる次数の数値微分を適用して結果を比較します。
計算されたデータはCSVファイルに出力され、コンソールにも表示されます。
最後に、元の関数と導関数、および様々な数値微分の結果をプロットします。
:returns: None
"""
print("Numerical differentiation using differnet approximations")
print("Write to [{}]".format(outfile))
xa, ya = make_data(x0, xstep, nx, noise)
# open outfile to write a csv file
f = open(outfile, 'w')
fout = csv.writer(f, lineterminator='\n')
fout.writerow([
'x', 'f(x)', 'f+noise', 'df/dx',
'Df/Dx(2-point)',
'Df/Dx(3-point)',
'Df/Dx(5-point)',
'Df/Dx(7-point)'
])
print("")
print("{:^6}:\t{:^10}\t{:^10}\t{:^10}\t"
"{:^10}\t{:^10}\t{:^10}\t{:^10}"
.format(
'x', 'f(x)', 'f+noise', 'df/dx',
'Df/Dx(2-point)',
'Df/Dx(3-point)',
'Df/Dx(5-point)',
'Df/Dx(7-point)'
))
xx = []
yf = []
yf1 = []
ydiff2 = []
ydiff3 = []
ydiff5 = []
ydiff7 = []
for ix in range(3, nx+3): # make_dataで前後3点ずつ多めに生成されているため、中心nx点のみを処理
idx = ix - 3
x = xa[ix]
y = ya[ix]
xx.append(x)
yf.append(func(x))
yf1.append(dfdx(x))
ydiff2.append(diff2forward(xa, ya, ix))
ydiff3.append(diff3(xa, ya, ix))
ydiff5.append(diff5(xa, ya, ix))
ydiff7.append(diff7(xa, ya, ix))
print("{:^5.3f}:\t{:^10.4f}\t{:^10.4f}\t{:^10.4f}\t".format(x, y, yf[idx], yf1[idx]), end = '')
print("{:^10.4f}\t{:^10.4f}\t{:^10.4f}\t{:^10.4f}\t".format(
ydiff2[idx], ydiff3[idx], ydiff5[idx], ydiff7[idx]), end = '')
print("")
fout.writerow([x, y, yf[idx], yf1[idx], ydiff2[idx], ydiff3[idx], ydiff5[idx], ydiff7[idx]])
f.close()
#=============================
# Plot graphs
#=============================
fig = plt.figure()
ax1 = fig.add_subplot(2, 1, 1)
ax2 = fig.add_subplot(2, 1, 2)
ax1.plot(xx, yf, label = 'f(x)(w/o noise)', linewidth = 0.5)
ax1.plot(xa, ya, label = 'f(x)(with noise)', linewidth = 0.5)
ax1.set_xlabel("x")
ax1.set_ylabel("f(x)")
ax1.legend(bbox_to_anchor=(1.05, 1.0), loc='upper left', borderaxespad = 0)
ax2.plot(xx, yf1, label = 'df/dx (w/o noise)')
ax2.plot(xx, ydiff2, label = r'$\Delta$f/$\Delta$x (2-points)', linewidth = 0.5)
ax2.plot(xx, ydiff3, label = r'$\Delta$f/$\Delta$x (3-points)', linewidth = 0.5)
ax2.plot(xx, ydiff5, label = r'$\Delta$f/$\Delta$x (5-points)', linewidth = 0.5)
ax2.plot(xx, ydiff7, label = r'$\Delta$f/$\Delta$x (7-points)', linewidth = 0.5)
ax2.set_xlabel("x")
ax2.set_ylabel("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()