matlab课后习题答案第四章.doc

《matlab课后习题答案第四章.doc》由会员分享，可在线阅读，更多相关《matlab课后习题答案第四章.doc（49页珍藏版）》请在得力文库 - 分享文档赚钱的网站上搜索。

1、第4章 数值运算习题 4 及解答1 根据题给的模拟实际测量数据的一组和 试用数值差分diff或数值梯度gradient指令计算，然后把和曲线绘制在同一张图上，观察数值求导的后果。（模拟数据从prob_data401.mat获得）目的l 强调：要非常慎用数值导数计算。l 练习mat数据文件中数据的获取。l 实验数据求导的后果l 把两条曲线绘制在同一图上的一种方法。解答（1）从数据文件获得数据的指令假如prob_data401.mat文件在当前目录或搜索路径上clearload prob_data401.mat （2）用diff求导的指令dt=t(2)-t(1);yc=diff(y)/dt;%注意

2、yc的长度将比y短1plot(t,y,b,t(2:end),yc,r)grid on （3）用gradent求导的指令（图形与上相似）dt=t(2)-t(1);yc=gradient(y)/dt;plot(t,y,b,t,yc,r)grid on 说明l 不到万不得已，不要进行数值求导。l 假若一定要计算数值导数，自变量增量dt 要取得比原有数据相对误差高1、2个量级以上。l 求导会使数据中原有的噪声放大。2 采用数值计算方法，画出在区间曲线，并计算。提示l 指定区间内的积分函数可用cumtrapz指令给出。l 在计算要求不太高的地方可用find指令算得。目的l 指定区间内的积分函数的数值计算

3、法和cumtrapz指令。l find指令的应用。解答dt=1e-4;t=0:dt:10;t=t+(t=0)*eps;f=sin(t)./t;s=cumtrapz(f)*dt;plot(t,s,LineWidth,3)ii=find(t=4.5);s45=s(ii) s45 = 1.6541 3 求函数的数值积分，并请采用符号计算尝试复算。提示l 数值积分均可尝试。l 符号积分的局限性。目的l 符号积分的局限性。解答dx=pi/2000;x=0:dx:pi;s=trapz(exp(sin(x).3)*dx s = 5.1370 符号复算的尝试syms xf=exp(sin(x)3);ss=in

4、t(f,x,0,pi) Warning: Explicit integral could not be found. In sym.int at 58ss =int(exp(sin(x)3),x = 0 . pi) 4 用quad求取的数值积分，并保证积分的绝对精度为。目的l quadl，精度可控，计算较快。l 近似积分指令trapz获得高精度积分的内存和时间代价较高。解答%精度可控的数值积分fx=(x)exp(-abs(x).*abs(sin(x);format longsq=quadl(fx,-10*pi,1.7*pi,1e-7) sq = 1.08784993815498 %近似积分算法

5、x=linspace(-10*pi,1.7*pi,1e7);dx=x(2)-x(1);st=trapz(exp(-abs(x).*abs(sin(x)*dx st = 1.08784949973430 %符号积分算法y=exp(-abs(x)*abs(sin(x)si=vpa(int(y,-10*pi,1.7*pi),16) y =exp(-abs(x)*abs(sin(x)si =1.087849499412911 5 求函数在区间中的最小值点。目的l 理解极值概念的邻域性。l 如何求最小值。l 学习运用作图法求极值或最小值。l 感受符号法的局限性。解答（1）采用fminbnd找极小值点在指

6、令窗中多次运行以下指令，观察在不同数目子区间分割下，进行的极小值搜索。然后从一系列极小值点中，确定最小值点。clearft=(t)sin(5*t).2.*exp(0.06*t.*t)+1.8*abs(t+0.5)-1.5*t.*cos(2*t);disp(计算中，把 -5,5 分成若干搜索子区间。)N=input( 请输入子区间数 N，注意使N=1 ？);%该指令只能在指令窗中运行tt=linspace(-5,5,N+1);for k=1:Ntmin(k),fobj(k)=fminbnd(ft,tt(k),tt(k+1);endfobj,ii=sort(fobj);%将目标值由小到大排列tmi

7、n=tmin(ii);%使极小值点做与目标值相应的重新排列fobj,tmin（2）最后确定的最小值点在的不同分割下，经观察，最后确定出最小值点是 -1.28498111480531相应目标值是-0.18604801006545（3）采用作图法近似确定最小值点（另一方法）（A）在指令窗中运行以下指令：clearft=(t)sin(5*t).2.*exp(0.06*t.*t)+1.8*abs(t+0.5)-1.5*t.*cos(2*t);t=-5:0.001:5;ff=ft(t);plot(t,ff)grid on,shg（B）经观察后，把最小值附近邻域放到足够大，然后运行以下指令，那放大图形被推

8、向前台，与此同时光标变为“十字线”，利用它点击极值点可得到最小值数据tmin2,fobj2=ginput(1) tmin2 = -1.28500000993975fobj2 = -0.18604799369136 出现具有相同数值的刻度区域表明已达最小可分辨状态（4）符号法求最小值的尝试syms tfts=sin(5*t)2*exp(0.06*t*t)-1.5*t*cos(2*t)+1.8*abs(t+0.5);dfdt=diff(fts,t);%求导函数tmin=solve(dfdt,t)%求导函数的零点fobj3=subs(fts,t,tmin)%得到一个具体的极值点 tmin =-.60

9、100931947716486053884417850955e-2fobj3 =.89909908144684551670208797723124 说明l 最小值是对整个区间而言的，极小值是对邻域而言的。l 在一个区间中寻找最小值点，对不同子区间分割进行多次搜索是必要的。这样可以避免把极小值点误作为最小值点。最小值点是从一系列极小值点和边界点的比较中确定的。l 作图法求最小值点，很直观。假若绘图时，自变量步长取得足够小，那么所求得的最小值点有相当好的精度。l 符号法在本例中，只求出一个极值点。其余很多极值点无法秋初，更不可能得到最小值。6 设，用数值法和符号法求。目的l 学习如何把高阶微分方程

10、写成一阶微分方程组。l ode45解算器的导数函数如何采用匿名函数形式构成。l 如何从ode45一组数值解点，求指定自变量对应的函数值。解答（1）改写高阶微分方程为一阶微分方程组令，于是据高阶微分方程可写出（2）运行以下指令求y(t)的数值解format longts=0,1;y0=1;0;dydt=(t,y)y(2);-2*y(1)+3*y(2)+1;%匿名函数写成的ode45所需得导数函数tt,yy=ode45(dydt,ts,y0); y_05=interp1(tt,yy(:,1),0.5,spline), %用一维插值求y(0.5) y_05 = 0.78958020790127 （3

11、）符号法求解syms t;ys=dsolve(D2y-3*Dy+2*y=1,y(0)=1,Dy(0)=0,t)ys_05=subs(ys,t,sym(0.5) ys =1/2-1/2*exp(2*t)+exp(t)ys_05 =.78958035647060552916850705213780 说明l 第条指令中的导数函数也可采用M函数文件表达，具体如下。function S=prob_DyDt(t,y)S=y(2);-2*y(1)+3*y(2)+1; 7 已知矩阵A=magic(8)，（1）求该矩阵的“值空间基阵”B ；（2）写出“A的任何列可用基向量线性表出”的验证程序（提示：利用rref

12、检验）。目的l 体验矩阵值空间的基向量组的不唯一性，但它们可以互为线性表出。l 利用rref检验两个矩阵能否互为表出。解答（1）A的值空间的三组不同“基”A=magic(8);%采用8阶魔方阵作为实验矩阵R,ci=rref(A);B1=A(:,ci)%直接从A中取基向量B2=orth(A)%求A值空间的正交基V,D=eig(A);rv=sum(sum(abs(D)1000*eps);%非零特征值数就是矩阵的秩B3=V(:,1:rv)%取A的非零特征值对应的特征向量作基 B1 = 64 2 3 9 55 54 17 47 46 40 26 27 32 34 35 41 23 22 49 15 1

13、4 8 58 59B2 = -0.3536 0.5401 0.3536 -0.3536 -0.3858 -0.3536 -0.3536 -0.2315 -0.3536 -0.3536 0.0772 0.3536 -0.3536 -0.0772 0.3536 -0.3536 0.2315 -0.3536 -0.3536 0.3858 -0.3536 -0.3536 -0.5401 0.3536B3 = 0.3536 0.6270 0.3913 0.3536 -0.4815 -0.2458 0.3536 -0.3361 -0.1004 0.3536 0.1906 -0.0451 0.3536 0.

14、0451 -0.1906 0.3536 0.1004 0.3361 0.3536 0.2458 0.4815 0.3536 -0.3913 -0.6270 （2）验证A的任何列可用B1线性表出B1_A=rref(B1,A)%若B1_A矩阵的下5行全为0，%就表明A可以被B1的3根基向量线性表出 B1_A = 1 0 0 1 0 0 1 1 0 0 1 0 1 0 0 1 0 3 4 -3 -4 7 0 0 1 0 0 1 -3 -4 4 5 -7 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

15、0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 B2_A=rref(B2,A) B2_A = Columns 1 through 7 1.0000 0 0 -91.9239 -91.9239 -91.9239 -91.9239 0 1.0000 0 51.8459 -51.8459 -51.8459 51.8459 0 0 1.0000 9.8995 -7.0711 -4.2426 1.4142 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Columns 8 through

16、11 -91.9239 -91.9239 -91.9239 -91.9239 51.8459 -51.8459 -51.8459 51.8459 -1.4142 4.2426 7.0711 -9.8995 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 B3_A=rref(B3,A) B3_A = Columns 1 through 7 1.0000 0 0 91.9239 91.9239 91.9239 91.9239 0 1.0000 0 42.3447 -38.1021 -33.8594 29.6168 0 0 1.0000 12.6462 -16.888

17、9 -21.1315 25.3741 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Columns 8 through 11 91.9239 91.9239 91.9239 91.9239 25.3741 -21.1315 -16.8889 12.6462 29.6168 -33.8594 -38.1021 42.3447 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 说明l magic(n)产生魔方阵。魔方阵具有很多特异的性质。就其秩而言，当n为奇数时，该矩阵满秩；

18、当n 为4的倍数时，矩阵的秩总是3；当 n 为偶数但不是4倍数时，则矩阵的秩等于 (n/2+2)。关于魔方阵的有关历史，请见第6.1.3节。8 已知由MATLAB指令创建的矩阵A=gallery(5)，试对该矩阵进行特征值分解，并通过验算观察发生的现象。目的l 展示特征值分解可能存在的数值问题。l condeig是比较严谨的特征值分解指令。l Jordan分解的作用。解答（1）特征值分解A=gallery(5)V,D=eig(A);diag(D)%为紧凑地显示特征值而写A = -9 11 -21 63 -252 70 -69 141 -421 1684 -575 575 -1149 3451

19、-13801 3891 -3891 7782 -23345 93365 1024 -1024 2048 -6144 24572ans = Columns 1 through 4 -0.0181 -0.0054 - 0.0171i -0.0054 + 0.0171i 0.0144 - 0.0104i Column 5 0.0144 + 0.0104i （2）验算表明相对误差较大AE=V*D/Ver_AE=norm(A-AE,fro)/norm(A,fro)%相对F范数 AE = 1.0e+004 * Columns 1 through 4 -0.0009 + 0.0000i 0.0011 - 0

20、.0000i -0.0021 + 0.0000i 0.0063 - 0.0000i 0.0070 - 0.0000i -0.0069 + 0.0000i 0.0141 - 0.0000i -0.0421 + 0.0000i -0.0575 + 0.0000i 0.0575 - 0.0000i -0.1149 + 0.0000i 0.3451 - 0.0000i 0.3891 - 0.0000i -0.3891 + 0.0000i 0.7781 - 0.0000i -2.3343 + 0.0000i 0.1024 - 0.0000i -0.1024 + 0.0000i 0.2048 - 0.00

21、00i -0.6144 + 0.0000i Column 5 -0.0252 + 0.0000i 0.1684 - 0.0000i -1.3800 + 0.0000i 9.3359 - 0.0001i 2.4570 - 0.0000ier_AE = 6.9310e-005 （3）一个更严谨的特征值分解指令Vc,Dc,eigc=condeig(A)%eigc中的高值时，说明相应的特征值不可信。 Vc = Columns 1 through 4 -0.0000 -0.0000 + 0.0000i -0.0000 - 0.0000i 0.0000 + 0.0000i 0.0206 0.0207 +

22、0.0000i 0.0207 - 0.0000i 0.0207 + 0.0000i -0.1397 -0.1397 + 0.0000i -0.1397 - 0.0000i -0.1397 + 0.0000i 0.9574 0.9574 0.9574 0.9574 0.2519 0.2519 - 0.0000i 0.2519 + 0.0000i 0.2519 - 0.0000i Column 5 0.0000 - 0.0000i 0.0207 - 0.0000i -0.1397 - 0.0000i 0.9574 0.2519 + 0.0000iDc = Columns 1 through 4 -

23、0.0181 0 0 0 0 -0.0054 + 0.0171i 0 0 0 0 -0.0054 - 0.0171i 0 0 0 0 0.0144 + 0.0104i 0 0 0 0 Column 5 0 0 0 0 0.0144 - 0.0104ieigc = 1.0e+011 * 5.2687 5.2313 5.2313 5.1725 5.1724 （4）对A采用Jordan分解并检验VJ,DJ=jordan(A);%求出准确的广义特征值，使A*VJ=VJ*D成立。DJAJ=VJ*DJ/VJer_AJ=norm(A-AJ,fro)/norm(A,fro) DJ = 0 1 0 0 0 0

24、0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0AJ = 1.0e+004 * -0.0009 0.0011 -0.0021 0.0063 -0.0252 0.0070 -0.0069 0.0141 -0.0421 0.1684 -0.0575 0.0575 -0.1149 0.3451 -1.3801 0.3891 -0.3891 0.7782 -2.3345 9.3365 0.1024 -0.1024 0.2048 -0.6144 2.4572er_AJ = 2.0500e-011 说明l 指令condeig的第3输出量eigc给出的是所谓的“矩阵特征值条件数”。

25、当特征条件数与相当时，就意味着矩阵A可能“退化”，即矩阵可能存在“代数重数”大于“几何重数”的特征值。此时，实施Jordan分解更适宜。l 顺便指出：借助condeig算得的特征值条件数与cond指令算得的矩阵条件数是两个不同概念。前者描述特征值的问题，后者描述矩阵逆的问题。l 本例矩阵A的特征值条件数很高，表明分解不可信。验算也表明，相对误差较大。l 当对矩阵A进行Jordan分解时，可以看到，A具有5重根。当对Jordan分解进行验算时，相对误差很小。9 求矩阵的解，A为3阶魔方阵，b是的全1列向量。提示l 了解magic指令l rref 用于方程求解。l 矩阵除法和逆阵法解方程。目的l

26、满秩方阵求解的一般过程。l rref 用于方程求解。l 矩阵除法和逆阵法解方程。解答A=magic(3);%产生3阶魔方阵b=ones(3,1);%（3*1）全1列向量R,C=rref(A,b)%Gauss Jordan消去法解方程，同时判断解的唯一性x=Ab%矩阵除解方程xx=inv(A)*b%逆阵法解方程 R = 1.0000 0 0 0.0667 0 1.0000 0 0.0667 0 0 1.0000 0.0667C = 1 2 3x = 0.0667 0.0667 0.0667xx = 0.0667 0.0667 0.0667 说明l rref指令通过对增广矩阵进行消去法操作完成解方

27、程。由分解得到的3根“坐标向量”和（或）C3指示的3根基向量，可见A3满秩，因此方程解唯一。l 在本例情况下，矩阵除、逆阵法、rref法所得解一致。10 求矩阵的解，A为4阶魔方阵，b是的全1列向量。提示l 用rref 可观察A不满秩，但b在A的值空间中，这类方程用无数解。l 矩阵除法能正确求得这类方程的特解。l 逆阵法不能求得这类方程的特解。l 注意特解和齐次解目的l A不满秩，但b在A的值空间中，这类方程的求解过程。l rref 用于方程求解。l 矩阵除法能正确求得这类方程的特解。l 逆阵法不能求得这类方程的特解。l 解的验证方法。l 齐次解的获取。l 全解的获得。解答（1）借助增广矩阵用

28、指令rref求解A=magic(4);%产生3阶魔方阵b=ones(4,1);%全1列向量R,C=rref(A,b)%求解，并判断解的唯一性 R = 1.0000 0 0 1.0000 0.0588 0 1.0000 0 3.0000 0.1176 0 0 1.0000 -3.0000 -0.0588 0 0 0 0 0C = 1 2 3 关于以上结果的说明：l R阶梯阵提供的信息n 前4列是原A阵经消元变换后的阶梯阵；而第5列是原b向量经相同变换后的结果。n R的前三列为“基”，说明原A阵秩为3；而第4列的前三个元素，表示原A阵的第4列由其前三列线性组合而成时的加权系数，即方程的一个解。n

29、R的第5列表明：b可由原A阵的前三列线性表出；b给出了方程的一个解；由于原A阵“缺秩”，所以方程的确解不唯一。l C数组提供的信息n 该数组中的三个元素表示变换取原A阵的第1，2，3列为基。n 该数组的元素总数就是“原A阵的秩”（2）矩阵除求得的解x=Ab Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.306145e-017.x = 0.0588 0.1176 -0.0588 0 运行结果指示：矩阵除法给出的解与rref解相同。（实际上，MATLAB在设计“除

30、法”程序时，针对“b在A值空间中”的情况，就是用rref求解的。）（3）逆阵法的解xx=inv(A)*b Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.306145e-017.xx = 0.0469 0.1875 -0.0625 -0.0156 （3）验证前面所得的解b_rref=A(:,C)*R(C,5)%验算rref的解b_d=A*x%验算矩阵除的解b_inv=A*xx%验算逆阵法的解 b_rref = 1 1 1 1b_d = 1 1 1 1b_inv

31、= 0.7344 1.5469 1.1719 1.8594 显然，在本例中，逆阵法的解是错误的。原因是：A不满秩，A的逆阵在理论上不存在。这里所给出的逆阵是不可信的。（4）求齐次解xg=null(A)%Ax=0的齐次解 xg = 0.2236 0.6708 -0.6708 -0.2236 （5）方程的全解齐次解的任何倍与特解之和就构成方程的全解。下面通过一组随机系数验证。rng default%为本书结果可被读者核对而设，并非必要。f=randn(1,6)%6个随机系数xx=repmat(x,1,6)+xg*f%产生6个不同的特解A*xx%所得结果的每列都应该是全1，即等于b. f = 0.5

32、377 1.8339 -2.2588 0.8622 0.3188 -1.3077xx = 0.1790 0.4689 -0.4463 0.2516 0.1301 -0.2336 0.4783 1.3479 -1.3976 0.6960 0.3315 -0.7596 -0.4195 -1.2890 1.4565 -0.6372 -0.2727 0.8184 -0.1202 -0.4101 0.5051 -0.1928 -0.0713 0.2924ans = 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000

33、1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 解答l （在用除法和逆阵法求解时出现）警告信息中RCOND = 1.306145e-017是矩阵A的估计条件倒数。该数愈接近0，A就愈“病态”；该数接近1时，A就愈“良态”。该条件数由rcond(A)给出。注意：rcond条件倒数与cond条件数的算法不同。11 求矩阵的解，A为4阶魔方阵，。提示l 由rref可以看出A不满秩，b不在A的值空间中，方程没有准确解。l 但可求最小二乘近似解。目的l A不

34、满秩，b不在A的值空间中，方程没有准确解。解答（1）借助增广矩阵用指令rref求解A=magic(4);%产生3阶魔方阵b=(1:4);R,C=rref(A,b)%求解，并判断解的唯一性 R = 1 0 0 1 0 0 1 0 3 0 0 0 1 -3 0 0 0 0 0 1C = 1 2 3 5 （2）用伪逆求最小二乘近似解（超出范围，仅供参考。）x=pinv(A)*b%非准确解b_pinv=A*x%验算 x = 0.0235 0.1235 0.1235 0.0235b_pinv = 1.3000 2.9000 2.1000 3.7000 解答l C表明，A的秩为3，A不满秩；R第5列第4元素非零，说明b不在A的值空