Numerial differentiation by Richardson extrapolation method
Download script from .\diff_richardson.py
import numpy as np
from math import exp
"""
Numerial differentiation by Richardson extrapolation method
"""
# define function to be differentiated
def func(x):
return exp(x)
# define analytical deviation of f(x)
def dfdx(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
