數(shù)值微分計(jì)算方法實(shí)驗(yàn)

1、數(shù)值微分計(jì)算方法實(shí)驗(yàn)【摘要】數(shù)值微分(numerical differentiation根據(jù)函數(shù)在一些離散點(diǎn)的函數(shù)值，推算它在某點(diǎn)的導(dǎo)數(shù)或高階導(dǎo)數(shù)的近似值的方法。通常用差商代替微商，或者用一個(gè)能夠近似代替該函數(shù)的較簡(jiǎn)單的可微函數(shù)（如多項(xiàng)式或樣條函數(shù)等）的相應(yīng)導(dǎo)數(shù)作為能求導(dǎo)數(shù)的近似值。例如一些常用的數(shù)值微分公式(如兩點(diǎn)公式、三點(diǎn)公式等就是在等距步長(zhǎng)情形下用插值多項(xiàng)式的導(dǎo)數(shù)作為近似值的。此外，還可以采用待定系數(shù)法建立各階導(dǎo)數(shù)的數(shù)值微分公式，并且用外推技術(shù)來(lái)提高所求近似值的精確度。當(dāng)函數(shù)可微性不太好時(shí)，利用樣條插值進(jìn)行數(shù)值微分要比多項(xiàng)式插值更適宜。如果離散點(diǎn)上的數(shù)據(jù)有不容忽視的隨機(jī)誤差，應(yīng)該用曲線擬

2、合代替函數(shù)插值，然后用擬合曲線的導(dǎo)數(shù)作為所求導(dǎo)數(shù)的近似值，這種做法可以起到減少隨機(jī)誤差的作用。數(shù)值微分公式還是微分方程數(shù)值解法的重要依據(jù)。關(guān) 鍵 詞 中心差分公式；理查森外推法；牛頓多項(xiàng)式微分；1、實(shí)驗(yàn)?zāi)康?.通過(guò)本次實(shí)驗(yàn)熟悉并掌握各種數(shù)值微分算法。2.掌握如何通過(guò)程序設(shè)計(jì)實(shí)現(xiàn)數(shù)值微分算法，從而更好地解決實(shí)際中的問(wèn)題。2、實(shí)驗(yàn)原理1. 精度為的中心差分公式：設(shè)，且，則 而且存在數(shù)，滿足其中項(xiàng)稱為截?cái)嗾`差。2. 精度為的中心差分公式：設(shè)，且，則 而且存在數(shù)，滿足其中項(xiàng)稱為截?cái)嗾`差。3.理查森外推法：利用低階公式推出高階求解數(shù)值微分的公式，定理如下：設(shè)的兩個(gè)精度為的近似值分別為和，而且它們滿足和這

3、樣可得到改進(jìn)的近似值表達(dá)式4. 牛頓多項(xiàng)式微分：利用的個(gè)插值點(diǎn)可得近似的次牛頓多項(xiàng)式，其可表示為則的導(dǎo)數(shù)可以表示為由此可知在處的導(dǎo)數(shù)值3、實(shí)驗(yàn)內(nèi)容1.P260 1用程序6.1求解下列函數(shù)在x處的倒數(shù)近似值，精度為小數(shù)點(diǎn)后13位。注：有必要改寫(xiě)程序中的max1的值和h的初始值。P260 2修改程序6.1，實(shí)現(xiàn)精度為O（h4）的中心差分公式（10）。用這個(gè)程序求解上題中給出的函數(shù)導(dǎo)數(shù)的近似值。精度為小數(shù)點(diǎn)后13位。P260 3使用程序6.2求解第1題中函數(shù)導(dǎo)數(shù)的近似值。精度為小數(shù)點(diǎn)后13位。注：有必要改變err，relerr和h的初始值。2.P270 1 修改程序6.3，使得可用它計(jì)算。四、實(shí)驗(yàn)結(jié)

4、果及分析（一）實(shí)驗(yàn)描述1.P260 1已知下列函數(shù)，用數(shù)值微分方法求解在處的導(dǎo)數(shù)近似值，精度為小數(shù)點(diǎn)后面的13位。（1）P260 1使用精度的中心差分公式求解（2）P260 2使用精度的中心差分公式求解（3）P260 3使用理查森外推法求解2.P270 1 計(jì)算階牛頓插值多項(xiàng)式的在插值點(diǎn)的導(dǎo)數(shù)值（二）實(shí)驗(yàn)算法1.P260 1計(jì)算精度為導(dǎo)數(shù)近似值的算法.依據(jù)精度為的中心差分公式，利用步長(zhǎng)序列。Step1：輸入函數(shù)，要計(jì)算的導(dǎo)數(shù)的橫坐標(biāo)及相對(duì)誤差容限toler，步長(zhǎng)序列長(zhǎng)度；Step2：設(shè)置初始參量，初始步長(zhǎng)，在處的近似導(dǎo)數(shù)值令，計(jì)算并令近似導(dǎo)數(shù)值序列前兩個(gè)初值，相對(duì)誤差序列前兩個(gè)初值。Step3

5、：for i=2 upto max1） 計(jì)算步長(zhǎng)時(shí)的近似導(dǎo)數(shù)值，；2） 近似導(dǎo)數(shù)值序列間隔，相對(duì)誤差序列。3） 若或，則停止計(jì)算循環(huán)；Step4：輸出即可得到符合精度要求的近似導(dǎo)數(shù)值，記錄。Step5：將儲(chǔ)存到向量中列表輸出。2.P260 2計(jì)算精度為導(dǎo)數(shù)近似值的算法.依據(jù)精度為的中心差分公式，利用步長(zhǎng)序列。Step1：輸入函數(shù)，要計(jì)算的導(dǎo)數(shù)的橫坐標(biāo)及相對(duì)誤差容限toler，步長(zhǎng)序列長(zhǎng)度；Step2：設(shè)置初始參量，初始步長(zhǎng)，在處的近似導(dǎo)數(shù)值令，計(jì)算并令近似導(dǎo)數(shù)值序列前兩個(gè)初值，相對(duì)誤差序列前兩個(gè)初值。Step3：for i=1 upto max1） 計(jì)算步長(zhǎng)時(shí)的近似導(dǎo)數(shù)值，；2） 近似導(dǎo)數(shù)值序

6、列間隔，相對(duì)誤差序列。3） 若或，則停止計(jì)算循環(huán)；Step4：輸出即可得到符合精度要求的近似導(dǎo)數(shù)值，記錄。Step5：將儲(chǔ)存到向量中輸出。3.P260 3利用理查森外推法設(shè)計(jì)計(jì)算滿足精度要求的近似導(dǎo)數(shù)值的算法，Step1：輸入函數(shù)，自變量，相對(duì)誤差容限toler，誤差容限delta，初始步長(zhǎng)，計(jì)算行數(shù)上限；Step2：計(jì)算初始化量1） 初始近似導(dǎo)數(shù)值；2） 初始誤差，初始相對(duì)誤差；3） ；Step3：計(jì)算出符合精度要求的近似導(dǎo)數(shù)值：若且且，則進(jìn)行下列計(jì)算，否則進(jìn)入step41）,;For k=1 upto j2 計(jì)算誤差3計(jì)算相對(duì)誤差4,返回step3Step4：輸出 ，為符合要求的導(dǎo)數(shù)近似值

7、。4.P270 1 利用牛頓多項(xiàng)式微分求解各插值點(diǎn)的近似的導(dǎo)數(shù)值算法流程如下：Step1：輸入及其對(duì)應(yīng)的函數(shù)值，Step2：for k=0 upto N，計(jì)算處的近似導(dǎo)數(shù)值1），2令，for j=1 upto N+1for i=j+1 upto N+13）p=1，for n=2 upto NStep3：向量即處的近似導(dǎo)數(shù)值，（3）實(shí)驗(yàn)結(jié)果1.P260 （1）利用精度為的中心差分公式求解結(jié)果（2）利用精度為的中心差分公式求解結(jié)果（3）利用理查森外推法求解結(jié)果3種方法所求的結(jié)果如下，結(jié)果的精度為13位有效數(shù)字表1 三種方法求解各個(gè)函數(shù)在其相應(yīng)x處近似導(dǎo)數(shù)值對(duì)應(yīng)函數(shù)及其求解點(diǎn)（以字母序號(hào)代替）理查森

8、外推法（a）75.1735024617575.1734947073175b）1.2286029410281.2285974234041.228733630836（c）1.9515627659521.9515596090241.951752228326（d）2.9655290467632.9655148293492.966366834548（e）1.0152010645781.0150610589881.015008921184由上表可知，三種方法的求解結(jié)果基本相同，但是因?yàn)楹瘮?shù)及求解點(diǎn)的設(shè)置的初始步長(zhǎng)和終止條件的不同，三種方法的求解結(jié)果沒(méi)有做到高度一致，因此可知求解結(jié)

9、果與初始步長(zhǎng)及終止條件的設(shè)定有較大的關(guān)系。五、實(shí)驗(yàn)結(jié)論本次數(shù)值實(shí)驗(yàn)利用了精度為、的中心差分公式及理查森外推公式計(jì)算了各個(gè)函數(shù)在其對(duì)應(yīng)點(diǎn)的近似導(dǎo)數(shù)值，三種方法的求解結(jié)果基本相同，但是因?yàn)楹瘮?shù)及求解點(diǎn)的設(shè)置的初始步長(zhǎng)和終止條件的不同，三種方法的求解結(jié)果沒(méi)有做到高度一致，因此可知求解結(jié)果與初始步長(zhǎng)及終止條件的設(shè)定有較大的關(guān)系。另外還設(shè)計(jì)了計(jì)算利用牛頓插值多項(xiàng)式的任意插值點(diǎn)導(dǎo)數(shù)計(jì)算算法，并且用相應(yīng)程序?qū)崿F(xiàn)，其中將插值點(diǎn)重新排列，使得所求點(diǎn)為第一個(gè)點(diǎn)是一個(gè)較大的創(chuàng)新，這樣使得各個(gè)插值點(diǎn)的導(dǎo)數(shù)計(jì)算更為快捷方便。附件（代碼）：1.P260 精度為近似導(dǎo)數(shù)值的計(jì)算子程序function df=cetr2df(

10、f,x,toler,max,h0%輸入輸出參數(shù)說(shuō)明：%df 輸出的近似導(dǎo)數(shù)值%f、x相應(yīng)函數(shù)及其點(diǎn)，toler誤差容忍上限，max最大計(jì)算次數(shù)，h0初始步長(zhǎng)for i=1:2h=10(-(i-1*h0;D(i=(feval(f,x+h-feval(f,x-h/(2*h;%計(jì)算導(dǎo)數(shù)近似值序列的前兩個(gè)值endE(1=0;E(2=abs(D(2-D(1;%計(jì)算導(dǎo)數(shù)近似值間隔序列的前兩個(gè)值R(1=1;R(2=2*E(2/(abs(D(1+abs(D(2;%計(jì)算導(dǎo)數(shù)近似值相對(duì)誤差序列的前兩個(gè)值for i=3:max%計(jì)算導(dǎo)數(shù)近似值序列h=h/10;D(i=(feval(f,x+h-feval(f,x-h

11、/(2*h;E(i=abs(D(i-D(i-1;R(i=2*E(i/(abs(D(i+abs(D(i-1;if (E(i>E(i-1|(R(i<=toler%計(jì)算導(dǎo)數(shù)近似值序列的終止條件break;endendn=i;df=D(n;%輸出滿足精度要求的數(shù)值導(dǎo)數(shù)近似值end精度為近似導(dǎo)數(shù)值的計(jì)算子程序function df=cetr4df(f,x,toler,max,h0%輸入輸出參數(shù)說(shuō)明：%df 輸出的近似導(dǎo)數(shù)值%f、x相應(yīng)函數(shù)及其點(diǎn)，toler誤差容忍上限，max最大計(jì)算次數(shù)，h0初始步長(zhǎng)for i=1:2%計(jì)算導(dǎo)數(shù)近似值序列的前兩個(gè)值h=10(-(i-1*h0;D(i=(-fe

12、val(f,x+2*h+8*feval(f,x+h-8*feval(f,x-h+feval(f,x-2*h/(12*h;endE(1=0;E(2=abs(D(2-D(1; %計(jì)算導(dǎo)數(shù)近似值間隔序列的前兩個(gè)值R(1=1;R(2=2*E(2/(abs(D(1+abs(D(2; %輸出滿足精度要求的數(shù)值導(dǎo)數(shù)近似值for i=3:max %計(jì)算導(dǎo)數(shù)近似值序列h=h/10;D(i=(-feval(f,x+2*h+8*feval(f,x+h-8*feval(f,x-h+feval(f,x-2*h/(12*h;E(i=abs(D(i-D(i-1;R(i=2*E(i/(abs(D(i+abs(D(i-1;if

13、 (E(i>E(i-1&(R(i<=tolerbreak;endendn=i;df=D(n; %輸出滿足精度要求的數(shù)值導(dǎo)數(shù)近似值err=E(n;relerr=R(n;end理查森外推法計(jì)算導(dǎo)數(shù)近似值的子函數(shù)function df=Richarson(f,x,toler,delta,max,h0%輸入輸出參數(shù)說(shuō)明：%df 輸出的近似導(dǎo)數(shù)值%f、x相應(yīng)函數(shù)及其點(diǎn)delta相對(duì)誤差容忍上限，toler誤差容忍上限，max最大計(jì)算次數(shù)，h0初始步長(zhǎng)h=h0;D(1,1=(feval(f,x+h-feval(f,x-h/(2*h;err=1;relerr=1;j=1;while(er

14、r>delta&(relerr>toler&(j<=max%進(jìn)行是否終止判斷h=h/2;D(j+1,1=(feval(f,x+h-feval(f,x-h/(2*h;for k=1:j %理查森外推公式計(jì)算高階精度的導(dǎo)數(shù)計(jì)算公式D(j+1,k+1=D(j+1,k+(D(j+1,k-D(j+1,k/(4(k-1;enderr=abs(D(j+1,j+1-D(j,j;relerr=2*err/(abs(D(j+1,j+1+abs(D(j,j;j=j+1;endn=jdf=D(n,n;%輸出得到滿足得到精度要求的近似導(dǎo)數(shù)值errrelerrEnd2.P270 1 function df=difnew(X,YN=length(X-1;A=Y;for M=1:N+1 %對(duì)X序列重新排列得到新的插值點(diǎn)序列if M=1t=X;else if M=N+1for i=2:N+1t(i=X(i-1;en