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

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

2、)/dt; %注意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 说明不到万不得已，不要进行数值求导。假若一定要计算数值导数，自变量增量dt 要取得比原有数据相对误差高1、2个量级以上。求导会使数据中原有的噪声放大。12 采用数值计算方法，画出在区间曲线，并计算。提示指定区间内的积分函数可用cumtrapz指令给出。在计算要求不太高的地方可用find指令算得。目的指定区间内的积分函数的数值计算法和

3、cumtrapz指令。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 13 求函数的数值积分，并请采用符号计算尝试复算。提示数值积分均可尝试。符号积分的局限性。目的符号积分的局限性。解答dx=pi/2000;x=0:dx:pi;s=trapz(exp(sin(x).3)*dx s = 5.1370 符号复算的尝试syms xf=exp(sin(x)3);ss=int(f,x,0,

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

5、-10*pi,1.7*pi,1e7);dx=x(2)-x(1);st=trapz(exp(-abs(x).*abs(sin(x)*dx st =.0878* %符号积分算法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 15 求函数在区间中的最小值点。目的理解极值概念的邻域性。如何求最小值。学习运用作图法求极值或最小值。感受符号法的局限性。解答（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:N tmin(k),fobj(k)=fminbnd(ft,tt(k),tt(k+1); end fobj,ii=sort(fobj); %将目标值由小到大排列 tmin=tmin(ii); %使极小值点做与目

7、标值相应的重新排列fobj,tmin（2）最后确定的最小值点 在的不同分割下，经观察，最后确定出最小值点是 -1.28498111480531相应目标值是 -0.186*545（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.186*136 出现具有相同数值的刻度区域表明已达最小可分辨状态（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 =-.60100931947716486053884417850

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

10、如何从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）符号法求解syms t;ys=dsolve(D2y-3*Dy+2*y=

11、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 说明第条指令中的导数函数也可采用M函数文件表达，具体如下。function S=prob_DyDt(t,y)S=y(2);-2*y(1)+3*y(2)+1; 17 已知矩阵A=magic(8)，（1）求该矩阵的“值空间基阵”B ；（2）写出“A的任何列可用基向量线性表出”的验证程序（提示：利用rref检验）。目的体验矩阵值空间的基向量组的不唯一性，但它们可以互为线性表出。利

12、用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 14 8 58 59B2 = -0.3536 0.5401 0.3536

13、 -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.0451 -0.1906 0.3536 0.1004 0.3361 0

14、.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 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

15、 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 11 -91.9239 -91.9239 -91.9239 -91

16、.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.8889 -21.1315 25.3741 0 0 0 0 0 0 0

17、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 说明magic(n)产生魔方阵。魔方阵具有很多特异的性质。就其秩而言，当n为奇数时，该矩阵满秩；当n 为4的倍数时，矩阵的秩总是3；当 n 为偶数但不是4倍数时，则矩

18、阵的秩等于 (n/2+2)。关于魔方阵的有关历史，请见第6.1.3节。18 已知由MATLAB指令创建的矩阵A=gallery(5)，试对该矩阵进行特征值分解，并通过验算观察发生的现象。目的展示特征值分解可能存在的数值问题。condeig是比较严谨的特征值分解指令。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 -13801 3891 -3891 7782 -23345 93365 102

19、4 -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.0000i -0.0021 + 0.0000i 0.0063 - 0.00

20、00i 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.0000i -0.6144 + 0.0000i Column 5 -0.0252

21、 + 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 + 0.0000i 0.0207 - 0.0000i 0.0207 + 0.0

22、000i -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 -0.0181 0 0 0 0 -0.0054 + 0.0171i 0 0

23、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 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0

24、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 说明指令condeig的第3输出量eigc给出的是所谓的“矩阵特征值条件数”。当特征条件数与相当时，就意味着矩阵A可能“退化”，即矩阵可能存在“代数重数”

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

26、产生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 说明rref指令通过对增广矩阵进行消去法操作完成解方程。由分解得到的3根“坐标向量”和（或）C3指示的3根基向量，可见A3满秩，因此方程解唯一。在本例情况下

27、，矩阵除、逆阵法、rref法所得解一致。110 求矩阵的解，A为4阶魔方阵，b是的全1列向量。提示用rref 可观察A不满秩，但b在A的值空间中，这类方程用无数解。矩阵除法能正确求得这类方程的特解。逆阵法不能求得这类方程的特解。注意特解和齐次解目的A不满秩，但b在A的值空间中，这类方程的求解过程。rref 用于方程求解。矩阵除法能正确求得这类方程的特解。逆阵法不能求得这类方程的特解。解的验证方法。齐次解的获取。全解的获得。解答（1）借助增广矩阵用指令rref求解A=magic(4); %产生3阶魔方阵b=ones(4,1); %全1列向量R,C=rref(A,b) %求解，并判断解的唯一性 R

28、 = 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 关于以上结果的说明：R阶梯阵提供的信息前4列是原A阵经消元变换后的阶梯阵；而第5列是原b向量经相同变换后的结果。R的前三列为“基”，说明原A阵秩为3；而第4列的前三个元素，表示原A阵的第4列由其前三列线性组合而成时的加权系数，即方程的一个解。R的第5列表明：b可由原A阵的前三列线性表出；b给出了方程的一个解；由于原A阵“缺秩”，所以方程的确解不唯一。C数组提供的信息该数组中的三个元素表示变换取原A阵的第1，2，3列为基。该数组的元素总数就是“原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在设计“除法”程序时，针对“b在A值空间中”的情况，就是用rref求解的。）（3）