cms.diff_richardson のソースコード

import numpy as np
from math import exp

"""
リチャードソン補外法による数値微分を実装します。

このスクリプトは、指数関数exp(x)の微分をリチャードソン補外法を用いて数値的に計算し、
その結果を解析解と比較します。計算過程でステップサイズhを徐々に小さくし、
D行列を構築して精度の高い近似値を求めます。
指定された収束条件を満たすまで計算を継続します。

:doc:`diff_richardson_usage`
"""

# define function to be differentiated


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

    この関数は、xの指数関数 `exp(x)` を返します。

    :param x: 微分対象の関数の入力値です。
    :returns: 入力xに対する `exp(x)` の値です。
    """
    return exp(x)


# define analytical deviation of f(x)


def dfdx(x):
    """微分対象の関数f(x)の解析的な導関数f'(x)を定義します。

    この関数は、`exp(x)` の導関数である `exp(x)` を返します。

    :param x: 微分対象の関数の入力値です。
    :returns: 入力xに対する `exp(x)` の解析的な導関数の値です。
    """
    return exp(x)


#===================
# parameters
#===================
h0    = 1.0
nhmax = 15
eps   = 1.0e-6
x0 = 1.0

#===================
# main routine
#===================
print("Numerical differentiation using by Richardson extrapolation method")

fx = func(x0)
f1x = dfdx(x0)
print("")
print("eps for convergnece: {}".format(eps))
print("Analytical values:")
print("  f({})={}".format(x0, fx))
print("  df/dx({})={}".format(x0, f1x))
print("")

# create 2D array by specifing the size to be nhmax * nhmax
D = []
for i in range(nhmax):
  D.append([])
  for j in range(nhmax):
    D[i].append([])

D[0][0] = (func(x0+h0) - func(x0-h0)) / 2.0 / h0
print("D{}{}: {}".format(0, 0, D[0][0]))

for m in range(1, nhmax+1):
    h0 *= 0.5
    D[0][m] = (func(x0+h0) - func(x0-h0)) / 2.0 / h0
    print("D0{}: {}".format(m, D[0][m]))
    for k in range(m-1, -1, -1):
        m1 = m - k
        D[m1][k] = (4.0**m1 * D[m1-1][k+1] - D[m1-1][k]) / (4.0**m1 - 1)
        print("D{}{}: {}".format(m1, k, D[m1][k]))
    diff = D[m][0] - D[m-1][0]
    print("  diff = D[{}][0] - D[{}][0] = {} - {} = {}".format(m, m-1, D[m][0], D[m-1][0], diff))
    if abs(diff) < eps:
        break