數(shù)值微分的計(jì)算方法

數(shù)值微分的計(jì)算方法內(nèi)容摘要求解數(shù)值微分問(wèn)題，就是通過(guò)測(cè)量函數(shù)在一些離散點(diǎn)上的值，求得函數(shù)的近似導(dǎo)數(shù)。本文就所學(xué)知識(shí)，歸納性地介紹了幾種常用的數(shù)值微分計(jì)算方法。并舉例說(shuō)明計(jì)算，實(shí)驗(yàn)結(jié)果表明了方法的有效性。關(guān)鍵詞數(shù)值微分Taylor展開(kāi)式Lagrange插值三對(duì)角矩陣引言：數(shù)值微分即根據(jù)函數(shù)在一些離散點(diǎn)的函數(shù)值，推算它在某點(diǎn)的導(dǎo)數(shù)或高階導(dǎo)數(shù)的近似值的方法。常見(jiàn)的可以用一個(gè)能夠近似代替該函數(shù)的較簡(jiǎn)單的可微函數(shù)(如多項(xiàng)式或樣條函數(shù)等)的相應(yīng)導(dǎo)數(shù)作為能求導(dǎo)數(shù)的近似值，由此也可導(dǎo)出多點(diǎn)數(shù)值微分計(jì)算公式。當(dāng)函數(shù)可微性不太好時(shí)，利用樣條插值進(jìn)行數(shù)值微分要比多項(xiàng)式插值更適宜。Taylor展開(kāi)式方法理論基礎(chǔ)：Taylor展開(kāi)式f(x)=f(x)+(x-x)f(x)+『2f〃(x)+...+￡，f(n)(x)+...TOC\o"1-5"\h\z0 0 0 2!o n! o我們借助Taylor展開(kāi)式，可以構(gòu)造函數(shù)fx在點(diǎn)x0的一階導(dǎo)數(shù)和二階導(dǎo)數(shù)的數(shù)值微分公式。取步長(zhǎng)h>0則f(x+h)=f(x)+hf'(x)+與f”&)&e(x,x+h) ⑴o o o2 1ioo所以f(x)=f(xo+h-f(xo)-hf惦)&e(x,x+h) (2)o h 2 1 1oo同理f(x—h)=f(x)—hf'(x)+gf”&)&e(x—h,x) (3)o oo2 2 2oof(x)—f(x—h)hf(x)= o~-~o +=f弋)&e(x—h,x) (4)o h 2 2 2oo式(2)和式(4)是計(jì)算fxo的數(shù)值微分公式，其截?cái)嗾`差為。h，為提高精度，將Taylor展開(kāi)式多寫(xiě)幾項(xiàng)f(x+h)=f(x)+hf，(x)+%f(x)+gf(x)+二f⑷(&)&e(x,x+h)o oo2o6o24iiooTOC\o"1-5"\h\zh2 h3 h4f(x—h)=f(x)—hf，(x)+=f”(x)—f”(x)+f(4)(&)&e(x—h,x)o oo2o6o24 2 2oo兩式相減得f(x+h)—f(x—h)h2 (5)f(x)= o o f'”(x)+O(h4) (5)o 2h 6o上式為計(jì)算f'(x°)的微分公式，其截?cái)嗾`差為O(h2),比式(2)和(4)精度高。兩式相加，如果f(4)(x)eC[xo—h,xo+h]，則有

f(x0)=f(X05)_2f(")+f(%-h)_豈⑷&)式(6)是計(jì)算f"(X0)的數(shù)值微分公式,其截?cái)嗾`差為O(h2)。例1.設(shè)函數(shù)fx=m%,x=2,h—o.i,試用數(shù)值微分公式計(jì)算f'2的值。0解由式⑵、式⑷和式⑸分別計(jì)算結(jié)果為f(2)a光嚴(yán)a048790.2ln2.1-ln1.9

f'(2)a a0.50040.2與真值f'2=05相比，式(5)計(jì)算的結(jié)果精度較高。數(shù)值微分的Lagrange插值方法設(shè)函數(shù)y=fX具有〃+1個(gè)實(shí)驗(yàn)數(shù)據(jù)：x「fXj i=0,1,2,...,〃，我們希望估計(jì)f'X的值，特別x=x,時(shí)，估計(jì)f，Xj的值?；诓逯捣椒ǖ臄?shù)值微分做法是，由已知X,fX j=0,1,???,n建立Lagrange插值多項(xiàng)式或Newton插值多項(xiàng)式(這里以Lagrange插值方法為例)，即fX=LX+RX于是f(x)=L(X)+R(x)當(dāng)X=X時(shí)，有jf(X.)=L(X.)+R〃(X.)(i=0,1,2,...,n)其中f(n+1)(&)R(x)=f——?、(x)

nj (n+1)! n+「/略去誤差項(xiàng)有f'(X.)aLn(X,)實(shí)際運(yùn)用中，等距節(jié)點(diǎn)更為常見(jiàn)。設(shè)h=~—,x.=a+ihj=0,1,2,-?，n, =th，nj于是有

t入 t(t—1)人 t(t—1)...(t—n+1)f(x)=y+ Ay+———A2y+...+ M+R(x)n!01! 0 2! 0 n! 0n!所以dt t(t—1)...(t—n+1)f‘(x)= {y+ Ay+...+ Any}idt01! 0 n! 0t=ihnf(n+1)(&)rf.、+——n n(i—/)(n+1)!j=0j歸多點(diǎn)數(shù)值微分公式由于高階插值的不穩(wěn)定性，實(shí)際應(yīng)用時(shí)多采用n=1,2,4的兩點(diǎn)、三點(diǎn)和五點(diǎn)等多點(diǎn)插值型求導(dǎo)公式。(1)兩點(diǎn)公式(n=1)f'(f'(氣)=1(七—匕)0h10f'(x)=1(y-y)1h1 0hR(x)=--f”(&)

0 2hR(x)=；f'電)(2)三點(diǎn)公式(n=2)f‘(x)=&-3y(2)三點(diǎn)公式(n=2)f‘(x)=&-3y+4y-y)

0 2h 0 1 2f'(x)=—(-y+y)2h0 2f'(x)=￡(y-4y+3y)2h0 1 2h2R(x)=3f⑶(&)03h2R(x)=-—f⑶(&)1 6h2R(x0)=—f⑶(&)(8)(3)五點(diǎn)公式(n=4)f'(x)=*-25y+48y-36y+16y-3y)0 12h 012341f'(x)=2-3y-10y+18y-6y+y)12h 0 1 2 3 4一、1 /cc 、<f'(x)=奇(y-8y+8y-y)12h0 1 3 4f(x3)=&(一*+6y1-18’2+10、3+3y「f'(x)=土(3y-16y+36y-48y+25y)L4 12h0 1 2 3 4h4R(x)=￡f⑸&)

0 5h4R(x)=——f(5)&)1 20h4R(x2)=-萬(wàn)f⑸&)h4R(x3)=-—f⑸&)R(x)=當(dāng)f(5)&)05(9)例2.設(shè)fx=lnx,取h=0.05,分別用三點(diǎn)公式和五點(diǎn)公式計(jì)算廣2的近似值。解：由式(8)有f'⑵=士(-3f⑵+4f"A一f(240))=0499802861f'⑵=士(一f(1-95)+f(2.05))=0.500104205f'⑵=d05(f(1-90)-4f(1-95)+3f(2))=0.499779376由式(9)有f'(2)二12x0.05(f'(2)二12x0.05(f(1.90)-8f(1.95)+6f(2.05)-f(2.10))=0.499999843數(shù)值微分的隱式格式前述的Taylor展式法和數(shù)值微分格式均稱(chēng)為顯式格式，即直接由已知的fXj(i=0,1,2,...,n)，經(jīng)過(guò)適當(dāng)?shù)乃阈g(shù)四則運(yùn)算，立即可得f'Xj的近似值。顯式格式優(yōu)點(diǎn)是計(jì)算方便，工作量小，缺點(diǎn)是數(shù)值不穩(wěn)定。為克服后一缺點(diǎn)，隱式格式常常具有數(shù)值穩(wěn)定性。數(shù)值微分的隱式格式建立方法常用通過(guò)Taylor展開(kāi)式方法或數(shù)值積分方法等不同途徑。首先我們用Taylor展開(kāi)式方法來(lái)推導(dǎo)數(shù)值微分的隱式格式。由式⑸和式⑹，我們用Xk代替x0得:TOC\o"1-5"\h\z° 1 T h2. ,f(x)=—[f(X+h)-f(x-h)]—乙f"'(x)+O(h4)

k 2h k k 6k1 _ h2f"(x)=—[f(X+h)-2f(X)+f(X-h)]-—f(4)&)

kh2k kk12所以也有1 一 h2f'"(X)=-[f(X+h)-2f(X)+f(X-h)]--f(5)&)

kh2k kk12將最后fXk表達(dá)式代入fXk表達(dá)式可得f'(X)=+[f(X+h)-f(X-h)]-![f，(X+h)-2f'(X)+f'(X-h)]+O(h4)

k2hk k6k kk略去誤差項(xiàng)Oh4，并且mk表示f'Xk的近似值，則有/ 3「m +4m+m=三[f(x )—f(x )] (k=1,2,...,n)同樣，用數(shù)值積分方法也可推導(dǎo)數(shù)值微分的隱式格式。L+1f(X)dx=JXk+1f'(X)dx (k=1,2,...,n)kXk-1 Xk-1根據(jù)simpson公式』bf(x)dx=b—a ab—a a+b^—[f(a)+4f(―^)+f(b)]有

6 2~ ~ 、2h「 / …… 、f(x)-f(x)=—[m+4m+m],(k=1,2,...,n)3=時(shí)f(Xk+1)-f(Xk-1)](k=1,2,...,n)上式是關(guān)于n+1個(gè)未知量m0,m^—m的n—1個(gè)方程，如果已知m0=八 。 ， 3…? 、―fX0,mn=f'Xn，記dk=-^[f(xk+1)一f牌「]，故有方程組:

「41 「41 ]md-m110141md2214 1md3=3???.." " "--?,,1 41mdn—2n-214md-m,——1u n-1」」n-1 n堂(10)這就是求m氣，…,七i的線性方程組。由于系數(shù)矩陣是嚴(yán)格對(duì)角占優(yōu)的三對(duì)角矩陣，因此非奇異，解存在且唯一，可由追趕法求解，且數(shù)值穩(wěn)定。例3.設(shè)函數(shù)fx=lnx在節(jié)點(diǎn)x=1.5,1.6,1.7,1.8,1.9,2.0上的函數(shù)值，及/'1-5,f'2.0如表1所示。試用數(shù)值微分隱式格式式(10)和式(11)求出相應(yīng)節(jié)點(diǎn)上一階1 一m141m2=141m314m—11- 4J解：由式(10)有導(dǎo)數(shù)的近似值。解：由式(10)有表1例3節(jié)點(diǎn)數(shù)據(jù)xf(x)f'(x)1.50.4054651080.6666666671.60.4700036291.70.5306282511.80.5877866641.90.6418538862.00.6931471820.5000000003.75489429-0.6666666673.533491053.336769053.16081554-0.500000000解方程組得:^1=0.62499828611483^2=0.58823447854067m3=0.55555484972249氣=0.52631517256938與f'(x)=1在各相應(yīng)節(jié)點(diǎn)上數(shù)值相比，上述結(jié)果約有五位有效數(shù)字，精度較高。h越小，x精度越高。附追趕法程序：%追趕法%定義三對(duì)角矩陣A的各組成單元。方程為Ax=d%a為主對(duì)角線下面的次對(duì)角線，b為主對(duì)角線，c為主對(duì)角線上面的次對(duì)角線，d為等式右端向量。%A=[4100% 1410% 0141% 0014]functionx=zhuiganfa(a,b,c,d)a=[0111];c=[111];b=[4444];d=[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);fori=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);fori=(n-1):-1:1x(i)=y(i)-u(i)*x(i+1);end運(yùn)彳丁結(jié)果：ans=0.62499828611483 0.58823447854067 0.55555484972249 0.52631517256938>>x=1.5:0.1:2.0;y=iog(x)；dy=[1/1.50.62499828611483 0.58823447854067 0.555554849722490.526315172569381/2.0];plot(x,y,'b-')holdonplot(x,1./x,'r-o')holdonplot(x,