数值分析18(数值微分)

1、1/32导数的数值计算方法导数的数值计算方法数值求导的外推方法数值求导的外推方法数值分析 1910:192/32重温微积分重温微积分微分微分(Differentiation)积分积分(Integration) 微分与积分构成了一对互逆的运算微分与积分构成了一对互逆的运算 ( ) ( ) dF xF xC ( )f x dxIntegrals as Sums and Derivatives as Difference凡線面體皆設為由小漸大凡線面體皆設為由小漸大,一剎那中所增之積即微分也。一剎那中所增之積即微分也。其全積即積分也其全積即積分也10:193/32重温微积分重温微积分微积分中蕴含的对立

2、统一思想微积分中蕴含的对立统一思想( ) ( )( ) ( )( )( )u x v xu x v x dxv x u x dx( ( ) ( )( ) ( )( ) ( )u x v xu x v xu x v x( )( )( )f bf afba ( ) ( )( )( )bbaaf x g x dxfg x dx ( )( )( )baf x dxF bF a 微微积积分分基基本本定定理理10:194/32Integrals as Sums and Derivatives as Difference x x1 x2 xm y y1 y2 ym如何计算导数如何计算导数(微分微分)或者梯度

3、或者梯度?3tan( )=sin (cosh )xf xxx函数复杂或仅仅给定函数复杂或仅仅给定离散的观察数据离散的观察数据(函数值函数值)？10:195/32回顾回顾0()( )( )limhf xhf xfxh 导数概念是精确刻画事物变化率的有力工具。导数概念是精确刻画事物变化率的有力工具。10:19 Science is knowledge which we understand so well that we can teach it to a computer; and if we dont fully understand something, it is an art to de

4、al with it. 计算机程序设计的艺术计算机程序设计的艺术 Donald Knuth(图灵奖得主图灵奖得主)6/32回顾回顾: : Taylor公式公式200000( )(1)0100()( )()()()()2!() ( )()(1)!()( )!nnnnfxf xf xfxxxxxfxxxnfxxn 20000(1)10()()()()2!() ( )!()(1)nnnnfxfxhfxfxhhfxhfnnh 10:197/32一阶前向差分公式一阶前向差分公式2()( )( )( )2hf xhf xhfxf ()( )( )=( )2f xhffxfxhh (),nOnh 前前向向差

5、差分分公公式式是是近近似似一一阶阶导导数数的的一一阶阶方方法法。一一般般地地如如果果误误差差是是我我们们就就称称公公式式是是 阶阶近近似似。0,0,( ),0( )( ),hhhhfxhffx 称称这这个个公公式式一一阶阶的的微微妙妙之之处处是是 与与 有有关关。一一阶阶的的概概念念是是当当时时误误差差应应正正比比于于 。当当时时是是移移动动的的 因因此此比比例例常常数数改改变变了了。但但只只要要连连续续 当当时时比比例例常常数数趋趋于于这这就就使使得得称称公公式式为为一一阶阶的的是是恰恰当当的的。10:198/32 两点中心差分公式两点中心差分公式231( )( )()()( )2+6ff

6、xhfxfxf xhhh 2( )()( )6=2f xhfhxfxfhh 232( )( )()()(2)6ff xhfxfxf xhhh 二阶导数近似二阶导数近似(4)2413( )( )()( )()(24+26)ffxf xhf xhfxfx hhh 22( )12()2 ( )()( )=f xhf xf xhhfhfx （4 4）23(4)24()( )()( )224( )( )6ffxf xhf xhfxfhx hh 10:199/32回顾回顾 介值定理介值定理 (Intermediate Value Theorem)10:19设函数设函数y=f(x)在闭区间在闭区间a,b上连

7、续上连续, ,且在这区间端且在这区间端点处取值不同时点处取值不同时, ,即：即：f(a)=A,f(b)=B,且且AB。那么。那么不论不论C是是A与与B之间的怎样一个数之间的怎样一个数, ,在闭区间在闭区间a,b内内至少有一点至少有一点, ,使得使得 f()=C 。10/32例例1. 取取h=0.1, ,分别用两点前向差分公式和两点中心差分别用两点前向差分公式和两点中心差分公式近似分公式近似f(x)=1/x在在x=2处的导数处的导数。112.12()( )( )=0.23810.1f xhf xfxh 2( )=,(2)0.2500fxxf 112.11.9()()( )=0.250620.2f

8、 xhf xhfxh 2( )()( )6=2f xhfhxfxfhh ()( )( )=( )2f xhffxfxhh 10:190.0108=0.0119= 0.01250.0005= 0.0006 = 0.000711/3212/32舍入误差舍入误差到目前为止到目前为止,所有公式都破坏了不要进行相近数所有公式都破坏了不要进行相近数相减的规则。这对于数值微分是一个极大的困相减的规则。这对于数值微分是一个极大的困难难,但是它在本质上式不可能避免的。但是它在本质上式不可能避免的。例例2. 求求f(x)=ex 在在x=0处导数的近似。处导数的近似。 h 前向差分前向差分 误差误差 中心差分中心差

9、分 误差误差10:1913/32例例3. 基于中心差分公式的研究更高阶近似公式基于中心差分公式的研究更高阶近似公式(2245)()()( )= ( )5!)2(6(hf xhf xhhF hfxfxhfx 235(54( )4( )+( )( )65()( )( )4!+( )2!hhf xhf xfxhhhffxfxxfx 24(4)35(5)()( )( )+( )2)( )( )5!4!6hhf xhf xfxhhhfxfxfxxf 4(5)2422(/ 2)(/ 2)( / 2)( )( )5! 2( )6 2hff xhf xhhF hfxfxxh 2224( )( )444 1 =

10、 ( )() hFFhFfxO h 10:19422 ( )4/ 3( / 2)1/ 3( )F hF hF h14/32松弛思想松弛思想 目标值目标值Q有两个精度相当的近似值有两个精度相当的近似值F1和和F2,如如果将这两个近似值加工成更高精度的结果呢果将这两个近似值加工成更高精度的结果呢?改善精度的一种简便而有效的办法是改善精度的一种简便而有效的办法是,取两者的取两者的某种加权平均值作为改进值某种加权平均值作为改进值,即令即令12112 (1() ) QFFQFFF 或或适当选取平均化系数适当选取平均化系数 调整校正量调整校正量 以将以将F1加工成某个更高精度结果。这种基于校加工成某个更高

11、精度结果。这种基于校正量的调整或松动的方法称之为松弛方法。正量的调整或松动的方法称之为松弛方法。 21()FF 10:1915/3211 ( )+ nnnnnnQnFQF hchch 近近似似给给定定的的量量 的的 阶阶公公式式24(5)()()( )( )( )265!f xhf xhhhfxfxfxh例例如如22( )( )21 Richardson nhnnnFFhQ 外外推推 2 / (21),11/ (21)nnn 外外推推2 ( )(1)( ) hnnQFF h10:1911 ( )+nnnnncQF hchh 111 ( / 2)+/ 2/2nnnnnnnc hQF hch 22

12、( )( )1211+ +nhnnnnFFhnQdh 16/32例例4. 推导中心差分公式的推导中心差分公式的Richardson外推公式外推公式2()()( )=( )26f xhf xhhfxfh 22( )( )21 Richardson nhnnnFFhQ 外外推推2222222222( )( )421()()()()2() 8() 8()()62, = 4/ 3 = hhhhhFFhf xf xf x hf x hhhf x hf xf xf x hhnF 10:19是否可以进一步地外推?17/32例例5. 用Richardson外推公式计算外推公式计算f(x)=x2e-x在在x=0

13、.5的导数。的导数。2(1/2)2(1/2)1( )(1/ 2)(1/ 2)2hhF hh eh eh解解10:190.1,0.05,0.025,h 当当时时0.450.45( )16049081( / 2)0761693( / 440.4)692624885F hF hF h 11( )0.45489(/ 2)0992.453181 591 248FhFh 2( )0.454897994Fh 0.454897994(0.5)f 精精 确确 值值18/321. n阶差分阶差分, diff(X,N,DIM) X = 3 7 5; 0 9 2, diff(X,1,1), diff(X,1,2)Ma

14、tlab微分微分I=imread(lena_color_512.tif);imshow(diff(I),)2. 符号微分符号微分diffsyms x; f=(sin(xtan(x)*cosh(x)3; f1=diff(f,1)10:1919/323. gradient x,y,z=peaks;contour(x,y,z,20)%contour plotsMatlab微分微分10:19 %Vector fieldsn=-2:.2:2;x,y,z=peaks(n);contour(x,y,z,10);U,V=gradient(z,.2);hold on,quiver(x,y,U,V)20/32应用

15、应用 1: : Image Stitching 21/32应用应用 2: Colorization22/32应用应用 3: High Dynamic Range23/32应用应用 4: Image Editting 24/32Adobe Photoshop: Healing Brush Todor Georgiev, senior research scientist at Adobe, asked me, Jeff, do you want to work on something interesting? We started to work on a scratch remover. T

16、hen we talked for days and days and days, and finally said, This will be solving a heat equation, right? And all of a sudden, I realized that it was something I knew when I was back in schoolI took a heat transfer class. In plain English, if you put hot rod and a cold rod together, what happens? The

17、 heat transfers from one to the other. That equation explains how the heat transfers, and its a very smooth diffusion. The same thing happens when you take a piece of an image from one to anotherit diffuses. It is like the pieces have different temperatures. You can see how they start to blend, and

18、that process is very smooth, very seamless. When you apply that equation to your image, you can fool your eye. You dont see the transition; its all smooth. In the imaging community, this essential goal is addressed by solving the Poisson equation. Jenchan Chien Adobe Photoshop首席科学家首席科学家25/3226/32Ima

19、ge Editing问题的数学模型问题的数学模型2*min s.t. |uuvuu 2222*( , )div , ( , )( , )|( , )uuxyx yvx yu x yux y 上述模型极小值点是上述模型极小值点是Poisson 方程方程 Dirichlet问题的解问题的解: :27/32( ),01(0)(1)0uf xxuu 一维一维Poisson方程方程 h = 1/(n+1) , xi= ih (i= 0,1, , n+1) ui= u(xi ), (i = 0,1, , n+1)121)2(iiiih f xuuu (4)2413( )( )()( )()(24+26)u

20、uxu xhu xhuxu x hhh 22( )12()2 ( )()( )=u xhu xu xhhuhux （4 4）23(4)24()( )()( )224( )( )6uuxu xhu xhuxuhx hh 28/32(1)( )( )211()/ 2ikkkiiiuuuh f x 10:1912122211222-1()-12-1()-12-1()-12()nnnnh fuuuxuh f xh f xh f x 迭代方法迭代方法: :直接方法直接方法: : 反斜线算符反斜线算符; 特别地对于稀疏矩阵更加高效特别地对于稀疏矩阵更加高效(1)(1)( )211()/ 2ikkkiiiuuuh f x 29/322222( , ),01,01(0, )(1, )( ,0)( ,1)0uuxyf x yxyuyuyu xu x 二维二维Poisson方程方程 h = 1/(n+1) , xi= ih (i = 0,1, , n+1)h = 1/(n+1) , yi= ih (i = 0,1, , n+1) ui,j= u(xi,yj ), (i = 0,1, , n+1, j = 0,1, , n+1)1,1,1,21(,4)iji ji ji jijiih f xy