数值积分的python实现——NewtonCotes、复化求积、Romberg、richardson递推

\qquad 使用牛顿-莱布尼兹公式直接去求一个函数 $f(x)$ 在区间 $[a,b]$ 上的定积分值 $I$，在数值方法上比较困难，因为需要先求出被积函数的原函数 $F(x)$。

$I={\int }_a^{b}f(x)dx=F(b)-F(a)$

\qquad 通过积分中值定理，在积分区间 $[a,b]$ 上存在某个点 $\xi$，使得

$I={\int }_a^{b}f(x)dx=(b-a)f(\xi ),\xi \in (a,b)$
\qquad

1. 机械求积

\qquad 根据积分中值定理，只需要合理选择 $\xi \in (a,b)$ 的值，就可以求出定积分，例如最基本的梯形公式和中矩形公式。
\qquad

用红色虚线矩形的面积，代替 $f(x)$ 与 $x$ 轴所围的面积。

• 梯形公式：　　取 $f(\xi )\approx \frac{f(a)+f(b)}{2}$

${\int }_a^{b}f(x)dx=\frac{b-a}{2}[f(a)+f(b)]$

• 中矩形公式：　取 $f(\xi )\approx f\left(\frac{a+b}{2}\right)$

${\int }_a^{b}f(x)dx=(b-a)f\left(\frac{a+b}{2}\right)$

\qquad 更一般地，可以在积分区间 $[a,b]$ 上选取一些求积节点 $x_k$，用 $f(x_k)$ 的加权平均（权值为 $A_k$）来近似 $f(\xi )$ 的值，也就是采用“机械求积”方法：

$f(\xi )\approx \sum _{k=0}^n{A}_{k}f(x_k)$$⟹{\int }_a^{b}f(x)dx=(b-a)\sum _{k=0}^n{A}_{k}f(x_k)$

梯形公式：　 $n=2,A_0=1,x_0=a,A_1=1,x_1=b$
中矩形公式： $n=1,A_0=1,x_0=\frac{a+b}{2}$

2. Newton-Cotes公式

\qquad 将积分区间 $[a,b]$ 划分为 $n$ 等分，步长为 $h=\frac{b-a}{n}$，选取等距求积节点 $x_k=a+kh$ 构造出牛顿-科特斯 $\text{(Newton-Cotes)}$ 公式：

$I_n=(b-a)\sum _{k=0}^n{C}_k^{(n)}f(x_k)$

\qquad 其中，$C_k^{(n)}$ 为科特斯系数，若 $x=a+th$，则

$C_k^{(n)}=\frac{h}{b-a}{\int }_0^n\prod _{\genfrac{}{}{0ex}{}{j=0}{j\ne k}}^n\frac{t-j}{k-j}dt=\frac{(-1{)}^{n-k}}{nk!(n-k)!}{\int }_0^n\prod _{\genfrac{}{}{0ex}{}{j=0}{j\ne k}}^n(t-j)dt$

\qquad

• 当 $n=1$ 时，$C_0^{(1)}=C_1^{(1)}=\frac{1}{2}$，也就是梯形公式（单区间，$x_0=a,x_1=b$）
  
  $T=I_1=\frac{b-a}{2}[f(a)+f(b)]$

• 当 $n=2$ 时，$C_0^{(2)}=\frac{1}{6},C_1^{(2)}=\frac{4}{6},C_2^{(2)}=\frac{1}{6}$，也就是辛普森$\text{(Simpson)}$公式（$2$ 等分区间）
  
  $S=I_2=\frac{b-a}{6}\left[f(a)+4f\left(\frac{a+b}{2}\right)+f(b)\right]$

• 当 $n=4$ 时，也就是科特斯公式（$4$ 等分区间）
  
  $C=I_4=\frac{b-a}{90}\left[7f(a)+32f(x_1)+12f(x_2)+32f(x_3)+7f(b)\right]$
  \qquad 其中， $x_0=a,x_4=b,x_k=a+kh,h=\frac{b-a}{4}$
  　
  【注意】当 $n\ge 8$ 时，科特斯系数 $C_k^{(n)}$ 出现负值，求积计算变得不稳定。

\qquad

3. 复化求积方法

\qquad 由于高阶的牛顿-科特斯公式不稳定，无法通过提高阶数 $n$ 来提高求积精度。

\qquad 为了提高精度，可以采用复化求积方法：将积分区间 $[a,b]$ 分成若干个子区间 $[x_k,x_{k+1}]$，在每个子区间上采用低阶 $(n<8)$ 的牛顿-科特斯公式。
\qquad

3.1 复化梯形公式

\qquad 将积分区间 $[a,b]$ 进行 $n$ 等分，步长为 $h=\frac{b-a}{n}$，等距求积节点 $x_k=a+kh,k=0,1,\cdots ,n$，在每个子区间 $[x_k,x_{k+1}],k=0,1,\cdots ,n-1$ 上采用梯形公式：

$T^{(k)}=\frac{x_{k+1}-x_k}{2}[f(x_k)+f(x_{k+1})]=\frac{h}{2}[f(x_k)+f(x_{k+1})]$

\qquad 记 $n$ 等分时求得积分的结果为 $T_n$，那么

${\int }_a^{b}f(x)dx\approx T_n=\sum _{k=0}^{n-1}{T}^{(k)}=\frac{h}{2}\sum _{k=0}^{n-1}[f(x_k)+f(x_{k+1})]$

误差为 $R_n(f)=-\frac{b-a}{12}{h}^2{f}^{{}^{\prime \prime }}(\eta ),\eta \in (a,b)$ 为 $h^2$ 阶

\qquad

3.2 复化辛普森公式

\qquad 将积分区间 $[a,b]$ 进行 $n$ 等分，步长为 $h=\frac{b-a}{n}$，等距求积节点 $x_k=a+kh,k=0,1,\cdots ,n$，在每个子区间 $[x_k,x_{k+1}],k=0,1,\cdots ,n-1$ 上采用辛普森公式：

$\begin{array}{rl}{S}^{(k)}& =\frac{x_{k+1}-x_k}{6}\left[f(x_k)+4f\left(x_{k+\frac{1}{2}}\right)+f(x_{k+1})\right]\\ & =\frac{h}{6}\left[f(x_k)+4f\left(x_{k+\frac{1}{2}}\right)+f(x_{k+1})\right]\end{array}$

\qquad 记 $n$ 等分时求得积分的结果为 $S_n$，那么

${\int }_a^{b}f(x)dx\approx S_n=\sum _{k=0}^{n-1}{S}^{(k)}=\frac{h}{6}\sum _{k=0}^{n-1}\left[f(x_k)+4f\left(x_{k+\frac{1}{2}}\right)+f(x_{k+1})\right]$

误差 $R_n(f)=-\frac{b-a}{180}{\left(\frac{h}{2}\right)}^4{f}^{(4)}(\eta ),\eta \in (a,b)$ 为 $h^4$ 阶

实现代码

例：求积分 $I={\int }_0^1\frac{\sin (x)}{x}dx$，$7$ 位有效数字准确值为 $0.9460831$

import numpy as np

def trapezoid(f,a,b):    # trapezoid(n=1)
    return (b-a)*(f(b)+f(a))/2
    
def simpson(f,a,b):    # simpson(n=2)
    t = (a+b)/2
    val = (b-a)*(f(a)+4*f(t)+f(b))/6
    return val

def compositeInt(f,a,b,n,mode='trapezoid'):
    h = (b-a)/n
    val = 0
    eps = 1e-10
    if mode=='trapezoid':
        for k in range(n):
            xk = a+k*h+eps
            xk1 = xk+h
            val += trapezoid(f,xk,xk1)        
    if mode=='simpson':
        for k in range(n):
            xk = a+k*h+eps
            xk1 = xk+h
            val += simpson(f,xk,xk1)            
    return val

if __name__ == '__main__':
    a = 0
    b = 1
    f = lambda x: np.sin(x)/x    
    print('====== trapezoid: newton-cotes(1st-order) ======')
    val = compositeInt(f, a, b, 8)
    print('trapezoid: %.7f' % (val))     
    print('====== simpson: newton-cotes(2nd-order) ======')
    val = compositeInt(f, a, b, 4,'simpson')
    print('simpson: %.7f' % (val))

运行结果：
====== trapezoid: newton-cotes(1st-order) ======
　　　trapezoid: 0.9456909　　（2位有效数字）
====== simpson: newton-cotes(2nd-order) ======
　　　simpson: 0.9460833　　 （6位有效数字）
　　　
$(1)$ 由于 $4$ 等分的复化辛普森方法在每个子区间上还需要再计算一个二分点，因此区间的等分数实际上是 $8$ 等分。
$(2)$ 然而，$4$ 等分的复化辛普森方法的精度高于 $8$ 等分的复化梯形方法。由误差公式也可计算出 $|R_8(f)|\le 0.000434$ 以及 $|R_4(f)|\le 0.271\times 10^{-6}$。
\qquad

4. 求积公式的递推化

\qquad 复化求积方法如果要提高精度，就必然增加求积子区间的数量。例如，将 $[a,b]$ 进行 $n$ 等分时，求积节点的数量为 $n+1$；若将每个求积区间 $[x_k,x_{k+1}]$ 再一分为二，求积节点的数量变为 $2n+1$，每次等分之后，求积区间的数量就要翻倍。如果仍然在每个子区间上使用复化公式，无疑会极大地增加计算量。

\qquad 为了提高计算效率，可以将等分前后的积分值与进行比较，从而形成“递推关系”。
\qquad

4.1 复化梯形公式的递推化

\qquad 为了提高复化求积公式的精度，可以选择将每个求积子区间 $[x_k,x_{k+1}]$ 一分为二的方式，也就在区间 $[x_k,x_{k+1}]$ 增加一个等距的中间点 $x_{k+\frac{1}{2}}=\frac{1}{2}(x_k+x_{k+1})$。

\qquad

• 假设将 $[a,b]$ 等分 $n$ 个点 $a,x_1,\cdots ,x_{n-2},b$，步长 $h=\frac{b-a}{n}$，考虑复化梯形公式：

$1)$ 以求积子区间 $[x_k,x_{k+1}]$ 为例，$n$ 等分时的梯形公式为

$T_1=\frac{h}{2}[f(x_k)+f(x_{k+1})]$

$2)$ 进行 $2n$ 等分时，需要将每个求积子区间 $[x_k,x_{k+1}]$ 一分为二，变成了两个等距的子区间 $[x_k,x_{k+\frac{1}{2}}]$ 和 $[x_{k+\frac{1}{2}},x_{k+1}]$，在两个子区间上分别采用梯形公式。此时，$2n$ 等分时的复化梯形公式为是：

$\begin{array}{rl}{T}_2& =\frac{h}{4}[f(x_k)+f(x_{k+\frac{1}{2}})]+\frac{h}{4}[f(x_{k+\frac{1}{2}})+f(x_{k+1})]\\ & =\frac{h}{4}[f(x_k)+2f(x_{k+\frac{1}{2}})+f(x_{k+1})]\\ & =\frac{h}{4}[f(x_k)+f(x_{k+1})]+\frac{h}{2}f(x_{k+\frac{1}{2}})\\ & =\frac{1}{2}{T}_2+\frac{h}{2}f(x_{k+\frac{1}{2}})\end{array}$

• 考虑整个求积区间 $[a,b]$，$n$ 等分时的复化梯形求积为：

$T_n=\frac{h}{2}\sum _{k=0}^{n-1}[f(x_k)+f(x_{k+1})]$

\qquad 按照上述规则，$2n$ 等分时的复化梯形求积为：

$T_{2n}=\frac{h}{4}\sum _{k=0}^{n-1}[f(x_k)+f(x_{k+1})]+\frac{h}{2}\sum _{k=0}^{n-1}f(x_{k+\frac{1}{2}})$

\qquad
\qquad
\qquad 对比这两个式子，可得到 $n$ 等分与 $2n$ 等分时采用复化梯形公式求积结果的关系：

$T_{2n}=\frac{1}{2}{T}_n+\frac{h}{2}\sum _{k=0}^{n-1}f(x_{k+\frac{1}{2}})$

\qquad 也就是说，从 $n$ 等分到 $2n$ 等分，不需要再在 $2n$ 个新的求积区间上使用复化梯形公式重新计算，只需要在 $n$ 等分求积结果 $T_n$ 的基础上，加上 $n$ 个新的求积节点的值 $f(x_{k+\frac{1}{2}})$ 即可。

实现代码

仍以求积分 $I={\int }_0^1\frac{\sin (x)}{x}dx$为例：

import numpy as np
def trapezoid(f,a,b):    			# trapezoid(n=1)
    return (b-a)*(f(b)+f(a))/2

def compositeInt_trape(f,a,b,n):    # recursive composite integration
    eps = 1e-10
    val = trapezoid(f,a+eps,b)
    print('T1: %.7f' % (val))
    i = 1
    while i<=np.log2(n):
        h = (b-a)/(2**i)
        t = 0
        for k in range(1,2**(i-1)+1):
            xkc = a + (2*k-1)*h
            t += f(xkc)*h    
        val = val/2 + t
        print('T%d: %.7f' % (2**i,val))
        i += 1
    return val    

if __name__ == '__main__':
    a = 0
    b = 1
    f = lambda x: np.sin(x)/x    
    n = 3           # 区间等分次数（3次等分后分为8个等距区间）
    val = compositeInt_trape(f, a, b, 2**n)
    print('recursive trapezoid: %.7f' % (val)) 

运行结果：
　　　T1: 0.9207355
　　　T2: 0.9397933
　　　T4: 0.9445135
　　　T8: 0.9456909
　　recursive trapezoid: 0.9456909
　　　
\qquad 本例中等分了 $3$ 次之后的子区间数为 $2^3=8$，和第 $3$ 节例子中复化梯形公式的结果相同，但是计算量减少了。　　　
\qquad

4.2 龙贝格(Romberg)算法

$\text{Romberg}$算法，是另一种减少计算量的算法，也是基于递推的方式。

• 分析梯形公式的误差
  \qquad
  $R_n(f)=-\frac{b-a}{12}{h}^2{f}^{{}^{\prime \prime }}(\eta ),\eta \in (a,b)$
  \qquad
  \qquad 可以得到：$\{\begin{array}{rl}I-T_n& =-\frac{b-a}{12}{h}^2{f}^{{}^{\prime \prime }}(\eta ),\eta \in (a,b)\\ & \\ I-T_{2n}& =-\frac{b-a}{12}{\left(\frac{h}{2}\right)}^2{f}^{{}^{\prime \prime }}(\bar{\eta }),\bar{\eta }\in (a,b)\end{array}$

\qquad
\qquad 若假设 $f^{{}^{\prime \prime }}(\eta )\approx f^{{}^{\prime \prime }}(\bar{\eta })$，那么
\qquad
$\frac{I-T_{2n}}{I-T_n}\approx \frac{1}{4}⟹I-T_{2n}\approx \frac{1}{3}(T_{2n}-T_n)$
\qquad
\qquad 从而可以通过递推关系得到积分值估计：$\overline{T}=T_{2n}+\frac{1}{3}(T_{2n}-T_n)=\frac{4}{3}{T}_{2n}-\frac{1}{3}{T}_n$

\qquad 又由　$T_n=\frac{h}{2}\sum _{k=0}^{n-1}[f(x_k)+f(x_{k+1})]$
\qquad 　和　$T_{2n}=\frac{h}{4}\sum _{k=0}^{n-1}[f(x_k)+f(x_{k+1})]+\frac{h}{2}\sum _{k=0}^{n-1}f(x_{k+\frac{1}{2}})$

\qquad 可得到积分值估计值 $\overline{T}$ 实际上就是辛普森公式的积分值 $S_n$：

$\begin{array}{rl}\overline{T}& =\frac{4}{3}{T}_{2n}-\frac{1}{3}{T}_n\\ & =\frac{h}{6}\sum _{k=0}^{n-1}[f(x_k+f(x_{k+1})]+\frac{2h}{3}\sum _{k=0}^{n-1}f(x_{k+\frac{1}{2}})\\ & =\frac{h}{6}\sum _{k=0}^{n-1}[f(x_k)+4f(x_{k+\frac{1}{2}})+f(x_{k+1})]=S_n\end{array}$
\qquad

• 依次类推，考虑（2等分）辛普森公式的误差

$\frac{I-S_{2n}}{I-S_n}\approx \frac{1}{16}⟹I-S_{2n}\approx \frac{1}{16}(S_{2n}-S_n)$
\qquad
\qquad 从而可以通过递推关系得到积分值估计：$\overline{S}=\frac{16}{15}{S}_{2n}-\frac{1}{15}{S}_n$

\qquad 验证可知积分值估计值 $\overline{S}$ 就是科特斯公式的积分值：$C_n=\frac{16}{15}{S}_{2n}-\frac{1}{15}{S}_n$
\qquad

• 依次类推，考虑（4等分）科特斯公式的误差

$\frac{I-C_{2n}}{I-C_n}\approx \frac{1}{64}⟹I-C_{2n}\approx \frac{1}{64}(C_{2n}-C_n)$
\qquad
\qquad 从而可以通过递推关系得到积分值估计：$\overline{C}=\frac{64}{63}{C}_{2n}-\frac{1}{63}{C}_n$

\qquad 也就得到了龙贝格$(\text{Romberg})$公式：

$R_n=\frac{64}{63}{C}_{2n}-\frac{1}{63}{C}_n$

\qquad 使用龙贝格公式时，需要首先将求积区间 $[a,b]$ 等分 $2^3=8$ 次。在这些子区间上，采用复化梯形公式进行求积运算。其具体过程可以表示成：

$k$	$h$	$T_{2^k}$	$S_{2^{k-1}}$	$C_{2^{k-2}}$	$R_{2^{k-3}}$
$0$	$b-a$	$T_1=0.92073549$
$1$	$\frac{b-a}{2}$	$T_2=0.93979328$	$S_1$
$2$	$\frac{b-a}{4}$	$T_4=0.94451352$	$S_2$	$C_1$
$3$	$\frac{b-a}{8}$	$T_8=0.94569086$	$S_4$	$C_2$	$R_1$

表中，$k$ 表示等分次数，$h$ 表示子区间长度，等分后的子区间数为 $2^k$
　　　$T_{2^k}$ 表示 $k$ 次等分时，采用复化梯形公式求积分的值（详见4.1节代码运行结果）
从第 $4$ 列开始，就是采用龙贝格递推算法的加速过程，不需要再进行复化求积。

4.3 理查森(Richardson)外推算法

\qquad 龙贝格算法实际上是理查森$\text{(Richardson)}$外推算法的一种特殊情况。

梯形公式的余项，可以展开成级数形式(${\alpha }_l$ 与 $h$ 无关)：
$T(h)=I+{\alpha }_1{h}^2+{\alpha }_2{h}^4+\cdots +{\alpha }_l{h}^{2l}+\cdots$
$T(\frac{h}{2})=I+{\alpha }_1\frac{h^2}{4}+{\alpha }_2\frac{h^4}{16}+\cdots +{\alpha }_l(\frac{h}{2}{)}^{2l}+\cdots$
$⟹T_1(h)=\frac{4T(\frac{h}{2})-T(h)}{3}=I+{\beta }_1{h}^4+{\beta }_2{h}^6+\cdots$ （辛普森公式序列）
$T_1(\frac{h}{2})=I+{\beta }_1\frac{h^4}{16}+{\beta }_2\frac{h^6}{64}+\cdots$
$⟹T_2(h)=\frac{16T_1(\frac{h}{2})-T_1(h)}{15}=I+{\gamma }_1{h}^6+{\gamma }_2{h}^8+\cdots$ （科特斯公式序列）
$\cdots \cdots$

\qquad 将龙贝格算法进行扩展，并进行公式化，可得到理查森$\text{(Richardson)}$递推关系：

$T_m^{(k)}=\frac{4^m}{4^m-1}{T}_{m-1}^{(k+1)}-\frac{1}{4^m-1}{T}_{m-1}^{(k)},k=1,2,\cdots$

\qquad 特别地，当选择 $m=3$ 时，即为龙贝格算法：

$(1)$　$m=1$ 时，$T_1^{(k)}=\frac{4}{3}{T}_0^{(k+1)}-\frac{1}{3}{T}_0^{(k)}$

$(2)$　$m=2$ 时，$T_2^{(k)}=\frac{16}{15}{T}_1^{(k+1)}-\frac{1}{15}{T}_1^{(k)}$

$(3)$　$m=3$ 时，$T_3^{(k)}=\frac{64}{63}{T}_2^{(k+1)}-\frac{1}{63}{T}_2^{(k)}$

\qquad 当 $m>3$ 时，即为理查森$\text{(Richardson)}$外推算法，如下表所示：

$k$	$h$	$T_0^{(k)}(T_{2^k})$	$T_1^{(k)}(S_{2^{k-1}})$	$T_2^{(k)}(C_{2^{k-2}})$	$T_3^{(k)}(R_{2^{k-3}})$	$T_4^{(k)}$
$0$	$b-a$	$T_0^{(1)}:T_1=0.92073549$
$1$	$\frac{b-a}{2}$	$T_0^{(2)}:T_2=0.93979328$	$T_1^{(1)}:S_1$
$2$	$\frac{b-a}{4}$	$T_0^{(3)}:T_4=0.94451352$	$T_1^{(2)}:S_2$	$T_2^{(1)}:C_1$
$3$	$\frac{b-a}{8}$	$T_0^{(4)}:T_8=0.94569086$	$T_1^{(3)}:S_4$	$T_2^{(2)}:C_2$	$T_3^{(1)}:R_1$
$4$	$\frac{b-a}{16}$	$T_0^{(5)}:T_{16}=0.94598503$	$T_1^{(4)}:S_8$	$T_2^{(3)}:C_4$	$T_3^{(2)}:R_2$	$T_4^{(4)}$
$⋮$	$⋮$	$⋮$	$⋮$	$⋮$	$⋮$	$⋮$

从第 $4$ 列开始，每个表项的值满足其左侧、左上侧的递推关系。

实现代码

仍以求积分 $I={\int }_0^1\frac{\sin (x)}{x}dx$为例：

import numpy as np

def trapezoid(f,a,b):    # trapezoid(n=1)
    return (b-a)*(f(b)+f(a))/2

def compositeInt_trape(f,a,b,n):    #recursive
    eps = 1e-10
    num = np.log2(n).astype(int)
    val = np.zeros(num+1)
    val[0] = trapezoid(f,a+eps,b)
    i = 1
    while i<=np.log2(n):
        h = (b-a)/(2**i)
        t = 0
        for k in range(1,2**(i-1)+1):
            xkc = a + (2*k-1)*h
            t += f(xkc)*h    
        val[i] = val[i-1]/2 + t
        i += 1
    return val

def romberg(f,a,b):
    # 只需要等分3次
    val_t = compositeInt_trape(f, a, b, 2**3)    
    t1 = np.concatenate((np.ones(1),val_t[0:-1]))
    t2 = (4*val_t - t1)/3
    val_s = t2[1:]    
    t1 = np.concatenate((np.ones(1),val_s[0:-1]))
    t2 = (16*val_s - t1)/15
    val_c = t2[1:]    
    t1 = np.concatenate((np.ones(1),val_c[0:-1]))
    t2 = (64*val_c - t1)/63
    val_r = t2[1]    
    print(val_t[-1::-1].flatten())
    print(val_s[-1::-1].flatten())
    print(val_c[-1::-1].flatten())
    print('[%.8f]' % (val_r))    
    return val_r

def richardson(f,a,b,k):
    n = k+2
    val_t = compositeInt_trape(f, a, b, 2**n)
    print(val_t[-1::-1].flatten())    
    for i in range(1,n):
        t0 = 4**i
        t1 = np.concatenate((np.ones(1),val_t[0:-1]))
        t2 = (t0*val_t - t1)/(t0-1)
        val_t = t2[1:]
        print(val_t[-1::-1].flatten())
    print('[%.8f]' % (val_t[1]))    
    return val_t[1]

if __name__ == '__main__':
    a = 0
    b = 1
    f = lambda x: np.sin(x)/x    
    val = romberg(f, a, b)
    print('romberg: %.7f' % (val))    
    val = richardson(f, a, b, 2)	# 比Romberg多等分一次
    print('richardson: %.7f' % (val))

运行结果：
[0.94569086 0.94451352 0.93979328 0.92073549]
[0.94608331 0.94608693 0.94614588]
[0.94608307 0.946083 ]
[0.94608307]
romberg: 0.9460831

只需要进行 $2^3=8$ 等分，并分别计算复合梯形公式的结果$T_1,T_2,T_4,T_8$ ($T_8$ 只具有2位有效数字)，就可以递推出具有更高精度的结果 $T_3^{(1)}$ (具有7位有效数字)。

[0.94598503 0.94569086 0.94451352 0.93979328 0.92073549]
[0.94608309 0.94608331 0.94608693 0.94614588]
[0.94608307 0.94608307 0.946083 ]
[0.94608307 0.94608307]
[0.94608307]
richardson: 0.9460831

参考：
数值微分的python实现