MATLAB数值微积分

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),foravectorX,is[X(2)-X(1)X(3)-X(2)...X(n)-X(n-1)].DIFF(X),foramatrixX,isthematrixofrowdifferences,（结果缺少一行）[X(2:n,:)-X(1:n-1,:)].DIFF(X,N,DIM)istheNthdifferencefunctionalongdimensionDIM.IfN=size(X,DIM),DIFFreturnsanemptyarray（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=123456789Dx=333333Dx_2=111111FX=111111111（列方向）FY=333333333（行方向）FX_2=222222222FY_2=666666666【例4.1-1】设xxxxfsin2cos1)(1，xxxfsin)(2，试用机器零阈值eps替代理论0计算极限)(lim)0(101xfLx，)(lim)0(202xfLx。x=eps;L1=(1-cos(2*x))/(x*sin(x)),L2=sin(x)/x,L1=0(注意错误！数值法求极限得到错误结果)L2=1symstf1=(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】已知)sin(tx，求该函数在区间]2,0[中的近似导函数。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)holdonplot(t,dxdt_eps)holdofflegend('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)holdonplot(t,dxdt_d)holdofflegend('x(t)','dx/dt')xlabel('t')图4.1-2增量适当所得导函数比较光滑【例4.1-3】已知)sin(tx，采用diff和gradient计算该函数在区间]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')holdonplot(t,dxdt_grad,'m','LineWidth',8)plot(t(1:end-1),dxdt_diff,'.k','MarkerSize',8)axis([0,2*pi,-1.1,1.1])title('[0,2\pi]')legend('x(t)','dxdt_{grad}','dxdt_{diff}','Location','North')xlabel('t'),boxoffholdoffsubplot(1,2,2)kk=(length(t)-10):length(t);holdonplot(t(kk),dxdt_grad(kk),'om','MarkerSize',8)plot(t(kk-1),dxdt_diff(kk-1),'.k','MarkerSize',8)title('[end-10,end]')legend('dxdt_{grad}','dxdt_{diff}','Location','SouthEast')xlabel('t'),boxoffholdoff图4.1-3diff和gradient求数值近似导数的异同比较4.1.2数值求和与近似数值积分【例4.1-4】求积分2/0)()(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')holdonplot(t,y,'r','LineWidth',3)h=stem(t,y,'LineWidth',2);set(h(1),'MarkerSize',10)axis([0,pi/2+d,0,1.5])holdoffshgsum求得积分trapz求得积分1.57621.3013图4.1-4sum和trapz求积模式示意4.1.3计算精度可控的数值积分一元函数的数值积分函数1quad、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函数2trapz功能梯形法数值积分格式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】求dxeIx102。symsxIsym=vpa(int(exp(-x^2),x,0,1))Isym=.74682413281242702539946743613185formatlongd=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】求dxdyxsy2110。symsxys=vpa(int(int(x^y,x,0,1),y,1,2))s=.40546510810816438197801311546432formatlongs_n=dblquad(@(x,y)x.^y,0,1,1,2)s_n=0.4054662672435084.1.4函数极值的数值求解【例4.1-7】已知)sin()(xexy，在2/2/x区间，求函数的极小值。（1）symsxy=(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,solutionsmayhavebeenlostxs0=-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:'goldensectionsearch,parabolicinterpolation'message:[1x112char]（3）xx=-pi/2:pi/200:pi/2;yxx=(xx+pi).*exp(abs(sin(xx+pi)));plot(xx,yxx)xlabel('x'),gridon-2-1.5-1-0.500.511.52345678910111213x图4.1-5在[-pi/2,pi/2]区间中的函数曲线[xx,yy]=ginput(1)xx=1.5054e-008yy=3.1416图4.1-6函数极值点附近的局部放大和交互式取值【例4.1-8】求222)1()(100),(xxyyxf的极小值点。它即是著名的Rosenbrock'sBanana测试函数，它的理论极小值是1,1yx。（1）ff=@(x)(100*(x(2)-x(1).^2)^2+(1-x(1))^2);（2）x0=[-5,-2,2,5;-5,-2,2,5];[sx,sfval,sexit,soutput]=fminsearch(ff,x0)sx=1.0000-0.68970.41518.08861.0000-1.91684.96437.8004sfval=2.4112e-010sexit=1soutput=iterations:384funcCount:615algorithm:[1x33char]message:[1x196char]（3）检查目标函数值formatshortedisp([ff(sx(:,1)),ff(sx(:,2)),ff(sx(:,3)),ff(sx(:,4))])Columns1through32.4112e-0105.7525e+0022.2967e+003Column43.3211e+0054.1.5常微分方程的数值解指令：[t,y]=ode45(fun,tspan,x0)【例4.1-9】求微分方程0)1(222xdtdxxdtxd，2，在初始条件0)0(,1)0(dtd