数值微分

数值微分

张瑞
中国科学技术大学数学科学学院

数值微分

问题. 已知函数$f(x)$在离散点处的函数值

\[f(x_i),i=0,1,\cdots,n \]

如何求$f(x)$的导数？

插值函数的微分

  • 由插值理论，利用函数在离散处的函数值$f(x_i),i=0,1,\cdots,n$，可以得到插值多项式$L_n(x)$。
  • 利用插值多项式的性质，近似可以得到$f(x)$的特性。
  • 取插值多项式的微分为数值微分公式

    \[f'(x)\approx L'_n(x) \]

    并且，误差为插值多项式的误差的微分，

    \[\begin{aligned} R(x)&=f'(x)-L'_n(x)=\left(\frac{f^{(n+1)}(\xi)}{(n+1)!}\omega_n(x) \right)' \\ &=\frac{f^{(n+1)}(\xi)}{(n+1)!}\omega'_n(x)+\frac{\omega_n(x)}{(n+1)!}\frac{d}{dx}f^{(n+1)}(\xi) \end{aligned} \]

    其中$\omega_n(x)=(x-x_0)(x-x_1)\cdots(x-x_n)$

注意到$\omega_n(x_i)=0$，当数值微分公式在$x_i$处取值时，有

\[\begin{aligned} R_n(x_i)&=f'(x_i)-L'_n(x_i)=\frac{f^{(n+1)}(\xi)}{(n+1)!}\omega'_n(x_i) \\ &=\frac{f^{(n+1)}(\xi)}{(n+1)!}(x_i-x_0)\cdots(x_i-x_{i-1})(x_i-x_{i+1})\cdots(x_i-x_n) \end{aligned} \]

注. 做误差分析的时候，同样可以利用表达式在$x_i$处的Taylor展开来求

两点公式

做$f(x)$的线性插值公式

\[L_1(x)=\frac{x-x_0}{x_1-x_0}f(x_1)+\frac{x-x_1}{x_0-x_1}f(x_0) \]

对上式求导

\[L'_1(x)=\frac{f(x_1)}{x_1-x_0}+\frac{f(x_0)}{x_0-x_1} =\frac{f(x_0)-f(x_1)}{x_0-x_1} \]

记$h=x_1-x_0>0$，则有公式

\[f'(x_0)=\frac{f(x_1)-f(x_0)}{h}-\frac{h}2f''(\xi), \xi \in(x_0,x_1) \]

定义 1.
公式

\[f'(x)\approx \frac{f(x+h)-f(x)}h \]

是点$x$与$x+h$的差商，称为向前差商公式。

在$x_1$处取值，可以得到

\[f'(x_1)=\frac{f(x_1)-f(x_0)}{h}+\frac{h}2f''(\xi), \xi \in(x_0,x_1) \]

定义 2.
公式

\[f'(x)\approx \frac{f(x)-f(x-h)}h \]

是点$x$与$x-h$的差商，称为向后差商公式。

三点公式

在等距节点$x_i=x_0+ih,i=0,1,2$上，插值多项式为

\[\begin{aligned} L_2(x)=&\frac{(x-x_1)(x-x_2)}{(x_0-x_1)(x_0-x_2)}f(x_0) +\frac{(x-x_0)(x-x_2)}{(x_1-x_0)(x_1-x_2)}f(x_1) \\ &+\frac{(x-x_0)(x-x_1)}{(x_2-x_0)(x_2-x_1)}f(x_2) \\ \end{aligned} \]

令$x=x_0+th$，则有

\[\begin{aligned} L_2(x_0+th)=&\frac12(t-1)(t-2)f(x_0)-t(t-2)f(x_1)+\frac12t(t-1)f(x_2) \\ \frac{d}{dx}L_2(x_0+th)=&\frac1{2h}[(2t-3)f(x_0)-4(t-1)f(x_1)+(2t-1)f(x_2)] \\ \end{aligned} \]

这样，可以得到如下的三点公式

\[\begin{aligned} f'(x_0)&=\frac1{2h}(-3f(x_0)+4f(x_1)-f(x_2))+\frac{h^2}3f'''(\xi) \\ f'(x_1)&=\frac1{2h}(-f(x_0)+f(x_2))-\frac{h^2}6f'''(\xi)\ \ \ \ {\color{green}(**)} \\ f'(x_2)&=\frac1{2h}(f(x_0)-4f(x_1)+3f(x_2))+\frac{h^2}3f'''(\xi) \\ \end{aligned} \]

其中，(**)式，又称为中心差商公式。可以写成

\[f'(x)\approx \frac{f(x+h)-f(x-h)}{2h} \]

五点公式

记$f_i=f(x_i),i=0,1,2,3,4$，则有如下公式

\[\begin{split} f'(x_0)&=\frac1{12h}(-25f_0+48f_1-36f_2+16f_3-3f_4)+\frac{h^4}5f^{(5)}(\xi) \\ f'(x_1)&=\frac1{12h}(-3f_0-10f_1+18f_2-6f_3+f_4)-\frac{h^4}{20}f^{(5)}(\xi) \\ f'(x_2)&=\frac1{12h}(f_0-8f_1+8f_3-f_4)-\frac{h^4}{30}f^{(5)}(\xi) \\ f'(x_3)&=\frac1{12h}(-f_0+6f_1-18f_2+10f_3+3f_4)-\frac{h^4}{20}f^{(5)}(\xi) \\ f'(x_4)&=\frac1{12h}(3f_0-16f_1+36f_2-48f_3+25f_4)+\frac{h^4}5f^{(5)}(\xi) \\ \end{split} \]

例 1. 函数$f(x)=e^x$在$x=2.7$附近有如下数据。求$f'(2.7)$

x 2.5 2.6 2.7 2.8 2.9
y 12.1825 13.4637 14.8797 16.4446 18.1741

解. 运行结果， 

二点公式=============
Error for 2.6,2.7 -0.7197948261613831
Error for 2.7,2.8  0.7694187373692749
Error for 2.5,2.7 -1.3935429040260185
Error for 2.7,2.9  1.5923364979782804
三点公式============
Error for 2.5,2.6,2.7 -0.04604674829675659
Error for 2.6,2.7,2.8  0.02481195560394589
Error for 2.7,2.8,2.9 -0.053499023239712784
Error for 2.5,2.7,2.9  0.09939679697613002
五点公式==============
Error for 2.5,2.6,2.7,2.8,2.9 -4.965818678748235e-05

从误差的数据，可以得出哪些结论？

  • $h$越小，误差越小
  • 点越多，误差越小（五点公式比三点公式误差小，三点公式比二点公式误差小）
  • 在三点公式中，中心差商格式的误差是另外两个公式误差的近似一半。
  • 在三点公式中，$h=0.2$时的误差大概是$h=0.1$时误差的$4$倍。

例 2. 使用递减的步长$h$来计算函数$f(x)=e^x$在$x=1.15$处的微分。

步长$h=0.1$，每次减半，计算结果如下 

有效数字Decimal:6 
0  :   diff:  3.16345  err:  -0.00526
1  :   diff:  3.1595   err:  -0.00131
2  :   diff:  3.1584   err:  -0.00021
3  :   diff:  3.1584   err:  -0.00021
4  :   diff:  3.1576   err:   0.00059
5  :   diff:  3.1536   err:   0.00459
6  :   diff:  3.152    err:   0.00619
7  :   diff:  3.1552   err:   0.00299
8  :   diff:  3.1488   err:   0.00939
9  :   diff:  3.22561  err:  -0.06742
10  :   diff:  3.22561 err:  -0.06742
11  :   diff:  3.17441 err:  -0.01622
12  :   diff:  2.66241 err:   0.49578
13  :   diff:  2.45761 err:   0.70058
14  :   diff:  4.91521 err:  -1.75702
15  :   diff:  0E+5  err:  3.15819
16  :   diff:  0E+5  err:  3.15819
17  :   diff:  0E+6  err:  3.15819
18  :   diff:  0E+6  err:  3.15819
有效数字Decimal:16 
0: 1.15, 3.158192909689768 
0: 3.163459196993385, -0.005266287303617 
1: 3.15950898790114,  -0.001316078211372 
2: 3.15852189839860,  -0.000328988708832 
3: 3.15827515493932,  -0.000082245249552 
4: 3.15821347088168,  -0.000020561191912 
5: 3.15819804998016,  -0.000005140290392 
6: 3.15819419476192,  -0.000001285072152 
7: 3.15819323095808,  -3.21268312E-7 
8: 3.15819299000704,  -8.0317272E-8 
9: 3.15819292976896,  -2.0079192E-8 
10: 3.15819291470848, -5.018712E-9 
11: 3.15819291095040, -1.260632E-9 
12: 3.15819290998784, -2.98072E-10 
13: 3.15819290976256, -7.2792E-11 
14: 3.15819290976256, -7.2792E-11 
15: 3.15819290918912,  5.00648E-10 
16: 3.15819290918912,  5.00648E-10 
17: 3.15819290918912,  5.00648E-10 
18: 3.15819291312128, -3.431512E-9
...
46: 2.462906046218241,  0.695286863471527 
47: 4.925812092436481, -1.767619182746713 
48: 0E-15, 3.158192909689768 
49: 0E-15, 3.158192909689768

前面的例子可以看到，

  • 随着$h$越小，误差越小。
  • 但当$h$越过某个值后，再减小它，会增大误差。这是因为随着$h$越来越小，舍入误差变得越来越大。

问题. 如何确定一个合适的步长$h$，使得误差在要求的范围？

需要知道当前$h$时，误差有多大。

问题. 如何判定当前步长$h$时，误差有多大？

事后误差估计

设步长为$h$的近似$F(h)$与真值$Q$的误差为

\[Q-F(h)=C(h) h^p \]

则对步长$\frac{h}2$，有

\[Q-F(\frac{h}2)=C(\frac{h}2) (\frac{h}2)^p \]

近似地认为$C(h)\approx C(\frac{h}2)$，有

\[2^p(Q-F(\frac{h}2))\approx Q-F(h) \]

即有，

\[\begin{aligned} (2^p-1)Q-2^pF(\frac{h}2)\approx & -F(h) \\ (2^p-1)(Q-F(\frac{h}2))\approx & F(\frac{h}2)-F(h) \end{aligned} \]

这样，得到了事后误差估计式

\[Q-F(\frac{h}2)\approx \frac{F(\frac{h}2)-F(h)}{2^p-1} \]

即，可以用$\frac{F(\frac{h}2)-F(h)}{2^p-1}$来近似作为$F(\frac{h}2)$的误差。

例 3. 中心差商格式的误差是$O(h^2)$，这样，可以用

\[\frac{F(\frac{h}2)-F(h)}{2^2-1}=\frac13(F(\frac{h}2)-F(h)) \]

来近似计算误差的值。

例 4. 向前差商的误差是$O(h)$，可以用$F(\frac{h}2)-F(h)$近似来估计误差。

给定误差要求$\epsilon$，如何确定一个步长$h$，使得数值微分的误差小于$\epsilon$？ 

h=0.1
df1=center_euler(h)  #中心差商公式计算微分
df2=center_euler(h/2)
while |df1-df2|>3*eps:  #事后误差估计公式
    df1=df2
    h=h/2                 # 步长变为原来的一半
    df2=center_euler(h/2)

return df2  # df2的误差满足要求

Richardson外推

假定步长$h$时的数值计算公式$F(h)$的误差表达式为

\[Q-F(h)=c\cdot h^p+O(h^r) \ \ \ \ \ {\color{green}(*)} \]

其中$r>p$。 将步长减半，则有

\[Q-F(\frac{h}2)=c (\frac{h}2)^p+O(h^r) \]

得到

\[2^p\left(Q-F(\frac{h}2)\right)+O(h^r)=c\cdot h^p \]

代入(*)式，有

\[\begin{aligned} Q-F(h)&=2^p(Q-F(\frac{h}2))+O(h^r) \\ \end{aligned} \]

由

\[\begin{aligned} Q-F(h)&=2^p(Q-F(\frac{h}2))+O(h^r) \\ \end{aligned} \]

得到

\[(1-2^p)Q -F(h)+2^pF(\frac{h}2)=O(h^r) \]

进而

\[Q - \frac{F(h)-2^pF(\frac{h}2)}{1-2^p}=O(h^r) \]

可以得到

\[Q-\frac{2^pF(\frac{h}2)-F(h)}{2^p-1}=O(h^r) \]

即公式

\[\frac{2^pF(\frac{h}2)-F(h)}{2^p-1} \]

与$Q$的误差为$O(h^r)$，是比$F(h)$（误差为$O(h^p)$）更高阶的近似。 这个公式称为Richardson外推公式。

可以看到

\[\frac{2^pF(\frac{h}2)-F(h)}{2^p-1}=F(\frac{h}2)+\frac{F(\frac{h}2)-F(h)}{2^p-1} \]

即Richardson外推公式就是近似值加上事后误差估计值。

例 5. 对中心差商格式，

\[\frac{1}{2h}(f(x+h)-f(x-h)) \]

在$x$处做Taylor展开，得到

\[f(x+h)=f(x)+f'(x)h+f''(x)\frac{h^2}{2!}+f'''(x)\frac{h^3}{3!}+f^{(4)}(x)\frac{h^4}{4!}+f^{(5)}(\xi_1)\frac{h^5}{5!} \]
\[f(x-h)=f(x)-f'(x)h+f''(x)\frac{h^2}{2!}-f'''(x)\frac{h^3}{3!}+f^{(4)}(x)\frac{h^4}{4!}-f^{(5)}(\xi_2)\frac{h^5}{5!} \]

因而， 有误差公式

\[f'(x)=\frac{1}{2h}(f(x+h)-f(x-h))-\frac{h^2}6f'''(x)-\frac{h^4}{5!}f^{(5)}(\xi) \]

用Richardson外推公式，可以得到

\[\begin{aligned} F_4(x)&=\frac{2^2F(\frac{h}2)-F(h)}{2^2-1} \\ &=\frac13\left(4\frac{f(x+\frac{h}2)-f(x-\frac{h}2)}{h}-\frac{f(x+h)-f(x-h)}{2h}\right) \\ &=\frac{1}{6h}\left(f(x-h)-8f(x-\frac{h}2)+8f(x+\frac{h}2)-f(x+h)\right) \end{aligned} \]

为五点公式。具有$O(h^4)$的精度。

高阶导数

利用插值多项式，同样可以近似高阶导数。

如，三点公式中，求2次插值多项式的2阶导数，有

\[L_2''(x)=\frac{f(x_0)-2f(x_1)+f(x_2)}{h^2} \]

可以得到2阶导数的数值近似

\[\begin{aligned} f''(x_1) & \approx \frac{f(x_0)-2f(x_1)+f(x_2)}{h^2} \\ &\approx \frac{f(x_1-h)-2f(x_1)+f(x_1+h)}{h^2} \end{aligned} \]

问题. 求函数的三阶导数，至少需要几个点？

误差分析

将公式在$x_1$处做Taylor展开

\[\begin{aligned} f(x_1+h)=f(x_1)+hf'(x_1)+\frac{h^2}2f''(x_1)+\frac{h^3}6f'''(x_1)+\frac{h^4}{4!}f^{(4)}(\xi_1) \\ f(x_1-h)=f(x_1)-hf'(x_1)+\frac{h^2}2f''(x_1)-\frac{h^3}6f'''(x_1)+\frac{h^4}{4!}f^{(4)}(\xi_2) \\ \end{aligned} \]

则有，

\[\begin{aligned} f''(x_1)&-\frac{f(x_1-h)-2f(x_1)+f(x_1+h)}{h^2} \\ &=-(\frac{h^2}{4!}f^{(4)}(\xi_1)+\frac{h^2}{4!}f^{(4)}(\xi_2)) =-\frac{h^2}{12}f^{(4)}(\xi) \\ &=-\frac{h^2}{12}f^{(4)}(x_1)+\frac{h^4}{2\times 6!}f^{(6)}(\xi) \end{aligned} \]

谢谢

vertical slide 2

目录

本节读完

例 6.

6.