4.1数值微积分4.1.1近似数值极限及导数Matlab 数值计算中，没有求极限指令，也没有求导指令，而是利用差分指令：用一个简单矩阵表现diff和gradient指令计算方式。

差分： Dx=diff(X)对向量： Dx=X(2:n)-X(1:n-1)对矩阵： DX=X(2:n,:)-X(1:n-1,:)长度小1.DIFF(X), for a vector X, is[X(2)-X(1) X(3)-X(2) ... X(n)-X(n-1)]. DIFF(X), for a matrix X, is the matrix of row differences, （结果缺少一行）[X(2:n,:) - X(1:n-1,:)].DIFF(X,N,DIM) is the Nth difference function along dimension DIM.If N >= size(X,DIM), DIFF returns an empty array （N阶差分）梯度：FX=gradient(F)Fx(1)=Fx(2)-Fx(1);F=[1,2,3;4,5,6;7,8,9]Dx=diff(F) （按行）Dx_2=diff(F,1,2) （按列）[FX,FY]=gradient(F)Fx(1)=Fx(2)-Fx(1),Fx(end)=F(end)-F(end-1)FX与F维数相同。

[FX_2,FY_2]=gradient(F,0.5)%采样间隔0.5 即：Fx(1)=(Fx(2)-Fx(1))/2F =1 2 34 5 67 8 9Dx =3 3 33 3 3Dx_2 =1 11 11 1FX =1 1 11 1 1 1 1 1 （列方向） FY =3 3 3 3 3 3 3 3 3 （行方向） FX_2 =2 2 2 2 2 2 2 2 2 FY_2 =6 6 6 6 6 6 6 6 6【例4.1-1】设xx xx f sin 2cos 1)(1-=，xxx fsin )(2=，试用机器零阈值eps 替代理论0计算极限)(lim )0(11x f L x →=，)(lim )0(22x fL x →=。

x=eps;L1=(1-cos(2*x))/(x*sin(x)), L2=sin(x)/x, L1 =(注意错误！ 数值法求极限得到错误结果) L2 = 1syms tf1=(1-cos(2*t))/(t*sin(t));f2=sin(t)/t;Ls1=limit(f1,t,0)Ls2=limit(f2,t,0)Ls1 =2Ls2 =1【例4.1-2】已知)2,0[π中的近x=，求该函数在区间]sin(t似导函数。

d=pi/100;t=0:d:2*pi;x=sin(t);dt=5*eps;x_eps=sin(t+dt);dxdt_eps=(x_eps-x)/dt;plot(t,x,'LineWidth',5)hold onplot(t,dxdt_eps)hold offlegend('x(t)','dx/dt')xlabel('t')图 4.1-1 增量过小引起有效数字严重丢失后的毛刺曲线x_d=sin(t+d);dxdt_d=(x_d-x)/d;plot(t,x,'LineWidth',5)hold onplot(t,dxdt_d)hold offlegend('x(t)','dx/dt')xlabel('t')图 4.1-2 增量适当所得导函数比较光滑【例4.1-3】已知)x=，采用diff和gradient计算该sin(t函数在区间]2,0[π中的近似导函数。

clfd=pi/100;t=0:d:2*pi;x=sin(t);dxdt_diff=diff(x)/d;dxdt_grad=gradient(x)/d;subplot(1,2,1)plot(t,x,'b')hold onplot(t,dxdt_grad,'m','LineWidth',8) plot(t(1:end-1),dxdt_diff,'.k','Marker Size',8)axis([0,2*pi,-1.1,1.1])title('[0, 2\pi]')legend('x(t)','dxdt_{grad}','dxdt_{dif f}','Location','North')xlabel('t'),box offhold offsubplot(1,2,2)kk=(length(t)-10):length(t);hold onplot(t(kk),dxdt_grad(kk),'om','MarkerS ize',8)plot(t(kk-1),dxdt_diff(kk-1),'.k','Mar kerSize',8)title('[end-10, end]')legend('dxdt_{grad}','dxdt_{diff}','Lo cation','SouthEast')xlabel('t'),box offhold off图 4.1-3 diff和gradient求数值近似导数的异同比较4.1.2数值求和与近似数值积分【例 4.1-4】求积分⎰=2/) ()(πdttyxs，其中)sin(2.0ty+=。

cleard=pi/8;t=0:d:pi/2;y=0.2+sin(t);s=sum(y);s_sa=d*s;s_ta=d*trapz(y);disp(['sum求得积分',blanks(3),'trapz求得积分'])disp([s_sa, s_ta])t2=[t,t(end)+d];y2=[y,nan];stairs(t2,y2,':k')hold onplot(t,y,'r','LineWidth',3) h=stem(t,y,'LineWidth',2); set(h(1),'MarkerSize',10) axis([0,pi/2+d,0,1.5])hold offshgsum求得积分 trapz求得积分1.5762 1.3013图 4.1-4 sum 和trapz求积模式示意4.1.3计算精度可控的数值积分一元函数的数值积分函数1 quad、quadl、quad8功能数值定积分，自适应Simpleson积分法。

格式q = quad(fun,a,b) %近似地从a到b计算函数fun的数值积分，误差为10-6。

若给fun输入向量x，应返回向量y，即fun是一单值函数。

q = quad(fun,a,b,tol) %用指定的绝对误差tol代替缺省误差。

tol越大，函数计算的次数越少，速度越快，但结果精度变小。

[q,n] = quad(fun,a,b,…) %同时返回函数计算的次数n…= quadl(fun,a,b,…) %用高精度进行计算，效率可能比quad更好。

…= quad8(fun,a,b,…) %该命令是将废弃的命令，用quadl代替。

例：>>fun = inline(‘’3*x.^2./(x.^3-2*x.^2+3)’);>>Q1 = quad(fun,0,2)>>Q2 = quadl(fun,0,2)计算结果为：Q1 =3.7224Q2 =3.7224函数2 trapz功能梯形法数值积分格式T = trapz(Y)%用等距梯形法近似计算Y的积分。

若Y是一向量，则trapz(Y)为Y的积分；若Y是一矩阵，则trapz(Y)为Y的每一列的积分；若Y是一多维阵列，则trapz(Y)沿着Y的第一个非单元集的方向进行计算。

T = trapz(X,Y)%用梯形法计算Y在X点上的积分。

若X为一列向量，Y为矩阵，且size(Y,1) =length(X)，则trapz(X,Y)通过Y的第一个非单元集方向进行计算。

T = trapz(…,dim) %沿着dim指定的方向对Y进行积分。

若参量中包含X，则应有length(X)=size(Y ,dim)。

例：>>X = -1:.1:1;>>Y = 1./(1+25*X.^2);>>T = trapz(X,Y)计算结果为：T =0.5492二重积分S2=ablquad(fun,xmin,xmax,ymin,ymax,tol)Tol 来控制绝对误差。

默认10^{-6}。

【例 4.1-5】求dx e I x ⎰-=102 。

syms xIsym=vpa(int(exp(-x^2),x,0,1)) Isym =.74682413281242702539946743613185format longd=0.001;x=0:d:1;Itrapz=d*trapz(exp(-x.*x))Itrapz =0.746824071499185fx='exp(-x.^2)';Ic=quad(fx,0,1,1e-8)Ic =0.746824132854452【例 4.1-6】求dxdy x s y ⎰⎰=21 1 0 。

syms x ys=vpa(int(int(x^y,x,0,1),y,1,2))s =.40546510810816438197801311546432format longs_n=dblquad(@(x,y)x.^y,0,1,1,2)s_n =0.4054662672435084.1.4函数极值的数值求解【例4.1-7】已知)sin()(ππ+⋅+=x ex y ，在2/2/ππ≤≤-x 区间，求函数的极小值。

（1）syms x y=(x+pi)*exp(abs(sin(x+pi)));yd=diff(y,x);xs0=solve(yd) yd_xs0=vpa(subs(yd,x,xs0),6)y_xs0=vpa(subs(y,x,xs0),6)y_m_pi=vpa(subs(y,x,-pi/2),6)y_p_pi=vpa(subs(y,x,pi/2),6)Warning: Warning, solutions may have been lostxs0 =-1.0676598444985783372948670854801yd_xs0 =0.y_xs0 =4.98043y_m_pi =4.26987y_p_pi =12.8096（2）x1=-pi/2;x2=pi/2;yx=@(x)(x+pi)*exp(abs(sin(x+pi))); [xn0,fval,exitflag,output]=fminbnd(yx, x1,x2)xn0 =-1.2999e-005fval =3.1416exitflag =1output =iterations: 21funcCount: 22algorithm: 'golden section search, parabolic interpolation'message: [1x112 char]（3）xx=-pi/2:pi/200:pi/2;yxx=(xx+pi).*exp(abs(sin(xx+pi)));plot(xx,yxx)xlabel('x'),grid on图 4.1-5 在[-pi/2,pi/2]区间中的函数曲线[xx,yy]=ginput(1)xx=1.5054e-008yy=3.1416图 4.1-6 函数极值点附近的局部放大和交互式取值【例4.1-8】求22)2xyf--=的极小值点。