数值计算day7-数值微分

上一节课主要介绍了曲线拟合与插值，曲线拟合主要包括线性拟合（单特征线性回归和非线性拟合（非线性方程特征变换、高阶多项式拟合），插值包括多项式插值（拉格朗日形式、牛顿形式）、样条插值（线性插值、二次样条插值、三次样条插值），其中三次样条插值还有一个便于求解的拉格朗日形式。这里的曲线拟合与机器学习中的回归问题非常相似，具有很大的参考意义。本节课主要介绍几种求解微分的数值方法。

1. 有限差分法

给定一个函数(f(x))，(f(x))在(x=a)处的微分(f'(x))定义为：$$frac{df(x)}{dx}|{x=a}=f'(a)=lim{x ightarrow} frac{f(x)-f(a)}{x-a}$$ 图上的解释是(f(x))在(x=a)处的斜率：

前向、后向以及中心差分法是最简单的有限差分法：

• 前向差分：$$frac{df}{dx}|{x=x_i}=frac{f(x{i+1})-f(x_i)}{x_{i+1}-x_i}$$

• 后向差分：$$frac{df}{dx}|{x=x_i}=frac{f(x{i})-f(x_{i-1})}{x_{i}-x_{i-1}}$$

• 中心差分：$$frac{df}{dx}|{x=x_i}=frac{f(x{i+1})-f(x_{i-1})}{x_{i+1}-x_{i-1}}$$

例：(f(x)=x^3)，计算(f(x))在(x=3)处的微分
(a) (x=2,x=3,x=4)
真实值：(f'(x)=3x^2),(f'(3)=27)
前向差分：(f'(3)=frac{f(4)-f(3)}{4-3}=64-27=37)
后向差分：(f'(3)=frac{f(3)-f(2)}{3-2}=27-8=19)
中心差分：(f'(3)=frac{f(4)-f(2)}{4-2}=frac{64-8}{2}=28)
(b) (x=2.75,x=3,x=3.25)
前向差分：(f'(3)=frac{f(3.25)-f(3)}{0.25}=29.3125)
后向差分：(f'(3)=frac{f(3)-f(2.75)}{0.25}=24.8125)
中心差分：(f'(3)=27.0625)

可以看到中心差分最为准确，且两点间距变小时，差分计算会更为准确。
MATLAB实现：

function dx = derivative(x,y)
% derivative calculates the derivative of a function that is given by a set
% of points. The derivative at the first and last points are calculated by
% using the forward and backward finite difference formula, respectively.
% The derivative at all the other points is calculated by the central
% finite difference formula.
% Input variables:
% x  A vector with the coordinates x of the data points.
% y  A vector with the coordinates y of the data points.
% Output variable:
% dx  A vector with the value of the derivative at each point.

n = length(x);
dx(1)=(y(2)-y(1))/(x(2)-x(1));
for i=2:n-1
    dx(i)=(y(i+1)-y(i-1))/(x(i+1)-x(i-1));
end
dx(n)=(y(n)-y(n-1))/(x(n)-x(n-1));

2. 泰勒公式有限差分法

2.1 一阶微分（两点法）

• 前向展开：
  (f(x))在点(x_{i+1})处的值可以使用如下泰勒公式来逼近：$$f(x_{i+1})=f(x_i)+f'(x_i)h+frac{1}{2!}f''(x_i)h^{2+frac{1}{3!}f'''(xi_1)h}3$$ 其中(h=x_{i+1}-x_i),(xi_1)是(x_{i+1})与(x_i)之间的数。求解该公式，有：$$f'(x_i)=frac{f(x_{i+1})-f(x_i)}{h}-frac{1}{2!}f''(x_i)h-frac{1}{3!}f'''(xi_1)h^2=frac{f(x_{i+1})-f(x_i)}{h}+O(h)$$ 等同于之前的前向差分，具有一阶准确度

• 后向展开：
  (f(x))在点(x_{i-1})处的值可以使用如下泰勒公式来逼近：$$f(x_{i-1})=f(x_i)-f'(x_i)h+frac{1}{2!}f''(x_i)h^{2-frac{1}{3!}f'''(xi_2)h}3$$ 其中(h=x_{i}-x_{i-1}),(xi_1)是(x_{i-1})与(x_i)之间的数。求解该公式，有：$$f'(x_i)=frac{f(x_{i})-f(x_{i-1})}{h}+frac{1}{2!}f''(x_i)h-frac{1}{3!}f'''(xi_2)h^2=frac{f(x_{i})-f(x_{i-1})}{h}+O(h)$$ 等同于之前的后向差分，具有一阶准确度

• 中心展开(假设间距相同)
  结合上述两种展开，可以得到：$$f(x_{i+1})-f(x_{i-1})=2hf'(x_i)+frac{1}{3!}(f'''(xi_1)+f'''(xi_2))h^3$$ 因此有 $$f'(x_i)=frac{f(x_{i+1})-f(x_{i-1})}{2h}+O(h^2)$$ 等同于之前的中心差分，可以看到，具有二阶准确度

2.2 一阶微分（三点法）

分别写出(x_{i-2},x_{i-1},x_{i+1},x_{i+2})四点的泰勒展开：

• (f(x_{i+1})=f(x_i)+f'(x_i)h+frac{1}{2!}f''(x_i)h^2+frac{1}{3!}f'''(xi_1)h^3)

• (f(x_{i-1})=f(x_i)-f'(x_i)h+frac{1}{2!}f''(x_i)h^2-frac{1}{3!}f'''(xi_2)h^3)

• (f(x_{i+2})=f(x_i)+f'(x_i)2h+frac{1}{2!}f''(x_i)(2h)^2+frac{1}{3!}f'''(xi_3)(2h)^3)

• (f(x_{i-2})=f(x_i)-f'(x_i)2h+frac{1}{2!}f''(x_i)(2h)^2-frac{1}{3!}f'''(xi_4)(2h)^3)

可以看到：$$f(x_{i+2})-4f(x_{i+1})=-3f(x_i)-2f'(x_i)h+frac{f'''(xi_3)}{3!}(2h)^{3-frac{4}{3!}f'''(xi_1)h}3$$ 进一步得：$$f'(x_i)=frac{-3f(x_i)+4f(x_{i+1})-f(x_{i+2})}{2h}+O(h^2)$$ 这是一阶微分的三点前向公式，具有二阶准确度，类似地，可以得到如下具有二阶准确度的三点后向公式：$$f'(x_i)=frac{f(x_{i-2})-4f(x_{i-1})+3f(x_i)}{2h}+O(h^2)$$

2.3 二阶微分（三点法）

注意：

• (f(x_{i+1})=f(x_i)+f'(x_i)h+frac{1}{2!}f''(x_i)h^2+frac{1}{3!}f'''(x_i)h^3+frac{1}{4!}f^{(4)}(xi_1)h^4)

• (f(x_{i-1})=f(x_i)-f'(x_i)h+frac{1}{2!}f''(x_i)h^2-frac{1}{3!}f'''(x_i)h^3+frac{1}{4!}f^{(4)}(xi_2)h^4)

• (f(x_{i+2})=f(x_i)+f'(x_i)2h+frac{1}{2!}f''(x_i)(2h)^2+frac{1}{3!}f'''(x_i)(2h)^3+frac{1}{4!}f^{(4)}(xi_3)(2h)^4)

• (f(x_{i-2})=f(x_i)-f'(x_i)2h+frac{1}{2!}f''(x_i)(2h)^2-frac{1}{3!}f'''(x_i)(2h)^3+frac{1}{4!}f^{(4)}(xi_4)(2h)^4)

前面两式相加，得 $$f(x_{i+1})+f(x_{i-1})=2f(x_i)+f''(x_i)h^{2+frac{1}{4!}f}{(4)}(xi_1)h^{4+frac{1}{4!}f}{(4)}(xi_2)h^4$$ 可得$$f''(x_i)=frac{f(x_{i+1})-2f(x_i)+f(x_{i-1})}{h^2}+O(h2)$$ 此即为三点中心差分公式，具有二阶准确度。类似地，推导可得如下五点中心差分公式：$$ f''(x_i)=frac{-f(x_{i-2})+16f(x_{i-1})-30f(x_i)+16f(x_{i+1})-f(x_{i+2})}{12h^2}+O(h4) $$ 具有四阶准确度。

另一方面：$$f(x_{i+2})-2f(x_{i+1})=-f(x_i)+f''(x_i)h^{2+frac{6f'''(x_i)}{3!}h}3+frac{1}{4!}f^{{(4)}(xi_3)(2h)}4-frac{2}{4!}f^{(4)}(xi_1)h4$$ 可得 $$f''(x_i)=frac{f(x_i)-2f(x_{i+1})+f(x_{i+2})}{h^2}+O(h)$$ 此即为三点前向差分公式，具有一阶准确度，类似可得如下具有一阶准确度三点后向差分公式： $$f''(x_i)=frac{f(x_{i-2})-2f(x_{i-1})+f(x_{i})}{h^2}+O(h)$$

MATLAB实现：

function [yd,ydd] = FirstScndDerivPt(x,y)
n = length(x);
h = x(2)-x(1);
%  首个数据，一阶导数使用三点前向差分,二阶导数使用4点前向差分
yd(1) = (-3*y(1)+4*y(2)-y(3))/(2*h);
ydd(1) = (2*y(1)-5*y(2)+4*y(3)-y(4))/(h^2);

% 中间数据，一阶导数使用两点中心差分，二阶导数使用三点中心差分
for i=2:n-1
    yd(i)=(y(i+1)-y(i-1))/(2*h);
    ydd(i)=(y(i-1)-2*y(i)+y(i+1))/(h^2);
end

% 末尾数据，一阶导数使用三点后向差分，二阶导数使用4点后向差分
yd(n) = (y(n-2)-4*y(n-1)+3*y(n))/(2*h);
ydd(n) = (-y(n-3)+4*y(n-2)-5*y(n-1)+2*y(n))/(h^2);
figure
subplot(3,1,1)
plot(x,y)
subplot(3,1,2)
plot(x,yd)
subplot(3,1,3)
plot(x,ydd)
end

3. 拉格朗日多项式求导公式

对点((x_i,y_i),(x_{i+1},y_{i+1}),(x_{i+2},y_{i+2}))进行拉格朗日多项式插值，有 $$f(x)=frac{(x-x_{i+1})(x-x_{i+2})}{(x_i-x_{i+1})(x_i-x_{i+2})}y_i+frac{(x-x_{i})(x-x_{i+2})}{(x_{i+1}-x_{i})(x_{i+1}-x_{i+2})}y_{i+1}$$ $$+frac{(x-x_{i})(x-x_{i+1})}{(x_{i+2}-x_{i})(x_{i+2}-x_{i+1})}y_{i+2}$$
此时 $$f'(x)=frac{2x-x_{i+1}-x_{i+2}}{(x_i-x_{i+1})(x_i-x_{i+2})}y_i+frac{2x-x_{i}-x_{i+2}}{(x_{i+1}-x_{i})(x_{i+1}-x_{i+2})}y_{i+1}$$ $$+frac{2x-x_{i}-x_{i+1}}{(x_{i+2}-x_{i})(x_{i+2}-x_{i+1})}y_{i+2}$$ 因此$$f'(x_i)=frac{2x_i-x_{i+1}-x_{i+2}}{(x_i-x_{i+1})(x_i-x_{i+2})}y_i+frac{x_{i}-x_{i+2}}{(x_{i+1}-x_{i})(x_{i+1}-x_{i+2})}y_{i+1}$$ $$+frac{x_{i}-x_{i+1}}{(x_{i+2}-x_{i})(x_{i+2}-x_{i+1})}y_{i+2}$$ 当(x_{i+1}-x_i=x_{i+2}-x_{i+1}=h)时，有 $$f'(x_i)=frac{-3y_i+4y_{i+1}-y_{i+2}}{2h}$$ 此式与三点前向差分公式一致。

4. 数值偏微分

对二元函数(f(x,y)) ，其在点((a,b))处的偏微分定义为：$$frac{partial f}{partial x}|{(a,b)}=lim{x ightarrow a} frac{f(x,b)-f(a,b)}{x-a}$$ $$frac{partial f}{partial y}|{(a,b)}=lim{x ightarrow b} frac{f(a,y)-f(a,b)}{ y-b}$$ 对一阶偏微分(frac{partial f}{partial x})，(frac{partial f}{partial y})，两点前向公式为：$$frac{partial f}{partial x}|{(x_i,y_i)}=frac{f(x{i+1},y_i)-f(x_i,y_i)}{x_{i+1}-x_i}$$ $$frac{partial f}{partial y}|{(x_i,y_i)}=frac{f(x{i},y_{i+1})-f(x_i,y_i)}{y_{i+1}-y_i}$$ 两点后向公式为：$$frac{partial f}{partial x}|{(x_i,y_i)}=frac{f(x{i},y_i)-f(x_{i-1},y_i)}{x_{i}-x_{i-1}}$$ $$frac{partial f}{partial y}|{(x_i,y_i)}=frac{f(x{i},y_{i})-f(x_i,y_{i-1})}{y_{i}-y_{i-1}}$$
两点中心差分公式：$$frac{partial f}{partial x}|{(x_i,y_i)}=frac{f(x{i+1},y_i)-f(x_{i-1},y_i)}{2h_x}$$ $$frac{partial f}{partial y}|{(x_i,y_i)}=frac{f(x{i},y_{i+1})-f(x_i,y_{i-1})}{2h_y}$$ 对二阶偏微分(frac{partial^2 f}{partial x^2},frac{partial^2 f}{partial y^2})，三点中心差分公式为：$$frac{partial^2 f}{partial x^{2}|_{(x_i,y_i)}=frac{f(x_{i-1},y_i)-2f(x_i,y_i)+f(x_{i+1},y_i)}{h}2_x}$$ $$frac{partial^2 f}{partial y^{2}|_{(x_i,y_i)}=frac{f(x_{i},y_{i-1})-2f(x_i,y_i)+f(x_{i},y_{i+1})}{h}2_y}$$ 对二阶偏微分(frac{partial^2 f}{partial x partial y}=frac{partial}{partial x}(frac{partial f}{partial y})=frac{partial}{partial y}(frac{partial f}{partial x}))，逐步计算即可：$$frac{partial^2 f}{partial x partial y}|{(x_i,y_i)}=frac{frac{f(x{i+1},y_{i+1})-f(x_{i-1},y_{i+1})}{2h_x}-frac{f(x_{i+1},y_{i-1})-f(x_{i-1},y_{i-1})}{2h_x}}{2h_y} ext{((frac{partial f}{partial x})取中心差分)}$$
MATLAB实现：

function [dfdx,dfdy] = ParDer(x,y,f)
n = length(x);
m = length(y);
hx = x(2)-x(1);
hy = y(2)-y(1);

% 首位数据使用三点前向
for j = 1:m
    dfdx(1,j) = (-3*f(1,j)+4*f(2,j)-f(3,j))/(2*hx);
end 
for i = 1:n
    dfdy(i,1) = (-3*f(i,1)+4*f(i,2)-f(i,3))/(2*hy);
end

% 两点中心差分
for i = 2:n-1
    for j = 1:m
        dfdx(i,j) = (f(i+1,j)-f(i-1,j))/(2*hx);
    end
end

for j = 2:m-1
    for i = 1:n
        dfdy(i,j) = (f(i,j+1)-f(i,j-1))/(2*hy);
    end
end

% 末尾数据使用三点后向
for j = 1:m
dfdx(n,j) = (f(n-2,j)-4*f(n-1,j)+3*f(n,j))/(2*hx);
end

for i = 1:n
    dfdy(i,m) = (f(i,m-2)-4*f(i,m-1)+3*f(i,m))/(2*hy);
end

end

5. Richardson外推加速算法

Richardson外推加速算法能够把两个低精度的方法组合成一个高精度的计算方法，假设$$D=D(h)+k_2h^2+k_4h4$$是一种数值微分计算方法，(D(h))是估计的微分，(k_2h^2)和(k_4h^4)是估计误差，可以看到，具有二阶精度，如果把间距调整为(frac{h}{2})，则$$D=D(frac{h}{2})+k_2(frac{h}{2})^{2+k_4(frac{h}{2})}4$$ 其精度仍是二阶的。但是有：$$4D-D=4D(frac{h}{2})-D(h)-frac{3}{4}k_4h^4$$ 因此，$$D=frac{1}{3}(4D(frac{h}{2})-D(h))+O(h^4)$$ 可以看到，具有四阶精度。

举个例子，考虑一阶微分的两点中心差分法（二阶精度）：$$f'(x_i)=frac{f(x_{i+1})-f(x_{i-1})}{2h}+frac{1}{3!}f'''(x_i)h^2+O(h4)$$ $$=frac{f(x_{i}+h)-f(x_{i}-h)}{2h}+frac{1}{3!}f'''(x_i)h^2+O(h4)$$ 缩短间距，有 $$f'(x_i)=frac{f(x_{i+1})-f(x_{i-1})}{h}+frac{1}{3!}f'''(x_i)(frac{h}{2})^2+O(h4)$$ $$=frac{f(x_{i}+frac{h}{2})-f(x_{i}-frac{h}{2})}{h}+frac{1}{3!}f'''(x_i)(frac{h}{2})^2+O(h4)$$ 按照Richardson外推加速算法，有 $$ 4f'(x_i)=4frac{f(x_{i}+frac{h}{2})-f(x_{i}-frac{h}{2})}{h}+frac{4}{3!}f'''(x_i)(frac{h}{2})^2+4O(h4) $$ 进一步 $$ 3f'(x_i)=[4frac{f(x_{i}+frac{h}{2})-f(x_{i}-frac{h}{2})}{h}-frac{f(x_{i}+h)-f(x_{i}-h)}{2h}]+3O(h^4) $$ 因此 $$ f'(x_i)=frac{1}{3}[4frac{f(x_{i}+frac{h}{2})-f(x_{i}-frac{h}{2})}{h}-frac{f(x_{i}+h)-f(x_{i}-h)}{2h}]+O(h^4) $$ 可以看到，精度提高到了四阶。

6. 总结

本节课主要介绍了一些数值微分算法，对于一阶微分，最常用的有两点前向差分（精度为(O(h))）、两点后向差分（精度为(O(h))）以及两点中心差分算法（精度为(O(h^2))），其表达式均可以通过处理泰勒展开式来得到。通过处理泰勒展开式还可以得到一阶微分的三点前向差分和三点后向差分算法，精度与两点中心差分一致。对于二阶微分，同样可以利用泰勒展开，得到三点前向差分、三点后向差分以及三点中心差分算法。另一方面，还可以通过拉格朗日插值公式，得到相应的微分计算公式。这些公式又都可以很容易推广到多元函数的数值微分中去。最后，对于两个精度不高的微分算法，可以通过Richardson外推加速算法得到一个精度更高的算法，在实际问题中具有很广泛的应用。