ch4_例题.doc

【例4.1-1】设，试用机器零阈值eps替代理论0计算极限，。本例演示：除非数值近似法求的极限经过理论验证，否则绝不要借助数值法求取极限。%不可信的“极限的数值近似计算”x=eps;L1=(1-cos(2*x)/(x*sin(x),%将得到一个错误的极限值L2=sin(x)/x,%恰巧与理论值一致 %可信的“极限的符号计算”syms tf1=(1-cos(2*t)/(t*sin(t);f2=sin(t)/t;Ls1=limit(f1,t,0)Ls2=limit(f2,t,0) 说明l 理论分析表明，而且。l 借助符号计算所求的极限与理论值一致。l 而用数值法近似计算的极限与理论不一致。在此，再次提醒读者：不要企图借助数值计算求取极限。【例4.1-2】已知，求。dfdx_p=abs(eps)/epsdfdx_n=abs(-eps)/(-eps) 【例4.1-3】已知，求该函数在区间 中的近似导函数。本例演示：自变量增量的适当取值对数值导函数的精度影响极大。（1）自变量的增量取eps数量级d=pi/100;t=0:d:2*pi;x=sin(t);dt=5*eps;%增量为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 增量过小引起有效数字严重丢失后的毛刺曲线（2）自变量的增量取得适当x_d=sin(t+d);dxdt_d=(x_d-x)/d;%以d=pi/100为增量算得的数值导数plot(t,x,LineWidth,5)hold onplot(t,dxdt_d)hold offlegend(x(t),dx/dt)xlabel(t) 图 4.1-2 增量适当所得导函数比较光滑说明l 当自变量增量太小时，与的数值十分接近，它们的高位有效数字完全相同。这样的高位有效数字消失，导致精度急剧变差。讨论题分别利用增量dt=1,3,5,20,100,500*eps，通过的方式，计算在处的数值导数。并对所得结果给于评述。请观察：计算处的数值导数，与步长的关系。dt=1,3,5,20,100,500*eps;dsin_dt_pi=(sin(pi+dt)-sin(pi)./dtplot(dsin_dt_pi),shgl 在例4.1-1和例4.1-2中，由于是完全准确值0，增量计算没有任何精度损失。l 再次强调指出：数值导数的使用应十分谨慎。l(3)采用diff和gradient指令。clfd=pi/100; %自变量增量 t=0:d:2*pi;x=sin(t);dxdt_diff=diff(x)/d; %diff求得的近似数值导数 dxdt_grad=gradient(x)/d; %gradient求得的近似数值导数subplot(1,2,1)plot(t,x,b)hold onplot(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, 2pi)legend(x(t),dxdt_grad,dxdt_diff,Location,North)xlabel(t),box offhold offsubplot(1,2,2)kk=(length(t)-10):length(t);%数组中最后11个数据的下标hold onplot(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),box offhold off 【例 4.1-5】求积分，其中。本例演示：trapz用于数值积分时的基本原理；sum的用法及注意事项。 cleard=pi/8;%子区间间隔t=0:d:pi/2;%包含5个采样点的一维数组y=0.2+sin(t);%5个点处的函数值数组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;%因采用stairs绘图需要而写y2=y,nan;%因采用stairs绘图需要而写stairs(t2,y2,:k)%虚线下面积表示d*sum的几何意义hold onplot(t,y,r,LineWidth,3)%红折线下面积表示d*trapz的几何意义h=stem(t,y,LineWidth,2);%用蓝空心杆线表示函数采样值set(h(1),MarkerSize,10)axis(0,pi/2+d,0,1.5)%使横坐标恰好是0, 5*dhold offshg sum求得积分 trapz求得积分 1.5762 1.3013 图 4.1-4 sum 和trapz求积模式示意【例 4.1-6】求 。本例演示：精度不可控和可控积分指令及结果的差别；泛函指令中函数的字符串表达法。（1）采用符号计算获得具有32位精度的积分值syms xIsym=vpa(int(exp(-x2),x,0,1) （2）采用trapz计算积分format long%采用15位数字显示计算结果d=0.001;x=0:d:1;Itrapz=d*trapz(exp(-x.*x) （3）采用字符串表达被积函数fx=exp(-x.2);Ic=quad(fx,0,1,1e-8)%控制绝对精度达1e-8 【例 4.1-7】求。本例演示：dblquad指令的调用格式；泛函指令中函数的“匿名函数”表达法。（1）符号计算法syms x ys=vpa(int(int(xy,x,0,1),y,1,2)% （2）数值积分法format longs_n=dblquad(x,y)x.y,0,1,1,2)% 说明l dblquad指令第1个输入量采用“匿名函数”(x,y)x.y写成。注意：被积函数一定要写成“数组运算”格式。l 对于本例而言，下列三条指令的作用相同。s_n=dblquad(x,y)x.y,0,1,1,2)%匿名函数表示s_n=dblquad(x.y,0,1,1,2)%字符串表示s_n=dblquad(inline(x.y),0,1,1,2)%内联对象表示【例4.1-8】已知，在区间，求函数的极小值。本例演示：fminbnd极小值求取法；教科书方法在本例失败；采用图形法求极值。（1）采用高等数学教科书方法求极小值失败syms xy=(x+pi)*exp(abs(sin(x+pi);yd=diff(y,x);%求导函数xs0=solve(yd)%求导函数为0的自变量值xs0yd_xs0=vpa(subs(yd,x,xs0),6)%验算用：导函数在xs0处为0吗？y_xs0=vpa(subs(y,x,xs0),6)%计算y(xs0)，可以发现它大于左边界点处的函数值y_m_pi=vpa(subs(y,x,-pi/2),6)%计算左边界点函数值y(-pi/2)y_p_pi=vpa(subs(y,x,pi/2),6)%计算右边界点函数值y(pi/2) （2）采用优化算法求极小值x1=-pi/2;x2=pi/2;%搜索区间的边界yx=(x)(x+pi)*exp(abs(sin(x+pi);%采用匿名函数形式定义被求极小值的函数y(x)xn0,fval,exitflag,output=fminbnd(yx,x1,x2)%xn0,fval分别是极值点和函数极小值 （3）图形法求极小值%在-pi/2,pi/2区间绘制函数y(x)曲线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)%从局部放大图上取出极值点及相应函数值 图 4.1-6 函数极值点附近的局部放大和交互式取值说明l 本例函数的特点：一，导函数不连续；二，在x = -1附近，导数存在，且为0，是一个极大值。【例4.1-9】求的极小值点。它即是著名的Rosenbrocks Banana 测试函数，它的理论极小值是。本例演示：二元函数极值点的求取；多搜索起点。（1）本例采用匿名函数表示测试函数如下ff=(x)(100*(x(2)-x(1).2)2+(1-x(1)2); 注意：在编写目标函数时，自变量不是采用x, y表示，而是采用一个名为x 的“二元向量”表示的。（2）用单纯形法求极小值点x0=-5,-2,2,5;-5,-2,2,5;%提供4个搜索起点sx,sfval=fminsearch(ff,x0)%sx给出一组使优化函数值非减的局部极小点 （3）检查目标函数值format short edisp(ff(sx(:,1),ff(sx(:,2),ff(sx(:,3),ff(sx(:,4) 【例 4.1-10】求微分方程在初始条件情况下的解，并图示。（见图4.1-7和4.1-8）（1）把高阶微分方程改写成一阶微分方程组令，于是原二阶方程可改写成如下一阶方程组，且设。这是有名的van der Pol Stiff方程，mu值大小决定“刚度”。如mu=1000时，要用ode15, ode23去解。（2）根据上述一阶微分方程组编写M函数文件DyDt.mfunction ydot=DyDt(t,y)mu=2;ydot=y(2);mu*(1-y(1)2)*y(2)-y(1);%注意一阶导数ydot是(2*1)列向量（3）解算微分方程tspan=0,30;%求解的时间区间y0=1;0;%初值向量应与DyDt.m文件中y形式一致。tt,yy=ode45(DyDt,tspan,y0);%plot(tt,yy(:,1)xlabel(t),title(x(t) 图 4.1-7 微分方程解（4）画相平面图（图4.1-8）plot(yy(:,1),yy(:,2)%函数和其导函数勾画的曲线称为“相轨迹”xlabel(位移),ylabel(速度) 图 4.1-8 平面相轨迹说明l 注意第条指令中 DyDt 就是函数文件DyDt.m的函数句柄。注意：DyDt.m文件所在目录必须设置成当前目录，或设置在MATLAB的搜索路径上。【例 4.2-1】矩阵乘法算例。本例演示：矩阵乘法的含义；如何减少和避免循环；如何判断两个双精度矩阵是否相等。（1）实现矩阵乘的“标量程式”clearrand(state,12)A=rand(2,4);B=rand(4,3);%-用三重循环体实现矩阵乘法-C1=zeros(size(A,1),size(B,2);%为C预分配内存，可加快计算for ii=1:size(A,1)for jj=1:size(B,2)for k=1:size(A,2)C1(ii,jj)=C1(ii,jj)+A(ii,k)*B(k,jj);endendend% （2）实现矩阵乘的“saxpy程式”（saxpy是LAPACK术语）%-用二重循环体实现矩阵乘法-C2=zeros(size(A,1),size(B,2);for jj=1:size(B,2)for k=1:size(B,1)C2(:,jj)=C2(:,jj)+