【正文】
E(1)=0。E(2)=abs(D(2)D(1))。 %計算導(dǎo)數(shù)近似值間隔序列的前兩個值 R(1)=1。R(2)=2*E(2)/(abs(D(1))+abs(D(2)))。 %輸出滿足精度要求的數(shù)值導(dǎo)數(shù)近似值 for i=3:max %計算導(dǎo)數(shù)近似值序列 h=h/10。 D(i)=(feval(f,x+2*h)+8*feval(f,x+h)8*feval(f,xh)+feval(f,x2*h))/(12*h)。 E(i)=abs(D(i)D(i1))。 R(i)=2*E(i)/(abs(D(i))+abs(D(i1)))。 if ((E(i)E(i1))amp。(R(i)=toler)) break。 end end n=i。 df=D(n)。 %輸出滿足精度要求的數(shù)值導(dǎo)數(shù)近似值 err=E(n)。 relerr=R(n)。 end 理查森外推 法 計算 導(dǎo)數(shù) 近似值的子函數(shù) function df=Richarson(f,x,toler,delta,max,h0) %輸入 輸出參數(shù) 說明 : %df 輸出 的 近似 導(dǎo)數(shù)值 %f、 x 相應(yīng)函數(shù)及其點 delta 相對誤差容 忍上限, toler 誤差容忍上限, max 最大計算次數(shù), h0 初始 步長 h=h0。 D(1,1)=(feval(f,x+h)feval(f,xh))/(2*h)。 err=1。 《數(shù)值方法》實驗報告 11 relerr=1。 j=1。 while((errdelta)amp。(relerrtoler)amp。(j=max))%進行 是否終止判斷 h=h/2。 D(j+1,1)=(feval(f,x+h)feval(f,xh))/(2*h)。 for k=1:j %理查森 外推 公式 計 算高階 精度 的導(dǎo)數(shù)計算公式 D(j+1,k+1)=D(j+1,k)+(D(j+1,k)D(j+1,k))/(4^(k)1)。 end err=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。 end n=j df=D(n,n)。%輸出 得到滿足得到 精度 要求的近似導(dǎo)數(shù)值 err relerr End 1 function df=difnew(X,Y) N=length(X)1。 A=Y。 for M=1:N+1 %對 X 序列重新排列 得到 新的 插值點 序列 if M==1 t=X。 else if M==N+1 for i=2:N+1 t(i)=X(i1)。 end else t(1)=X(M)。 for i=2:M t(i)=X(i1)。 end for i=M+1:N+1 t(i)=X(i)。 end end end for i=1:N+1 %依據(jù)重新 排列的插值點序列,計算差商 表 以 得到 牛頓 插 《數(shù)值方法》實驗報告 12 值多項式 for i=j+1:N+1 A(i)=(A(i)A(i1))/(t(i)t(ij+1))。 end p=1。 df(M)=A(2)。 for k=2:N %累加 求和求解 插值點 導(dǎo)數(shù)值 df(M)=df(M)+p*(t(1)t(k))*A(k+1)。 end end end end