注解

此笔记本可在此处下载： 03_Numerical_Derivation.ipynb

数值数据的推导

代码作者：Emile Roux emile.roux@univ-smb.fr

RISE 幻灯片

范围

有限数 \(N\) 数据点数量 \((x,y)\) 是可用的，它描述了一个函数 \(f(x)\) ：计算导数 \(\dfrac{{df}}{{dx}}\)

https://en.wikipedia.org/wiki/Numerical_differentiation

#setup
%load_ext autoreload
%matplotlib nbagg
%autoreload 2
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib as mpl

创建合成数据

对于本笔记本，我们使用来自已知函数的数据混合。这样我们就可以检查结果的准确性。

第一个函数

# The function
def f(x):
    return (x**2-3*x+2*np.sin(5*x))

# The derivative
def df(x):
    return 2*x - 3 + 10 * np.cos(5*x)

第二功能

# The function
def f(x):
    return (x**2)

# The derivative
def df(x):
    return 2*x

绘制函数及其嘲笑：

N = 400
xmin, xmax = 0., 4
x = np.linspace(xmin, xmax, N)

fig = plt.figure()
plt.plot(x,f(x),'b',label = "$f(x)$")
plt.legend()
plt.grid()

<IPython.core.display.Javascript object>

fig = plt.figure()
plt.plot(x,df(x),'b',label = "$\dfrac{df}{dx}$")
plt.legend()
plt.grid()

<IPython.core.display.Javascript object>

现在我们假设这个函数只有40个点是已知的：

N = 40
xmin, xmax = 0., 4
xi = np.linspace(xmin, xmax, N)
fi = f(xi)

fig = plt.figure()
plt.plot(xi,fi,'o',label = "the discretize data $(x_i,f_i)$")
plt.legend()
plt.grid()

<IPython.core.display.Javascript object>

一点有限差分公式

\(\dfrac{df}{dx} = \dfrac{f(x_{i+1})-f(x_i)}{x_{i+1}-x_i}\)

numpy diff（）函数是计算此公式的快速方法：

df_1p = np.diff(fi)/np.diff(xi)

但是becarefull生成的数组的大小等于 \(n-1\) .

fig = plt.figure()
plt.plot(x,df(x),'b',label = "Analytical derivative")
plt.plot(xi[:-1],df_1p,'ro',label = "Finite difference derivative (1 point)")
plt.legend()
plt.grid()

<IPython.core.display.Javascript object>

两点有限差分公式

\(\dfrac{df}{dx} = \dfrac{f(x_{i+1})-f(x_{i-1})}{x_{i+1}-x_{i-1}}\)

numpy数组的切片方法是一种很好的操作方法

# Slicing exemple
a=np.linspace(0,6,7)
print(a[2:])
print(a[:-2])

[ 2.  3.  4.  5.  6.]
[ 0.  1.  2.  3.  4.]

使用numpy数组的切片方法，可以快速计算2点公式：

df_2p = (fi[2:] - fi[:-2])/(xi[2:] - xi[:-2])

fig = plt.figure()
plt.plot(x,df(x),'b',label = "Analytical derivative")
plt.plot(xi[:-1],df_1p,'ro',label = "Finite difference derivative (1 point)")
plt.plot(xi[1:-1],df_2p,'gs',label = "Finite difference derivative (2 points)")
plt.legend()
plt.grid()

<IPython.core.display.Javascript object>

我们来看看导数计算的误差

fig = plt.figure()
plt.plot(xi[:-1],(df_1p-df(xi[:-1])),'ro',label = "Finite difference derivative (1 point)")
plt.plot(xi[1:-1],(df_2p-df(xi[1:-1])),'gs',label = "Finite difference derivative (2 points)")
plt.legend()
plt.grid()

<IPython.core.display.Javascript object>

如果数据有噪声，附加什么？

fi_noise = f(xi) + 0.2*np.random.randn(fi.size)

fig = plt.figure()
plt.plot(x,f(x),':r',label = "$f(x)$")
plt.plot(xi,fi_noise,'o',label = "noisy data")
plt.legend()
plt.grid()

<IPython.core.display.Javascript object>

两点法更准确。让我们用它来推导噪声数据：

df_2p_n = (fi_noise[2:] - fi_noise[:-2])/(xi[2:] - xi[:-2])

fig = plt.figure()
plt.plot(x,df(x),':r',label = "Analytical derivative")
plt.plot(xi[1:-1],df_2p_n,'-g',label = "Finite difference derivative of noisy data")
plt.legend()
plt.grid()

<IPython.core.display.Javascript object>