python Numdifftool的有限差分导数（Fornberg方法）

Posted 2021-05-09

import numpy as np
import numdifftools.fornberg as ndf
import matplotlib.pyplot as plt
from random import random


def plot(outfile, x, fx=None, dfx=None, dfx_ref=None):
	plt.figure()
	if fx is not None:
		plt.plot(x, fx, label=r"$f(x)$", marker="*")
	if dfx is not None:
		plt.plot(x, dfx, label=r"$f'(x)$ (calculated)", marker="+")
	if dfx_ref is not None:
		plt.plot(x, dfx_ref, label=r"$f'(x)$ (reference)", marker="x")
	plt.legend()
	plt.xlabel("$x$")
	plt.savefig(outfile)
	plt.close()


def noisify(x, multiplier=1.):
	xmin, xmax = x.min(), x.max()
	deltax = xmax - xmin
	n = len(x)
	h = deltax / float(n-1)

	for i in range(len(x)):
		x[i] = x[i] * (1. + (2. * random() - 1.) / (float(n-1)) * multiplier)

	return x


# 1st order differentiation with even grid spacing
x = np.linspace(-1, 1, 21)
fx = np.tanh(x)
dfx = ndf.fd_derivative(fx, x, n=1)
dfx_ref = 1. - np.power(np.tanh(x), 2)
plot("even_tanh.png", x, fx, dfx, dfx_ref)

x = np.linspace(-np.pi, np.pi, 21)
fx = np.sin(x)
dfx = ndf.fd_derivative(fx, x, n=1)
dfx_ref = np.cos(x)
plot("even_sin.png", x, fx, dfx, dfx_ref)

x = np.linspace(0, 1, 21)
fx = np.exp(x)
dfx = ndf.fd_derivative(fx, x, n=1)
plot("even_exp.png", x, fx, dfx)

# 1st order differentiation with uneven (randomized) grid spacing
x = np.sort(noisify(np.linspace(-1, 1, 21), 5.))
fx = np.tanh(x)
dfx = ndf.fd_derivative(fx, x, n=1)
dfx_ref = 1. - np.power(np.tanh(x), 2)
plot("uneven_tanh.png", x, fx, dfx, dfx_ref)

x = np.sort(noisify(np.linspace(-np.pi, np.pi, 21), 5.))
fx = np.sin(x)
dfx = ndf.fd_derivative(fx, x, n=1)
dfx_ref = np.cos(x)
plot("uneven_sin.png", x, fx, dfx, dfx_ref)

x = np.sort(noisify(np.linspace(0, 1, 21), 5.))
fx = np.exp(x)
dfx = ndf.fd_derivative(fx, x, n=1)
plot("uneven_exp.png", x, fx, dfx)