一阶导数的数值计算及其MATLAB程序

第三篇 第八章 数值微分第八章 数值微分8.2 一阶导数的数值计算及其MATLAB程序8.2.1 差商求导及其MATLAB程序例 8.2.1 设.（1）分别利用前差公式和后差公式计算的近似值和误差，取4位小数点计算，其中步长分别取，80,.（2）将（1）中计算的的近似值分别与精确值比较.解 （1）编写计算的一阶导数计算的近似值和误差估计的MATLAB程序，并输入 x=0.79;h=0.1,0.01,0.001,0.0001;M=80;x1=x+h;x2=x-h; y=sin(5.*x.2-21);y1=sin(5.*x1.2-21); y2=sin(5.*x2.2-21); yq=(y1-y)./h, yh=(y-y2)./h,wu=abs(h.*M/2), syms x,f=sin(5.*x.2-21); yx=diff(f,x)运行后屏幕显示利用前差公式和后差公式计算的近似值yq，yh和误差估计wu，取4位小数点计算，其中步长分别取，M=80，导函数yxyq =1.46596380397978 4.22848550173043 4.44250759584697 4.46320955293622yh =5.96885352366536 4.68672022108227 4.48833808130555 4.46779260847907wu =4.00000000000000 0.40000000000000 0.04000000000000 0.00400000000000yx =10*cos(5*x2-21)*x（2）计算的值.输入程序 x=0.79; yx =10*cos(5*x2-21)*x,wuq=abs(yq-yx), wuh=abs(yh-yx)运行后屏幕显示利用前差公式和后差公式计算的近似值与精确值的绝对误差wuq，wuh和的精确值yx如下yx = 4.46550187104484wuq =2.99953806706506 0.23701636931441 0.02299427519787 0.00229231810861wuh =1.50335165262053 0.22121835003744 0.02283621026072 0.002290737434248.2.2 中心差商公式求导及其MATLAB程序利用精度为的三点公式计算的近似值和误差估计的MATLAB主程序function n,xi,yx,wuc=sandian(h,xi,fi,M)n=length(fi); yx=zeros(1,n); wuc=zeros(1,n); x1= xi(1); x2= xi(2); x3= xi(3);y1=fi(1); y2=fi(2); y3=fi(3); xn= xi(n); xn1= xi(n-1); xn2= xi(n-2); yn=fi(n); yn1=fi(n-1); yn2=fi(n-2);for k=2:n-1yx(1)=(-3*y1+4*y2-y3)/(2*h); yx(n)=(yn2-4*yn1+3*yn)/(2*h);yx(2)=( fi(3)- fi(1)/(2*h); yx(k)=( fi(k+1)- fi(k-1)./(2*h); wuc(1)=abs(h.2.*M./3); wuc(n)=abs(h.2.*M./3);wuc(2:n-1)=abs(-h.2.*M./6);end利用精度为的三点公式计算的近似值和误差估计的MATLAB主程序function x,yxj, wuc=sandian3(h,xi,fi,M)yxj=zeros(1,3); wuc=zeros(1,3); x1= xi(1);x2= xi(2); x3= xi(3); y1=fi(1); y2=fi(2); y3=fi(3);for t=1:3s(t)=(2*t-5)*y1-4*(t-2)*y2+(2*t-3)*y3)/(2*h); x=xi; y=s(t); yxj(t)=y;if t=2wuc(2)=abs(-h.2*M/6);elsewuc(1:2:3)=abs(h.2*M/3);endend例 8.2.3 设已给出的数据表85：表85x1.000 0 1.100 0 1.200 0 1.300 0 1.400 0 1.500 0 1.600 0 f(x)0.250 0 0.226 8 0.206 6 0.189 0 0.173 6 0.160 0 0.147 9 M= 0.750 2,试用三点公式计算下列问题：(1)当h=0.1时，在x=1.000 0，1.100 0，1.200 0，1.300 0，1.400 0，1.500 0，1.600 0处的一阶导数的近似值，并估计误差；(2)当h=0.2时，在x=1.000 0，1.200 0， 1.400 0，1.600 0处的一阶导数的近似值，并估计误差；(3)当h=0.3时，在x=1.000 0，1.300 0 ，1.600 0处的一阶导数的近似值，并估计误差；(4) 表85中的数据是函数在相应点的数值，将(1)至(3)计算的一阶导数的近似值与的一阶导数值比较，并求出它们的绝对误差.解 （1）保存M文件sandian.m，sandian3.m；（2）在MATLAB工作窗口输入如下程序 syms x,y=1/(1+x)2); yx=diff(y,x,1),yx3=diff(y,x,3),运行后将屏幕显示的结果为yx = yx3 =-2/(1+x)3 -24/(1+x)5（3）在MATLAB工作窗口输入如下程序h=0.1; xi=1.0000:h:1.6000;fi=0.2500 0.2268 0.2066 0.1890 0.1736 0.1600 0.1479;x=1:0.001:1.6; yx3 =-24./(1+x).5; M= max(abs(yx3);n1,x1,yx1,wuc1=sandian(h,xi,fi,M)yxj1=-2./(1+xi).3,wuyxj1=abs(yxj1- yx1)h=0.2; xi=1.0000:h:1.6000; fi=0.2500 0.2066 0.1736 0.1479;x=1:0.001:1.6; yx3 =-24./(1+x).5; M= max(abs(yx3);n2,x2,yx2,wuc2=sandian(h,xi,fi,M)yxj2=-2./(1+xi).3,wuyxj2=abs(yxj2- yx2)h=0.3; xi=1.0000:h:1.6000; fi=0.2500 0.1890 0.1479;x=1:0.001:1.6; yx3 =-24./(1+x).5; M=max(abs(yx3);x3,yx3, wuc3=sandian3(h,xi,fi,M)yxj3=-2./(1+xi).3,wuyxj3=abs(yxj3- yx3)或 h1=0.1,x=1.0000,1.1000,1.2000,1.3000,1.4000,1.5000,1.6000;f=0.2500,0.2268,0.2066,0.1890,0.1736,0.1600,0.1479;xi=x(1:3);f11=f(1:3); M= 0.7502;x11,yxj11,wuc11=sandian3(h1,xi,f11,M), xi= x(4:6);f12=f(4:6); x12,yxj12,wuc12=sandian3(h1,xi,f12,M), xi=x(5:7);f13=f(5:7); x13,yxj13,wuc13=sandian3(h1,xi,f13,M), h2=0.2, xi= x(1:2:5);f21= f(1:2:5);x21,yxj21,wuc21=sandian3(h2,xi,f21,M),xi= x(2:2:6);f22=f(2:2:6); x22,yxj22,wuc22=sandian3(h2,xi,f22,M),xi= x(3:2:7);f23=f(3:2:7); x23,yxj23,wuc23=sandian3(h2,xi,f23,M), h3=0.3, xi= x(1:3:7);f31= f(1:3:7); x31,yxj31,wuc31=sandian3(h3,xi,f31,M),将运行的结果（略）.8.2.3 理查森外推法求导及其MATLAB程序（一）一般形式的理查森外推法及其MATLAB程序利用理查森外推法计算的近似值和误差估计的MATLAB程序function Dy,dy,jdw,n=diffext1(fun,x0,jdwc,max1)h=1;j=1; n=1;jdW=1;xdW=1; x1=x0+h;x2=x0-h;Dy(1,1)=(feval(fun,x1)- feval(fun,x2)/(2*h); while(jdWjdwc)&(j x0=0.79;jdwc=0.0000001,max1=100;Dy,dy,jdw,n=diffext1(fun,x0,jdwc,max1),wu=4.46550187104484-dy运行后屏幕显示的近似值dy，dy与精确值的绝对误差wu，导数近似值的迭代矩阵Dy，jdW=|Dy(n,n)-Dy(n,n-1)|的值,最佳近似值dy的坐标n如下Dy = Columns 1 through 4 0.95036708207779 0 0 00.87447447334140 0.84917693709594 0 01.04543344913993 1.10241977440611 1.11930263022679 03.33291848491782 4.09541349684379 4.29494641167297 4.345353455822914.16327772677060 4.44006414072153 4.46304085031338 4.465709016006084.38872901334021 4.46387944219674 4.46546712896176 4.465505641321264.44623223629059 4.46539997727405 4.46550134627921 4.46550188941123 Columns 5 through 7 0 0 0 0 0 0 0 0 0 0 0 0 4.46618099859504 0 0 4.46550484377347 4.46550418282057 0 4.46550187469786 4.46550187179553 4.46550187123118dy = jdw = n= wu = 4.46550187123118 5.643530087695581e-010 7 -1.863398324530863e-010（二）精度为的中心差商公式和误差估计及其MATLAB程序用精度为的中心差商公式计算的近似值和误差估计的MATLAB程序function x0,yx,wuc=zxcs4(h,x0,fi,M)xi=x0-2*h,x0-h,x0,x0+h,x0+2*h;x1= xi(1); x2= xi(2); x3= xi(3); x4= xi(4);x5= xi(5);y1=fi(1); y2=fi(2); y3=fi(3);y4=fi(4); y5=fi(5);yx=(8*y4-8*y2-y5+y1)/(12*h); wuc=abs(h.4*M/30);用精度为的中心差商公式计算的近似值和误差估计的MATLAB程序function x,y,yx,wuc=zxcs5(fun,h,x0, M)x=zeros(1,5);y= zeros(1,5);for k=1:5x(k)=x0+(k-3).*h;y(k)= feval(fun, x(k);endx;y;yx=(8*y(4)-8*y(2)-y(5)+y(1)/(12*h); wuc=abs(h.4*M/30);例 8.2.6 设.（1）分别根据(8.7) 式和(8.24)式计算的近似值，并估计误差，取小数点后14位计算，其中步长分别取 ，M=872,.（2）将（1）中计算的的近似值分别与精确值比较.解 （1）计算的一阶导数的近似值和误差估计.方法1 输入程序x0=0.79;h1=0.1,M=872;x,y,yx1,wuc1=zxcs5(fun,h1,x0,M)h2=0.01, x,y,yx2,wuc2=zxcs5(fun,h2,x0, M), h3=0.001, x,y,yx3,wuc3=zxcs5(fun,h3,x0, M), h4=0.0001, x,y,yx4,wuc4=zxcs5(fun,h4,x0, M),运行后屏幕显示x=(-2h, -h, ,+h, +2h), 在x处函数值向量y，步长分别取时，根据 (8.24)式计算函数在处的导数的近似值yx1, yx2, yx3, yx4和误差估计wuc1，wuc2，wuc3，wuc4如下x = Columns 1 through 3 0.59000000000000 0.69000000000000 0.79000000000000 Columns 4 through 5 0.89000000000000 0.99000000000000y = Columns 1 through 3 -0.39855804055354 0.22803197210174 0.82491732446828 Columns 4 through Column 5 0.97151370486625 0.38160930304430yx1 = wuc1 = 4.30640543209856 0.00290666666667yx2 = wuc2 = 4.46548476000396 2.906666666666667e-007yx3 = wuc3 = 4.46550186933213 2.906666666666667e-011yx4 = wuc4 = 4.46550187103480 2.906666666666667e-015方法2编写计算的一阶导数的近似值和误差估计的MATLAB程序，并输入此程序 x=0.79;h=0.1,0.01,0.001,0.0001; M=872;x1=x+h;x2=x-h;x3=x+2*h;x4=x-2*h; y1=sin(5.*x1.2-21);y2=sin(5.*x2.2-21);y3=sin(5.*x3.2-21); y4=sin(5.*x4.2-21); yc2=(y1-y2)./(2*h),wuc2=abs(h.2*M/6),yc4=(8*y1-8*y2-y3+y4)./ (12*h), wuc4=abs(h.4*M/30),syms x,f=sin(5.*x.2-21); dy=diff(f,x)运行后屏幕显示步长分别取，M= 872，分别根据(8.7) 式和(8.24)式计算的近似值yc2，yc4及其误差估计值wuc2，wuc4，导函数dy如下yc2 =3.71740866382257 4.45760286140635 4.46542283857626 4.46550108070765wuc2 =1.45333333333333 0.01453333333333 0.00014533333333 0.00000145333333yc4 =4.30640543209856 4.46548476000396 4.46550186933213 4.46550187103480wuc4 =0.00290666666667 0.00000029066667 0.00000000002907 0.00000000000000dy =10*cos(5*x2-21)*x（2）计算的值.输入程序 x=0.79; dy =10*cos(5*x2-21)*x, wu2=abs(yc2-dy), wu4=abs(yc4-dy)运行后屏幕显示的近似值与精确值的绝对误差wuc2，wuc4和的精确值dy如下dy =4.46550187104484wuc2=0.74809320722227 0.00789900963848 0.00007903246857 0.00000079033719wuc4=0.15909643894627 0.00001711104088 0.00000000171271 0.00000000001003（三）变步长的中心差商公式及其MATLAB程序用变步长的中心差商公式计算的近似值和误差估计的MATLAB程序function n,H,Dy,W=difflim(fun,x0,h0,wu,max1)Dy=zeros(1,max1); W=zeros(1, max1); H=zeros(1, max1);h=h0;H(1)=h; E(1)=0; x1=x0+h;x2=x0-h; h1=h/10;x3=x0+h1;x4=x0-h1;Dy(1)=(feval(fun,x1)- feval(fun,x2)./(2*h);Dy(2)=(feval(fun,x3)- feval(fun,x4)./(2*h1);W(1)=abs(Dy(2)- Dy(1); k=1; while( (W(k)wu)&(k x0=0.79;wu=0.00001;max1=100;h0=0.2;n,H,Dy,W= difflim(fun,x0,h0,wu,max1)jdwc=4.46550187104484- Dy运行后屏幕显示迭代次数n, 变步长数组H，用变步长的中心差商公式计算的近似值Dy，误差估计值W，jdwc，精确值与近似值Dy的差如下n = 2H = 0.02000000000000 0.00020000000000 0.00000020000000Dy = 4.43395716561354 4.46549870972618 4.46550187410688W = 2.48353880661895 0.03154154411264 0.00000316438070jdwc = 0.03154470543130 0.00000316131866 -0.000000003062048.2.4 牛顿（Newton）多项式求导及其MATLAB程序利用牛顿插值多项式求导的MATLAB主程序function df,A,P=diffnew(X,Y)n=length(X);A=Y;for j=2:n for i=n:-1:j A(i)=(A(i)- A(i-1)/(X(i)-X(i-j+1); end endx0=X(1);df=A(2);chsh=1;m=length(A)-1; for k=2:m chsh=chsh*(x0-X(k); df=df+chsh*(A(k+1);end P=poly2sym(A);例8.2.9 根据下表给定的一组数据（X ,Y）写出 (8.30)式的具体形式及其精度和名称，并用它计算的近似值，取4位小数计算.X3.135 2 3.335 2 3.535 2 3.735 2 3.935 2 4.135 2 4.335 2 4.535 2 Y0.126 6 -0.060 2 -0.603 2 -0.998 0 -0.119 4 0.995 3 -0.654 2 0.158 1 解 因为表中所给数据（X ,Y）是等差数列，公差为h=0.2, 即x0=3.135 2, x1= x0+h, x2= x0+2h, x3= x0+3h, x4= x0+4h, x5= x0+5h, x6= x0+6h, x7= x0+7h.输入程序X=3.1352,3.3352,3.5352,3.7352,3.9352,4.1352,4.3352,4.5352 ;Y=0.1266,-0.0602,-0.6032,-0.9980,-0.1194,0.9953,-0.6542 0.1581;df,A,P=diffnew(X,Y)运行后屏幕显示的近似值df和阶牛顿多项式P及其系数向量A如下df = -0.2428A =0.1266 -0.9340 -4.4525 10.5083 16.1667 -72.4818 64.7309 108.6155P=633/5000*x7-467/500*x6-1781/400*x5+1478916440133901/140737488355328*x4+4550512123488957/281474976710656*x3-27833/384*x2+4555032337958703/70368744177664*x+7643132912526197/70368744177664例8.2.10 设函数.利用(8.23)式的精度为的的阶前差公式、阶后差公式和阶中心差商公式计算的近似值，并与精确值比较，取14位小数计算，h=0.001.解 （1）输入计算程序 x0=3.1352;h=0.001; Xq=3.1352, x0+h, x0+2*h, x0+3*h, x0+4*h;Xh=3.1352, x0-h, x0-2*h, x0-3*h, x0-4*h; Xz=3.1352, x0+h, x0-h, x0+2*h,x0-2*h;Yq=sin(5*Xq.2-21); Yh=sin(5*Xh.2-21); Yz=sin(5*Xz.2-21);dfq,Aq,Pq=diffnew(Xq,Yq), dfh,Ah,Ph=diffnew(Xh,Yh), dfz,Az,Pz=diffnew(Xz,Yz)syms x, f=sin(5.*x.2-21) ; dy=diff(f,x)运行后屏幕显示阶前差公式Pq和它的系数向量Aq及其的近似值dfq；阶后差公式Ph和它的系数向量Ah及其的近似值dfh；阶中心差商公式Pz和它的系数向量Az及其的近似值dfz，导数的符号形式分别如下dfq =-31.09973888240658Aq =1.0e+003 * Columns 1 through 4 0.00012659805395 -0.03116184195521 -0.05190522279697 5.11405450158927 Column 5 5.04316638931140Pq=4561175588748639/36028797018963968*x4-8771278738603423/281474976710656*x3-7305010688969305/140737488355328*x2+2811481194788801/549755813888*x+5545020085856995/1099511627776dfh =-31.09973905786443Ah =1.0e+003 * Columns 1 through 4 0.00012659805395 -0.03102749818686 -0.08234817664867 5.02372948127046 Column 5 8.97444729902917Ph=4561175588748639/36028797018963968*x4-8733464329536931/281474976710656*x3-5794737776087263/70368744177664*x2+345228061216123/68719476736*x+5483510392960745/549755813888dfz =-31.09974389667719Az =1.0e+003 * Columns 1 through 4 0.00012659805395 -0.03116184195521 -0.06717188417315 5.08888712539359 Column 5 7.53074172148027Pz=4561175588748639/36028797018963968*x4-8771278738603423/281474976710656*x3-4726801133312021/70368744177664*x2+5595290566809831/1099511627776*x+2070034522136353/274877906944dy =10*cos(5*x2-21)*x（2）为了将的三种近似值与精确值比较，输入程序 x=3.1352; dy =10*cos(5*x2-21)*x, wuq=dfq- dy , wuh=dfh- dy, wuz=dfz- dy,运行后屏幕显示值及其三种近似值分别与精确值的差如下dy = wuq =-31.09974488361178 6.001205203887139e-006wuh = wuz =5.825747354748501e-006 8.869345909407912e-0078.3 高阶导数的数值计算及其MATLAB程序8.3.1 插值或拟合高阶数值导数及其MATLAB程序利用拉格朗日插值多项式构造阶导数的数值计算公式的MATLAB主程序function C,L,dyk,k=ndaolag(X,Y,n)m=length(X); n1=m; L=ones(m,m);for k=1:m V=1; for i=1:m if k=i V=conv(V,poly(X(i)/(X(k)-X(i); endendL1(k,:)=V; l(k,:)=poly2sym (V);endC=Y*L1;L=Y* l; syms x dykfor k=1:n k;dyk=diff(L,x,k)end例 8.3.2 给出节点数据， ，.（1）作五次拉格朗日插值多项式L和的1至5阶数值导数公式；（2）利用此公式求在处的1至5阶导数的近似值.解 （1）保存名为ndaolag.m的M文件.（2）在MATLAB工作窗口输入程序 X=-2.15 -1.00 0.01 1.02 2.03 3.25;Y=17.03 7.24 1.05 2.03 17.06 23.05; C,L,dyk,k= ndaolag(X,Y,5)运行后屏幕显示五次拉格朗日插值多项式L及其系数向量C，的1至5阶数值导数公式dyk（略）.（3）为了求在处的1至5阶导数的近似值，输入程序 x=-1.2345;dy1=-26370266994304203933/5764607523034234880+18224487282009991221/2882303761517117440*x2+59786195406624056511/230584300921369395200*x3-12501150855594615669/11529215046068469760*x4+2505830415074824099257/368934881474191032320*xdy2=18224487282009991221/1441151880758558720*x+179358586219872169533/230584300921369395200*x2-12501150855594615669/2882303761517117440*x3+2505830415074824099257/368934881474191032320dy3=18224487282009991221/1441151880758558720+179358586219872169533/115292150460684697600*x-37503452566783847007/2882303761517117440*x2dy4=179358586219872169533/115292150460684697600-37503452566783847007/1441151880758558720*xdy5 =-37503452566783847007/1441151880758558720运行后屏幕显示在处的1至5阶导数的近似值如下dy1 = dy2 = dy3 = dy4 = dy5 =-6.3294 0.5262 -8.1043 33.6814 -26.0232例 8.3.3 已知，.（1）作四次拉格朗日插值多项式L和的1至4阶数值导数公式;（2）利用上面的公式求在处的1至4阶导数的精确值，近似值及其绝对误差，取小数点后4位和后14位计算.解 保存名为ndaolag.m的M文件.在MATLAB工作窗口输入程序 X= pi/6,pi/4,pi/3,5*pi/12,pi/2; Y=0.5,0.7071,0.8660,0.9659,1; C,L,dyk,k=ndaolag(X,Y,4), x=pi/6: pi/12: pi/2; y=sin(x); C1,L1,dyk1,k1= ndaolag(x,y,4),for i=1:4i,syms x,dyi=diff(sin(x),x,i)end运行后屏幕显示四次拉格朗日插值多项式L及其系数向量C，的1至4阶数值导数dyk和符号数dyi公式及其导数阶数k和i.输入程序 x=2*pi/9;dY1=685769833743917463/703687441776640000+1261982467915759/10995116277760000*x-24314514941446017/35184372088832000*x2+6241572970666541/43980465111040000*x3dY2=1261982467915759/10995116277760000-24314514941446017/17592186044416000*x+18724718911999623/43980465111040000*x2dY3=-24314514941446017/17592186044416000+18724718911999623/21990232555520000*xdY4=18724718911999623/21990232555520000dy1=-6644262608911145859863501570909/79228162514264337593543950336*x3+77172015474009052195694711806203/316912650057057350374175801344*x2-34408318598267032320107245880583/158456325028528675187087900672*x+37371347371255248186645713466863/633825300114114700748351602688+1/2*2(1/2)*(-2496629188266579/17592186044416*x3+16667207822766783/35184372088832*x2-8898045875479213/17592186044416*x+6002949942623721/35184372088832)+1/2*3(1/2)*(7489887564799735/35184372088832*x3-11765087874894195/17592186044416*x2+91239728215363/137438953472*x-7257297691828675/35184372088832)dy2=-19932787826733437579590504712727/79228162514264337593543950336*x2+77172015474009052195694711806203/158456325028528675187087900672*x-34408318598267032320107245880583/158456325028528675187087900672+1/2*2(1/2)*(-7489887564799737/17592186044416*x2+16667207822766783/17592186044416*x-8898045875479213/17592186044416)+1/2*3(1/2)*(22469662694399205/35184372088832*x2-11765087874894195/8796093022208*x+91239728215363/137438953472)dy3=-19932787826733437579590504712727/39614081257132168796771975168*x+77172015474009052195694711806203/158456325028528675187087900672+1/2*2(1/2)*(-7489887564799737/8796093022208*x+16667207822766783/17592186044416)+1/2*3(1/2)*(22469662694399205/17592186044416*x-11765087874894195/8796093022208)dy4=-19932787826733437579590504712727/39614081257132168796771975168-7489887564799737/17592186044416*2(1/2)+22469662694399205/35184372088832*3(1/2)dyi1=cos(x),dyi2 =-sin(x),dyi3 =-cos(x),dyi4 =sin(x)WuY1= dyi1-dY1, WuY2= dyi2-dY2, WuY3= dyi3-dY3, WuY4= dyi4-dY4,wuy1= dyi1-dy1, wuy2= dyi2-dy2, wuy3= dyi3-dy3, wuy4= dyi4-dy4,运行后将输出的在处的1至4阶导数 (取小数点后4位和后14位) 的近似值及其精确值，绝对误差的计算结果（略）.8.3.2 高阶泰勒（Taylor）数值导数及其MATLAB程序例 8.3.5 设.（1）分别利用（8.31）和（8.30）式计算的近似值和误差限，绝对误差和相对误差，精度为小数点后4位，其中步长分别取， 9 464, .（2）将（1）中计算的的近似值分别与精确值比较.解 （1）输入程序如下 x=0.79;h=0.1 0.01 0.001 0.0001; M=9464;x1=x+h;x2=x-h; y=sin(5.*x.2-21);y1=sin(5.*x1.2-21); y2=sin(5.*x2.2-21); yz=(y1+y2-2*y)./(h.2), wu=abs(-h.2.*M/24),syms x,f=sin(5.*x.2-21); dy2=diff(f,x,2)运行后屏幕显示利用中心差商公式计算的近似值yz和误差限wu，精度为小数点后4位，其中步长分别取，M=9 464，二阶导函数dy2 如下yz =-45.02889719685589 -45.82347193518466 -45.83048545869772 -45.83055543960768wu =3.94333333333333 0.03943333333333 0.00039433333333 0.00000394333333dy2 =-100*sin(5*x2-21)*x2+10*cos(5*x2-21)（2）输入程序 x=0.79; dy2 =-100*sin(5*x2-21)*x2+10*cos(5*x2-21)wuj=abs(yz-dy2), wux= wuj./abs(yz)运行后屏幕显示利用中心差商公式计算的近似值与精确值的绝对误差wuj，相对误差wux 和的精确值dy2如下dy2 =-45.83055620608437wuj=0.80165900922847 0.00708427089971 0.00007074738664 0.00000076647669wux=0.01780321214006 0.00015459917375 0.00000154367526 0.000000016724148.4 数值梯度和数值偏导数的计算及其MATLAB程序8.4.1 梯度和偏导数的数值计算及其MATLAB程序例8.4.3 设二元函数e，取区域为，和的步长为Hx= Hy=0.2，试求数值梯度向量，并画图.解 输入计算数值梯度向量和画图的程序 x,y = meshgrid(-2.1:.2:2.1, -2.1:.2:2.1); z = 3*x.2 .* exp(-x.2 - y.2);