数值计算day7-数值微分

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

1. 有限差分法

给定一个函数f(x)f(x)x=a处的微分f'(x)定义为:\frac{df(x)}{dx}|_{x=a}=f'(a)=\lim_{x\rightarrow} \frac{f(x)-f(a)}{x-a} 图上的解释是f(x)x=a处的斜率:

image.png

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

  • 前向差分:\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}}

image.png

例: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_1x_{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_1x_{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(h^2) 此即为三点中心差分公式,具有二阶准确度。类似地,推导可得如下五点中心差分公式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(h^4) 具有四阶准确度。

另一方面: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)h^4 可得 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. 数值偏微分

image.png

对二元函数 ,其在点处的偏微分定义为: 对一阶偏微分,,两点前向公式为: 两点后向公式为:
两点中心差分公式: 对二阶偏微分,三点中心差分公式为: 对二阶偏微分,逐步计算即可:
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_4h^4是一种数值微分计算方法,D(h)是估计的微分,k_2h^2k_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(h^4) =\frac{f(x_{i}+h)-f(x_{i}-h)}{2h}+\frac{1}{3!}f'''(x_i)h^2+O(h^4) 缩短间距,有 f'(x_i)=\frac{f(x_{i+1})-f(x_{i-1})}{h}+\frac{1}{3!}f'''(x_i)(\frac{h}{2})^2+O(h^4) =\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(h^4) 按照Richardson外推加速算法,有 4f'(x_i)=4\frac{f(x_{i}+\frac{h}{2})-f(x_{i}-\frac{h}{2})}{h}+\frac{4}{3!}f'''(x_i)(\frac{h}{2})^2+4O(h^4) 进一步 3f'(x_i)=[4\frac{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}[4\frac{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外推加速算法得到一个精度更高的算法,在实际问题中具有很广泛的应用。

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 205,386评论 6 479
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 87,939评论 2 381
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 151,851评论 0 341
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 54,953评论 1 278
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 63,971评论 5 369
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 48,784评论 1 283
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 38,126评论 3 399
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,765评论 0 258
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 43,148评论 1 300
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,744评论 2 323
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,858评论 1 333
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,479评论 4 322
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 39,080评论 3 307
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,053评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,278评论 1 260
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 45,245评论 2 352
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 42,590评论 2 343

推荐阅读更多精彩内容

  • 何谓数值分析? 众所周知现实生活中科学技术上面的问题大多数都能够通过建立对应的数学模型把实际问题转化成一个数学问题...
    罗泽坤阅读 8,774评论 0 14
  • 计算机配色发展了许多年,但仍不能满足使用需求,归根结底是由于其理论存在一些问题,使得其应用受到一定限制。如果可以找...
    申申申申申申阅读 1,806评论 0 2
  • 最近一段时间在复习数值计算相关内容,也恰逢简书断更了,不用每天督促着自己非要更出点什么东西才好,有更多的空间来打磨...
    抄书侠阅读 1,238评论 0 6
  • 写在前面:关于为何要写这个笔记? 关于学习与遗忘,在考完这门课后,我还能记得些什么呢?引用mbinary的文章:h...
    冬风十里Y阅读 1,399评论 0 1
  • 概率论与数理统计 无穷小阶数 无穷小量表述:线性逼近 相当于利用切线和斜率来理解误差和逼近。 泰勒级数:线性逼近 ...
    Babus阅读 808评论 0 1