数值微分的计算方法

1、数值微分的计算方法内容摘要 求解数值微分问题，就是通过测量函数在一些离散点上的值，求得函数的近似导数。本文就所学知识，归纳性地介绍了几种常用的数值微分计算方法。并举例说明计算，实验结果表明了方法的有效性。关键词数值微分 Taylor 展开式Lagrange 插值三对角矩阵引言：数值微分即根据函数在一些离散点的函数值，推算它在某点的导数或高阶导数的近似值的方法。常见的可以用一个能够近似代替该函数的较简单的可微函数(如多项式或样条函数等)的相应导数作为能求导数的近似值，由此也可导出多点数 值微分计算公式。当函数可微性不太好时，利用样条插值进行数值微分要比多项式插 值更适宜。I. Taylor 展开

2、式方法理论基础：Taylor展开式2nf(X)=f(Xo 广(X冷/ J )+2f"(X。)*+fCX°我们借助Taylor展开式，可以构造函数 f (x)在点x= X。的一阶导数和二阶导数的数值微分公式。取步长 h。则(1)(2)(3)(4)f (Xo h) = f (x。) hf'(x。) h f ( i) i (X。，X。 h) 2所以f (Xo h)f (Xo) h .f (X。)- f ( i) i(号。h)h2同理2 (Xo h,X。)hf (Xo -h) =f (X。)-hf (X。)万 f ( 2)f'(5。)-(尚-h) §f&#

3、39;'(侦 2 (X。一h,X。)h2式 和式(4)是计算f'(X。)的数值微分公式，其截断误差为0(h),为提高精度，将Taylor展开式多写几项f(xo h) = f(Xo) hf2 3 4(X。)3 f (X。)X f (X。)板 f ( 1)2624, h2f (X。-h) = f (X。)一hf (X。)23 4f (X。)- ： f (X。)h f(4)( 2)6241 (乂。，*。 h)2 (X。-h, X。)两式相减得-2(5)f ( x h) - f(Xnh) h'"4f (X。)= f(X° h) f(X° h)_ L

4、 f (xo) 0(h )2h62上式为计算f (X。)的微分公式，其截断误差为 0(h),比式(2)和(4)精度寄两式相加，如果f (4)(X)ECX。一 h,x。+h,则有2(6)f''(x)_f(x)h)_2f(X0)f(X0_h) h f(4)( f(Xo) 一h2一12 f J)E(xoh,x°+h)2式(6)是计算f (Xo)的数值微分公式，其截断误差为O(h2)。例1 .设函数f(x)= lnx,X0= 2,h= 0.1，试用数值微分公式计算f (2)的值。解 由式(2)、式和式(5)分别计算结果为ln 2.1ln 2f (2) 0.48790.1ln

5、2 - ln 1.9 - 一f (2)：、: 0.51290.1ln2.1_ln1.9f (2) 0.50040.2与真值f'(2)= 0.5相比，式(5)计算的结果精度较高。2. 数值微分的Lagrange插值方法设函数y= f (x)具有n+ 1个实验数据：(为，f (xi)(i = 0,1,2，n),我们希望估讨'，、f (x)的值，特别x =x时，估计f (x)的值。基于插值方法的数值微分做法是，由已知(x, f (x)(i = 0,1，,n)建立Lagrange插值多项式或Newton插值多项式(这里以 Lagrange插值方法为例)，即f (x)= Ln (x)+

6、尺(x)于是f x = Ln xRn xx= x时，有f (xi) =Ln(xi)Rn(xi)(i =0,1,2,., n)其中'f (n 1)( )，Rn(x) = n1(xi)(n 1)!略去误差项有f (x) ： Ln(xi)实际运用中，等距节点更为常见。设, b- ah =，为=na + ih (i = 0,1,2，,n), x = a + th ,于是有f (x) = y° 十；侦 +")M 十.EaM + Rn(X)所以d , tf (为)=?y0y0 -.-1!()n(1 -j)(n 1)!j=0jTdtn (n 1) h fn+ Ktf.Jt n 1

7、) ,：ny0tgn!3. 多点数值微分公式由于高阶插值的不稳定性, 型求导公式。(1)两点公式(n=1)实际应用时多采用n=1,2,4的两点、三点和五点等多点插值f'(X0) =：(yi -y0)h、1 ,、f (X1) =(y -y°) h_ ,h _ R (X0) 2 f ()R(X1)业 f''()2(2) 三点公式(n=2),',、1 ，、f (冷)=二(3y° 4乂 y?)2h一，,、1 .f (X1) = = (y。+y2)2h,',1f (X2) =(y° 4山 +3y2)2hR(X0)=§f(3)

8、()R(X1) =-项 f(3)()62R(X0) = §f(3)()3(8)五点公式(n=4)r,',、1,f (x°)(-25y° 48y -36y2 16y3 -3y4)12h.',、1,f (X1)(-3y0 -10y1 18y2 -6y3 - y4)12h，、1,f (X2)(y° -8y1 8y3 -y4)12h,'1f (X3)(T0 6贝18y2 10y3 3y4)12h，、1,f (X4)(3y° -16y" 36y48y3 25y4)12hR(X0)=gf(5)()5R(xJA* )20R(

9、X2h：f(5)() 30R(X3)=-£f(5)()20R(X°)=£f(5)()5(9)例2.设f (x) = ln x，取h = 005 ,分别用三点公式和五点公式计算f' (2)的近似值。解：由式(8)有1f (2). .(-3f (2) 4 f (2.05) f (2.10) =0.4998028612 0.05f'(2) =1(f (1.95) - f (2.05) =0.5001042052 0.05 1f (2). . (f (1.90) -4f (1.95) - 3f (2) =0.4997793762 0.05由式(9)有1f

10、(2) =1(f (1.90) _8 f (1.95) 6f (2.05) _f (2.10) =0.49999984312 0.05t . _f. t . 一' 一、_. .与真值f (2)= 0.5相比，三点公式已有相当满意精度，而五点公式的结果是十分满意的。4. 数值微分的隐式格式前述的Taylor 展式法和数值微分格式均称为显式格式，即直接由已知的 ， 、 . . '， 、f (Xi)(i = 0,1,2，n)，经过适当的算术四则运算，立即可得 f (X)的近似值。显式格式优点是计算方便，工作量小，缺点是数值不稳定。为克服后一缺点，隐式格式常常具有数值稳定性。数值微分的

11、隐式格式建立方法常用通过Taylor展开式方法或数值积分方法等不同途径。首先我们用Taylor展开式方法来推导数值微分的隐式格式。由式和式(6)，我们用&代替X0得：1h'”一 4f(XjF皿 h)-f(Xk 顶-百 f ("0(h)f''(Xk) =f(Xk h) -2f (Xk) f(Xk h) f (4)()h212所以也有'''1'''h2(5)-f同二仍心为-小口丛叫-"？()将最后f (xk)表达式代入f (Xk)表达式可得114f (Xk)Fg 枷-皿-h)飞仃(Xk +h) -2

12、f (f (Xk -h) +o(h)略去误差项o(h4)，并且mk表示f (Xk)的近似值，则有,3 、mk 1 - 4mk mk 1 f (Xk 1) - f (Xk 盘(k =1,2,., n)h同样，用数值积分方法也可推导数值微分的隐式格式。Xk 1Xk ，f(X)dX=f'(Xk )dX(k = 1,2,., n)'、X _Lbb -aa b.根据 simpson 公式f(X)dX =f (a) +4f( 一)+ f (b)有 a 、，622hf (Xk 1) - f (Xk)=;mk4m" mk 1, (k = 1,2,., n)6,3，mk4mk mk 1

13、 =f g 1) - f (Xj)(k =1,2,., n)h上式是关于n+ 1个未知量m0,m,mn的n- 1个方程，如果已知 m0 =3 f (X°),mn = f (Xn)，记 dk = f (Xf) f (Xj),故有万程组：h4 1一 m1-d m°1 4 1m2d2141m3d3+ + +a， > ， 14 1mudn:. 1 4 一：mn一idn-mn_(10)这就是求m1,典，mn-1的线性方程组。由于系数矩阵是严格对角占优的三对角矩阵，因此非奇异，解存在且唯一，可由追赶法求解，且数值稳定。例3.设函数f (x)= lnx在节点x= 1.5, 1.6,

14、 1.7, 1.8, 1.9上的函数值，及-'. f (1.5), f (2.0)如表1所示。试用数值微分隐式格式式(10)和式(11)求出相应节点上一阶导数的近似值。表1例3节点数据xf(x)f ' (x)1.50.4054651080.6666666671.60.4700036291.70.5306282511.80.5877866641.90.6418538862.00.6931471820.500000000解：由式(10)有面m21 m3七邮一3.75489429 -0.6666666673.533491053.33676905316081554 -0.5000000

15、00-解方程组得:m =0.62499828611483m2 =0.58823447854067ms =0.55555484972249m4 =0.52631517256938.1 与f (x)=在各相应节点上数值相比，上述结果约有五位有效数字，精度较高。h越小,精度越高。附追赶法程序：%追赶法%定义三对角矩阵 A的各组成单元。方程为 Ax=d% a为主对角线下面的次对角线，b为主对角线，c为主对角线上面的次对角线，d为等式右端向量。% A=4100%1410%0141%0014function x=zhuiganfa(a,b,c,d)a=0 1 1 1;c=1 1 1;b=4 4 4 4;d

16、=3.75489429-0.666666667,3.53349105,3.33676905,3.16081554-0.5;n=length(b);u0=0;y0=0;a(1)=0;L(1)=b(1)-a(1)*u0;y(1)=(d(1)-y0*a(1)/L(1);u(1)=c(1)/L(1);for i=2:(n-1)L(i)=b(i)-a(i)*u(i-1);y(i)=(d(i)-y(i-1)*a(i)/L(i);u(i)=c(i)/L(i);endL(n)=b(n)-a(n)*u(n-1);y(n)=(d(n)-y(n-1)*a(n)/L(n);x(n)=y(n);for i=(n-1):-1:1x(i)=y(i)-u(i)*x(i+1);end运行结果：ans =0.624998286114830.588234478540670.555554849722490.52631517256938>> x=1.5:0.1:2.0;0.588234478540670.55555484972249y=log(x);dy=1/1.50.624998286114830.52631517256938 1/2.0;plot(x,y,'b-'