cms.diff_h のソースコード

import csv
import numpy as np

"""
  数値微分におけるステップ幅hの影響を評価するスクリプト。

  詳細説明:
    さまざまなhの値を用いて、2点前方差分法による数値微分の精度を検証します。
    計算結果はCSVファイルに出力され、コンソールにも表示されます。

  関連リンク:
    :doc:`diff_h_usage`
"""

# define function to be differentiated


def func(x):
    """
    微分対象の関数 `f(x) = x^3` を定義します。

    :param x: float: 関数評価点。
    :returns: float: `x` における関数の値 `x^3`。
    """
    return x**3


# define analytical deviation of f(x)


def dfdx(x):
    """
    `f(x) = x^3` の厳密な導関数 `f'(x) = 3x^2` を定義します。

    :param x: float: 導関数評価点。
    :returns: float: `x` における関数の厳密な導関数の値 `3x^2`。
    """
    return 3.0 * x * x


# numerical differentiation by two-point forward difference method


def diff2forward(func, x, h):
    """
    2点前方差分法を用いて関数の数値微分を計算します。

    詳細説明:
      微分対象の関数 `func` の `x` における導関数を、
      ステップ幅 `h` を用いた前方差分近似 `(func(x+h) - func(x)) / h` で計算します。

    :param func: callable: 微分する関数。
    :param x: float: 微分を計算する点。
    :param h: float: 微分計算に使用するステップ幅。
    :returns: float: `x` における関数の数値微分値。
    """
    return (func(x+h) - func(x)) / h;


#===================
# parameters
#===================
outfile = 'diff_h.csv'
h = [1.0, 0.1, 0.01, 0.001, 1.0e-6]
x0     = 0.0
xstep  = 0.1
nx     = 21

#===================
# main routine
#===================
print("Numerical differentiation using two-point forward difference method")
print("Write to [{}]".format(outfile))

# open outfile to write a csv file
f = open(outfile, 'w')
fout = csv.writer(f, lineterminator='\n')
fout.writerow([
    'x', 'f(x)', 'df/dx', 
    'Df/Dx(h={})'.format(h[0]),
    'Df/Dx(h={})'.format(h[1]),
    'Df/Dx(h={})'.format(h[2]),
    'Df/Dx(h={})'.format(h[3]),
    'Df/Dx(h={})'.format(h[4])
    ])

print("")
print("{:^5}:\t{:^12}\t{:^12}\t"
    "{:^12}\t{:^12}\t{:^12}\t{:^12}\t{:^12}"
    .format(
    'x', 'f(x)', 'df/dx', 
    'h={}'.format(h[0]), 
    'h={}'.format(h[1]), 
    'h={}'.format(h[2]), 
    'h={}'.format(h[3]), 
    'h={}'.format(h[4])) 
    )

for ix in range(nx):
    x = x0 + ix * xstep
    fx = func(x)
    f1x = dfdx(x)
    print("{:^5.3f}:\t{:^10.6f}\t{:^10.6f}\t".format(x, fx, f1x), end = '')
    diff = []
    for i in range(len(h)):
        diff.append(diff2forward(func, x, h[i]))
        print("{:^12.8f}\t".format(diff[i]), end = '')
    fout.writerow([x, fx, f1x] + diff)
    print("")

f.close()