[matlab]经典数值微分计算方法的痛与痒

0 写在前面的话

数值微分和数值积分一对兄弟给人的第一印象是一件事的两面，但事实上前者是远比后者困难的一件事，数值微分作为一个不适定的问题，在人们认为十分自然的经典做法中却存在着很多坏的特性，接下来我们将看到数值微分无法通过加密网格尺度来提高收敛精度以及稍有扰动就会对数值微分结果产生非常恶劣的影响这两个坏消息是如何产生的。

1 问题

求解函数$f(x)$在$x_0$处的导数值：
$$f^{\prime }(x_0)=\underset{h\to 0}{\lim }\frac{f(x_0+h)-f(x_0)}{h}$$

2 格式简介

2.1 向前差分

设$f\in C^2$，则对问题有向前差分格式
$$f^{\prime }(x_0)\approx \frac{f(x_0+h)-f(x_0)}{h}$$
由Taylor展开可知
$$\frac{f(x_0+h)-f(x_0)}{h}-f^{\prime }(x_0)=\frac{h}{2}{f}^{\prime \prime }(x_0)+O(h^2)$$
理论上向前差分格式是一阶收敛的。

2.2 向后差分

设$f\in C^2$，则对问题有向后差分格式
$$f^{\prime }(x_0)\approx \frac{f(x_0)-f(x_0-h)}{h}$$
由Taylor展开可知
$$\frac{f(x_0)-f(x_0-h)}{h}-f^{\prime }(x_0)=-\frac{h}{2}{f}^{\prime \prime }(x_0)+O(h^2)$$
理论上向后差分格式是一阶收敛的。

2.3 中心差分

设$f\in C^3$，则对问题有中心差分格式
$$f^{\prime }(x_0)\approx \frac{f(x_0+h)-f(x_0-h)}{2h}$$
由Taylor展开可知
$$\frac{f(x_0+h)-f(x_0-h)}{2h}-f^{\prime }(x_0)=O(h^2)$$
理论上中心差分格式是二阶收敛的。

2.4 Richardson外推与中心差分结合

利用低阶格式，经过简单的变换可以得到高阶求解数值微分的格式：
当$f$的光滑性足够高时，由Taylor展开可知
$$f(x_0\pm h)=f(x_0)\pm h{f}^{\prime }(x_0)+\frac{h^2}{2}{f}^{\prime \prime }(x_0)\pm \frac{h^3}{6}{f}^{\prime \prime \prime }(x_0)+\frac{h^4}{24}{f}^{(4)}(x_0)\pm ...$$
于是对中心差分格式的误差可以精确表示为
$$e_{{}_h}^{(0)}=\frac{f(x_0+h)-f(x_0-h)}{2h}-f^{\prime }(x_0)=c_1{h}^2+c_2{h}^4+c_3{h}^6+...$$
$$e_{{}_{{}^h/{}_2}}^{(0)}=\frac{f(x_0+{}^h/{}_2)-f(x_0-{}^h/{}_2)}{h}-f^{\prime }(x_0)=c_1{(\frac{h}{2})}^2+c_2{(\frac{h}{2})}^4+c_3{(\frac{h}{2})}^6+...$$
将上面两式结合可以得到
$$e_h^{(1)}=e_h^{(0)}-4e_{{}^h/{}_2}^{(0)}=\widetilde{c_2}{h}^4+\widetilde{c_3}{h}^6+...=O(h^4)$$
由此可见每经过一层Richardson外推理论上可将差分格式的精度提高两阶。
这样可以得到改进的数值微分格式：
$$f^{\prime }(x_0)=D_k(h)+O(h^{2k+2})=D_{k-1}(h)+\frac{D_{k-1}(h)-D_{k-1}(2h)}{4^k-1}+O(h^{2k+2})$$
$$D_{\text{0}}(h)=\frac{f(x_0+h)-f(x_0-h)}{2h}$$

当光滑性足够好时，Richardson外推迭代理论上可以达到任意精度。

2.5 隐式方法

给定$f^{\prime }(x_0)$、$f^{\prime }(x_n)$，利用三次样条插值的思想可得到$f^{\prime }(x)$的近似。
若令$m_i=f^{\prime }(x_i)$，由Taylor展开式变换后有
$$m_i=\frac{f(x_{i+1})-f(x_{i-1})}{2h}-\frac{h^2}{6}{f}^{\prime \prime \prime }(x_i)+O(h^4)\text{\hspace{1em}(1)}$$
由Taylor展开可知，对足够光滑的函数$g$有
$$g(x\pm h)=g(x)\pm h{g}^{\prime }(x)+\frac{h^2}{2}{g}^{\prime \prime }(x)\pm \frac{h^3}{6}{g}^{\prime \prime \prime }(x)+O(h^4)\text{\hspace{1em}(2)}$$
由（1）容易得到$g$的二阶导数中心差分为
$$g^{\prime \prime }(x)=\frac{g(x+h)-2g(x)+g(x-h)}{h^2}+O(h^2)$$
令$g=f^{\prime }$，在$x_i$处有$f^{\prime \prime \prime }$的近似：
$$f^{\prime \prime \prime }(x_i)=\frac{m_{i+1}-2m_i+m_{i-1}}{h^2}+O(h^2)\text{\hspace{1em}(3)}$$
由（1）和（3）得到$f^{\prime }(x_i)$的隐式表达式：
$$m_i=\frac{f(x_{i+1})-f(x_{i-1})}{2h}-\frac{h^2}{6}\frac{m_{i+1}-2m_i+m_{i-1}}{h^2}+O(h^4)\text{\hspace{1em}(4)}$$
当取$m_i\approx \widetilde{m_i}=\frac{f(x_{i+1})-f(x_{i-1})}{2h}-\frac{h^2}{6}\frac{m_{i+1}-2m_i+m_{i-1}}{h^2}$ 时，由（4）显见得理论上隐式方法的误差为四阶。
经整理后，（4）可转化为如下求解线性方程组的问题
$$\frac{1}{6}\left(\begin{array}{ccccc}4& 1& & & \\ 1& 4& \ddots & & \\ & \ddots & \ddots & \ddots & \\ & & \ddots & 4& 1\\ & & & 1& 4\end{array}\right)\left(\begin{array}{c}{\tilde{m}}_1\\ {\tilde{m}}_2\\ ⋮\\ {\tilde{m}}_{n-2}\\ {\tilde{m}}_{n-1}\end{array}\right)=\frac{\text{1}}{2h}\left(\begin{array}{c}f(x_2)-f(x_0)\\ f(x_3)-f(x_1)\\ ⋮\\ f(x_{n-1})-f(x_{n-3})\\ f(x_n)-f(x_{n-2})\end{array}\right)-\frac{1}{6}\left(\begin{array}{c}{m}_0\\ 0\\ ⋮\\ 0\\ m_n\end{array}\right)$$
此方程组的系数矩阵总为三对角矩阵，利用追赶法可以在$O(n)$时间内求得方程组唯一解。

3 数值实验

3.1 算例1

$$f(x)=e^x;x=1$$

上图是五种方法在没有测量误差的前提下得出的。从左图可知五种方法的收敛性都是满足的，随着网格尺度$h$的减小，误差逐渐减小，呈现出了明显的正相关。右图是对左图的一个简单的变换，从右图可知向前差分、向后差分都是一阶收敛的，中心差分是二阶收敛的，隐式方法是四阶收敛的，进行两次外推的Richardson方法约为六至八阶收敛的。

上图是在没有测量误差的前提下，使用Richardson方法进行若干次外推和中心差分格式得出的。从图中可知Richardson外推每进行一次迭代收敛阶大约上升两阶，同时可以观察到相同的初始网格尺度$h$下Richardson方法外推层数越多越快达到机器精度，收敛曲线呈锯齿状，表明此时舍入误差的影响已经不能再忽略。

3.2 算例2

$$f(x)=\sin x;x=1$$

上图是五种方法在没有测量误差的前提下得出的，与算例1中相比较得到的结论是相同的。在将网格进一步加密以后，得到如下图：

从上图可知，随着网格的加密，受舍入误差的影响逐渐明显，最终所有方法都在某一点后不再继续收敛。

3.3 算例3

$$f(x)=\sin x;x=1$$
下图是五种方法在存在测量误差的前提下得出的，其中所加的测量误差是由白噪声模拟的。

从上两图可见，在存在噪声时，随着网格的加密，所有方法都会表现出不适定性，而且随着噪声强度$p$的加大，对网格尺度$h$的要求会越来越大。

4 总结

五种计算方法在以上三个算例中的表现与理论结果基本相符，其中向前差分、向后差分是一阶收敛的，中心差分是二阶收敛的，Richardson外推与中心差分结合后是$2k+2$阶收敛的，其中$k$是迭代次数，而隐式方法是四阶收敛的。

在对数值微分不适定性的观察中，舍入误差的影响随着网格尺寸$h$的减小逐渐明显，出现明显的锯齿状曲线先后顺序为：Richardson外推、隐式方法、中心差分、向前/后差分；噪声的影响随着噪声强度的增大而增大，往往在一个加上一个非常小的噪声后，结果就会变得非常糟糕；在相同噪声强度下，影响随着网格尺寸$h$的减小逐渐明显，出现明显的锯齿状曲线先后顺序为：Richardson外推、隐式方法、中心差分、向前/后差分。

从个人的理解中，笔者认为舍入误差和噪声对结果的影响在计算的初始阶段有着相似的特性，都是给精确值加上了一个小的误差$\delta$，由对向前差分的分析：
当$\left|f^{\delta }(x)-f(x)\right|\le \delta$时，
$$\left|\frac{f^{\delta }(x_0+h)-f^{\delta }(x_0)}{h}-f^{\prime }(x_0)\right|$$
$$=\left|\frac{f^{\delta }(x_0+h)-f^{\delta }(x_0)}{h}-\frac{f(x_0+h)-f(x_0)}{h}+\frac{f(x_0+h)-f(x_0)}{h}-f^{\prime }(x_0)\right|$$
$$\le ch+\frac{2\delta }{h}$$
以及同理对中心差分的分析：
$$\left|\frac{f^{\delta }(x_0+h)-f^{\delta }(x_0-h)}{2h}-f^{\prime }(x_0)\right|$$
$$=\left|\frac{f^{\delta }(x_0+h)-f^{\delta }(x_0-h)}{2h}-\frac{f(x_0+h)-f(x_0-h)}{2h}+\frac{f(x_0+h)-f(x_0-h)}{2h}-f^{\prime }(x_0)\right|$$
$$\le c{h}^2+\frac{\delta }{h}$$

故可知在存在一个小的误差$\delta$时，向前、向后、中心差分格式的网格尺度都不能无限小，在达到某个值后继续减小误差反而会增大，由于Richardson外推在第一层使用了中心差分，且在迭代中要不断将网格继续剖分，同时还有外推累积的误差，故Richardson外推法也会在网格尺度减小到某个值后误差会增大，不再继续收敛，而且这种现象出现的比前面几种低阶格式更早。
对于隐式方法，虽然对于整体不便于分析，但类比前面几种计算方法，误差估计式
$$\left|\widetilde{f^{\prime }(x_0)}-f^{\prime }(x_0)\right|\le O(h^n+\frac{\delta }{h});n>1$$
中$O(\frac{\delta }{h})$总来自于$f$在各点的误差除以网格尺度$h$，而$O(h^n)$则主要是截断Taylor展式高阶项带来的，于是由隐式方法的表达式
$$m_i=\frac{f(x_{i+1})-f(x_{i-1})}{2h}-\frac{h^2}{6}\frac{m_{i+1}-2m_i+m_{i-1}}{h^2}+O(h^4)$$
可知隐式方法应该也具有不适定性，即在有一个小误差存在时，网格尺度不能无限小，在达到某个值后继续减小误差反而会增大，以上的分析与实验结论也基本相符。

5 参考代码

接下来是皆大欢喜的MATLAB代码环节 ~

function main
%%
clear;clc;close all;
%%
%初始化
h = 10.^linspace(-0.25,-2.75);
L = length(h);
x0 = 1;%求导点
f=@func;%测试函数
E_Forward = zeros(1,L);
E_Backward = zeros(1,L);
E_Central = zeros(1,L);
E_Richardson = zeros(1,L);
E_Implicit = zeros(1,L);
%%
%计算误差
syms x,y=f(x);
dfunc=diff(y);
for i = 1:L
    Exact_solution=subs(dfunc,x,x0);
    E_Forward(i) = abs(forward(f,x0,h(i)) - Exact_solution);
    E_Backward(i) = abs(backward(f,x0,h(i)) - Exact_solution);
    E_Central(i) = abs(central(f,x0,h(i)) - Exact_solution);
    E_Richardson(i) = abs(richardson(f,x0,h(i),3) - Exact_solution);
    E_Implicit(i) = abs(implicit(f,x0,h(i),20) - Exact_solution);
end
%%
%绘图：log scale（这很重要）
figure;
subplot(1,2,1);
plot(log10(h),log10(E_Forward),'LineWidth',2);hold on;
plot(log10(h),log10(E_Backward),'--g','LineWidth',2);hold on;
plot(log10(h),log10(E_Central),'LineWidth',2);hold on;
plot(log10(h),log10(E_Richardson),'LineWidth',2);hold on;
plot(log10(h),log10(E_Implicit),'LineWidth',2);
legend("E Forward","E Backward","E Central","E Richardson","E Implicit");
title("收敛性分析图")
xlabel("网格：log10(h)")
ylabel("误差：log10(Error)")
grid minor
%%
end

function f = func(x)
%测试函数
%f = exp(x);%x=1
f = sin(x);%x=1
%f = x^(x^x);%x=0.0001
end
function df = forward(f,x0,h)
%向前差分
df = (f(x0+h)-f(x0)) / h;
end
function df = backward(f,x0,h)
%向后差分
df = (f(x0)-f(x0-h)) / h;
end
function df = central(f,x0,h)
%中心差分
df = (f(x0+h)-f(x0-h)) / (2*h);
end

function df = richardson(f,x0,h,N)
%中心差分+Richarson外推，N为外推层数,第一层为中心差分
df_tmp = zeros(1,N);
for i = 1:N
    df_tmp(i) = central(f,x0,2*h*0.5^i);
end
%disp("Richarson外推过程：")
%disp(df_tmp)
for i = 1:N
    j = N-i;
    df_tmp(1:j) = (df_tmp(2:j+1)*4^i - df_tmp(1:j)) / (4^i-1);
    %disp(df_tmp(1:j))
end
df = df_tmp(1);
end

function df=implicit(f,x0,h,N)
%隐格式，在x0附近计算2N+1个点的微分，默认最远两端点导数为已知
n=2*N+1;
F=zeros(n,1);
x1=x0-N*h;
for i=1:n
    F(i)=f(x1+i*h)-f(x1+(i-2)*h);
end
F=F/(2*h);
%求已知点导数
syms x,y=f(x);
dfun=diff(y);
F(1)=F(1)-subs(dfun,x,x1)/6;
F(n)=F(n)-subs(dfun,x,x0+N*h)/6;
%追赶法求解三对角系数矩阵的线性方程组
M=crout(2*ones(1,n)/3,ones(1,n-1)/6,ones(1,n-1)/6,F);
df=M(N+1);
end

function x=crout(a,c,d,b)
%追赶法 a：主对角线 c、d：次对角线
n=length(a);
n1=length(c);
n2=length(d);
p=1:n;
q=1:n-1;
x=1:n;
y=1:n;
p(1)=a(1);
for i=1:n-1
    q(i)=c(i)/p(i);
    p(i+1)=a(i+1)-d(i)*q(i);
end
y(1)=b(1)/p(1);
for i=2:n
    y(i)=(b(i)-d(i-1)*y(i-1))/p(i);
end
x(n)=y(n);
for i=(n-1):-1:1
    x(i)=y(i)-q(i)*x(i+1);
end
end