《matlab课后习题答案第四章.pdf》由会员分享，可在线阅读，更多相关《matlab课后习题答案第四章.pdf（25页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、精品-第4章 数值运算习题 4 及解答1 1根据题给的模拟实际测量数据的一组根据题给的模拟实际测量数据的一组t和和y(t)试用数值差分试用数值差分diffdiff或或数值梯度数值梯度 gradientgradient 指令计算指令计算y(t)，然后把，然后把y(t)和和y(t)曲线绘制在同曲线绘制在同一张图上，观察数值求导的后果。（模拟数据从一张图上，观察数值求导的后果。（模拟数据从prob_data401.matprob_data401.mat获得）获得）目的强调：要非常慎用数值导数计算。练习 mat 数据文件中数据的获取。实验数据求导的后果把两条曲线绘制在同一图上的一种方法。解答（1）从数
2、据文件获得数据的指令假如 prob_data401.mat 文件在当前目录或搜索路径上clearclearload prob_data401.matload prob_data401.mat（2）用 diff 求导的指令dt=t(2)-t(1);dt=t(2)-t(1);yc=diff(y)/dt;yc=diff(y)/dt;%注意注意 ycyc 的长度将比的长度将比 y y 短短 1 1plot(t,y,b,t(2:end),yc,r)plot(t,y,b,t(2:end),yc,r)grid ongrid on1.510.50-0.5-1-1.5-201234567-精品精品-（3）用 g
3、radent 求导的指令（图形与上相似）dt=t(2)-t(1);dt=t(2)-t(1);yc=gradient(y)/dt;yc=gradient(y)/dt;plot(t,y,b,t,yc,r)plot(t,y,b,t,yc,r)grid ongrid on说明不到万不得已，不要进行数值求导。假若一定要计算数值导数，自变量增量dt 要取得比原有数据相对误差高1、2 个量级以上。求导会使数据中原有的噪声放大。2 2采用数值计算方法，采用数值计算方法，画出画出y(x)计算计算y(4.5)。x0sin tdt在在0,10区间曲线，并区间曲线，并t提示指定区间内的积分函数可用cumtrapz 指
4、令给出。y(4.5)在计算要求不太高的地方可用find 指令算得。目的指定区间内的积分函数的数值计算法和cumtrapz 指令。find 指令的应用。解答dt=1e-4;dt=1e-4;t=0:dt:10;t=0:dt:10;t=t+(t=0)*eps;t=t+(t=0)*eps;f=sin(t)./t;f=sin(t)./t;s=cumtrapz(f)*dt;s=cumtrapz(f)*dt;plot(t,s,LineWidth,3)plot(t,s,LineWidth,3)ii=find(t=4.5);ii=find(t=4.5);s45=s(ii)s45=s(ii)s45=1.6541-
5、精品精品-21.81.61.41.210.80.60.40.2002468103 3求函数求函数f(x)e尝试复算。尝试复算。提示数值积分均可尝试。符号积分的局限性。目的符号积分的局限性。解答dx=pi/2000;dx=pi/2000;x=0:dx:pi;x=0:dx:pi;s=trapz(exp(sin(x).3)*dxs=trapz(exp(sin(x).3)*dxs=5.1370sin3x的数值积分的数值积分s 0f(x)dx，并请采用符号计算，并请采用符号计算符号复算的尝试syms xsyms xf=exp(sin(x)3);f=exp(sin(x)3);ss=int(f,x,0,pi
6、)ss=int(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 4用用 quadquad 求取求取1.7 5e xsin x dx的数值积分，的数值积分，并保证积分的绝对精度并保证积分的绝对精度-精品精品-为为109。目的quadl，精度可控，计算较快。近似积分指令 trapz 获得高精度积分的内存和时间代价较高。解答%精度可控的数值积分精度可控的数值积分fx=(x)exp(-abs(x).*abs(sin(x);fx=(x)exp(-abs(x
7、).*abs(sin(x);format longformat longsq=quadl(fx,-10*pi,1.7*pi,1e-7)sq=quadl(fx,-10*pi,1.7*pi,1e-7)sq=1.08784993815498%近似积分算法近似积分算法x=linspace(-10*pi,1.7*pi,1e7);x=linspace(-10*pi,1.7*pi,1e7);dx=x(2)-x(1);dx=x(2)-x(1);st=trapz(exp(-abs(x).*abs(sin(x)*dxst=trapz(exp(-abs(x).*abs(sin(x)*dxst=1.087849499
8、73430%符号积分算法符号积分算法y=exp(-abs(x)*abs(sin(x)y=exp(-abs(x)*abs(sin(x)si=vpa(int(y,-10*pi,1.7*pi),16)si=vpa(int(y,-10*pi,1.7*pi),16)y=exp(-abs(x)*abs(sin(x)si=1.0878494994129115 5求函数求函数f(t)(sin5t)值点。值点。2e0.06t1.5tcos2t 1.8t 0.5在区间在区间5,5中的最小中的最小2目的理解极值概念的邻域性。如何求最小值。学习运用作图法求极值或最小值。感受符号法的局限性。解答（1）采用 fminbn
9、d 找极小值点在指令窗中多次运行以下指令，观察在不同数目子区间分割下，进行的极小值搜索。然后从一系列极小值点中，确定最小值点。clearclearft=(t)sin(5*t).2.*exp(0.06*t.*t)+1.8*abs(t+0.5)-1.5*t.*cos(2*t);ft=(t)sin(5*t).2.*exp(0.06*t.*t)+1.8*abs(t+0.5)-1.5*t.*cos(2*t);disp(disp(计算中，把计算中，把-5,5-5,5 分成若干搜索子区间。分成若干搜索子区间。)N=input(N=input(请输入子区间数请输入子区间数 N N，注意使，注意使 N=1N=1
10、？);%);%该指令只能在指令窗中运行该指令只能在指令窗中运行t tt=linspace(-5,5,N+1);t=linspace(-5,5,N+1);f for k=1:Nor k=1:N-精品精品-tmin(k),fobj(k)=fminbnd(ft,tt(k),tt(k+1);tmin(k),fobj(k)=fminbnd(ft,tt(k),tt(k+1);e endnd fobj,ii=sort(fobj);%fobj,ii=sort(fobj);%将目标值由小到大排列将目标值由小到大排列t tmin=tmin(ii);min=tmin(ii);%使极小值点做与目标值相应的重新排列使极
11、小值点做与目标值相应的重新排列fobj,tminfobj,tmin（2）最后确定的最小值点在N1,2,10的不同分割下，经观察，最后确定出最小值点是-1.28498111480531-1.28498111480531相应目标值是-0.18604801006545-0.18604801006545（3）采用作图法近似确定最小值点（另一方法）（A）在指令窗中运行以下指令：clearclearft=(t)sin(5*t).2.*exp(0.06*t.*t)+1.8*abs(t+0.5)-1.5*t.*cos(2*t);ft=(t)sin(5*t).2.*exp(0.06*t.*t)+1.8*abs(
12、t+0.5)-1.5*t.*cos(2*t);t=-5:0.001:5;t=-5:0.001:5;ff=ft(t);ff=ft(t);plot(t,ff)plot(t,ff)grid on,shggrid on,shg（B）经观察后，把最小值附近邻域放到足够大，然后运行以下指令，那放大图形被推向前台，与此同时光标变为“十字线”，利用它点击极值点可得到最小值数据tmin2,fobj2=ginput(1)tmin2,fobj2=ginput(1)tmin2=-1.28500000993975fobj2=-0.18604799369136-精品精品-出现具有相同数值的刻度区域表明已达最小可分辨状态（
13、4）符号法求最小值的尝试syms tsyms tfts=sin(5*t)2*exp(0.06*t*t)-1.5*t*cos(2*t)+1.8*abs(t+0.5);fts=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);dfdt=diff(fts,t);%求导函数求导函数tmin=solve(dfdt,t)tmin=solve(dfdt,t)%求导函数的零点求导函数的零点fobj3=subs(fts,t,tmin)fobj3=subs(fts,t,tmin)%得到一个具体的极值点得到一个具体的极值点t
14、min=-.60100931947716486053884417850955e-2fobj3=.89909908144684551670208797723124说明最小值是对整个区间而言的，极小值是对邻域而言的。在一个区间中寻找最小值点，对不同子区间分割进行多次搜索是必要的。这样可以避免把极小值点误作为最小值点。最小值点是从一系列极小值点和边界点的比较中确定的。作图法求最小值点，很直观。假若绘图时，自变量步长取得足够小，那么所求得的最小值点有相当好的精度。符号法在本例中，只求出一个极值点。其余很多极值点无法秋初，更不可能得到最小值。-精品精品-d2y(t)6 6设设2 3dy(t)2y(t)1
15、,y(0)1,dy(0)0，用数值法和符号法求，用数值法和符号法求dtdtdty(t)t0.5。目的学习如何把高阶微分方程写成一阶微分方程组。ode45 解算器的导数函数如何采用匿名函数形式构成。如何从 ode45 一组数值解点，求指定自变量对应的函数值。解答（1）改写高阶微分方程为一阶微分方程组令y1(t)y(t),y2(t)dy(t)，于是据高阶微分方程可写出dtdy1(t)y2(t)dtdy(t)2 2y1(t)3y2(t)1dt（2）运行以下指令求 y(t)的数值解format longformat longts=0,1;ts=0,1;y0=1;0;y0=1;0;dydt=(t,y)y
16、(2);-2*y(1)+3*y(2)+1;dydt=(t,y)y(2);-2*y(1)+3*y(2)+1;%匿名函数写成的匿名函数写成的 ode45ode45 所需得导数函数所需得导数函数tt,yy=ode45(dydt,ts,y0);tt,yy=ode45(dydt,ts,y0);y_05=interp1(tt,yy(:,1),0.5,spline),%y_05=interp1(tt,yy(:,1),0.5,spline),%用一维插值求用一维插值求 y(0.5)y(0.5)y_05=0.78958020790127（3）符号法求解syms t;syms t;ys=dsolve(D2y-3*
17、Dy+2*y=1,y(0)=1,Dy(0)=0,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_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)function S=prob_DyDt(t,y)S=y(2);-2*y(1)+3*y(2)+1;S=y(2);-2*y(1)+
18、3*y(2)+1;7 7已知矩阵已知矩阵 A=magic(8)A=magic(8)，（，（1 1）求该矩阵的“值空间基阵”）求该矩阵的“值空间基阵”B B；（；（2 2）写出“写出“A A 的任何列可用基向量线性表出”的验证程序（提示：利的任何列可用基向量线性表出”的验证程序（提示：利用用 rrefrref 检验）。检验）。-精品精品-目的体验矩阵值空间的基向量组的不唯一性，但它们可以互为线性表出。利用 rref 检验两个矩阵能否互为表出。解答（1）A 的值空间的三组不同“基”A=magic(8);A=magic(8);%采用采用 8 8 阶魔方阵作为实验矩阵阶魔方阵作为实验矩阵R,ci=rr
19、ef(A);R,ci=rref(A);B1=A(:,ci)B1=A(:,ci)%直接从直接从 A A 中取基向量中取基向量B2=orth(A)B2=orth(A)%求求 A A 值空间的正交基值空间的正交基V,D=eig(A);V,D=eig(A);rv=sum(sum(abs(D)1000*eps);rv=sum(sum(abs(D)1000*eps);%非零特征值数就是矩阵的秩非零特征值数就是矩阵的秩B3=V(:,1:rv)B3=V(:,1:rv)%取取 A A 的非零特征值对应的特征向量作基的非零特征值对应的特征向量作基B1=64 2 3 9 55 54 17 47 46 40 26 2
20、7 32 34 35 41 23 22 49 15 14 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
21、.1906 -0.0451 0.3536 0.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=rref(B1,A)%若若 B1_AB1_A 矩阵的下矩阵的下 5 5 行全为行全为 0 0，%就表明就表明 A A 可以被可以被 B1B1 的的 3 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 -
22、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 0B2_A=rref(B2,A)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 -
23、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.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 0B3_A=rref(B3,A)B3_A=rref(B3,A)B3_A=Columns 1 through 7 1
24、.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 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
25、.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 倍数时，则矩阵的秩等于(n/2+2)。关于魔方阵的有关历史，请见第6.1.3 节。8 8已知由已知由 MATLABMATLAB 指令创建的矩阵指令创建的矩阵 A=gallery(5)A=gallery(5)，试对该矩阵进行，试对该矩阵进行特征值分解，并通过验算观察发生的现象。特征值分解，并通过验
26、算观察发生的现象。目的展示特征值分解可能存在的数值问题。condeig 是比较严谨的特征值分解指令。Jordan 分解的作用。解答（1）特征值分解A=gallery(5)A=gallery(5)V,D=eig(A);V,D=eig(A);diag(D)diag(D)%为紧凑地显示特征值而写为紧凑地显示特征值而写A=-9 11 -21 63 -252 70 -69 141 -421 1684 -575 575 -1149 3451 -13801 3891 -3891 7782 -23345 93365 1024 -1024 2048 -6144 24572ans=Columns 1 throug
27、h 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/VAE=V*D/Ver_AE=norm(A-AE,fro)/norm(A,fro)er_AE=norm(A-AE,fro)/norm(A,fro)%相对相对 F F 范数范数AE=1.0e+004*Columns 1 through 4 -0.0009+0.0000i 0.0011-0.0000i -0.0021+0.0000i 0.0063-0.0000i 0.0070-0.0000i -0
28、.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+0.0000i 0.1684-0.0000i -1.3800+0.0000i 9.3359-0.0001
29、i 2.4570-0.0000i-精品精品-er_AE=6.9310e-005（3）一个更严谨的特征值分解指令Vc,Dc,eigc=condeig(A)Vc,Dc,eigc=condeig(A)%eigc%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.0000i -0.1397 -0.1397+0.0000i -0
30、.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 0 0 -0.0054-0.0171i 0 0 0 0 0.0144+0.0104i 0 0 0 0 Colu
31、mn 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);VJ,DJ=jordan(A);%求出准确的广义特征值，使求出准确的广义特征值，使 A*VJ=VJ*DA*VJ=VJ*D 成立。成立。DJDJAJ=VJ*DJ/VJAJ=VJ*DJ/VJer_AJ=norm(A-AJ,fro)/norm(A,fro)er_AJ=norm(A-AJ,fro)/norm(A,fro)DJ=0 1 0 0 0 0 0 1 0 0 0 0 0 1
32、 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说明指令 condeig 的第 3 输出量 eigc 给出的是所谓的“矩阵特征值条件数”。当特征条件数与1eps相当
33、时，就意味着矩阵A 可能“退化”，即矩阵可能存在“代数重数”大于“几何重数”的特征值。此时，实施Jordan 分解更适宜。顺便指出：借助 condeig 算得的特征值条件数与cond 指令算得的矩阵条件数是两个不同概念。前者描述特征值的问题，后者描述矩阵逆的问题。本例矩阵 A 的特征值条件数很高，表明分解不可信。验算也表明，相对误差较大。当对矩阵 A 进行 Jordan 分解时，可以看到，A 具有 5 重根。当对 Jordan 分解进行验算时，相对误差很小。9 9求矩阵求矩阵AxAx b b的解，的解，A A 为为 3 3 阶魔方阵，阶魔方阵，b b 是是(31)的全的全 1 1 列向量。列向
34、量。提示了解 magic 指令rref 用于方程求解。矩阵除法和逆阵法解方程。目的满秩方阵求解的一般过程。rref 用于方程求解。矩阵除法和逆阵法解方程。解答A=magic(3);A=magic(3);%产生产生 3 3 阶魔方阵阶魔方阵b=ones(3,1);b=ones(3,1);%（3*13*1）全）全 1 1 列向量列向量R,C=rref(A,b)R,C=rref(A,b)%Gauss Jordan%Gauss Jordan 消去法解方程，同时判断解的唯一性消去法解方程，同时判断解的唯一性x=Abx=Ab%矩阵除解方程矩阵除解方程xx=inv(A)*bxx=inv(A)*b%逆阵法解方
35、程逆阵法解方程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 满秩，因此方程解唯一。在本例情况下，矩阵除、逆阵法、rref 法所得解一致。1010求矩阵求矩阵AxAx b b的解，的解，A A 为为 4 4 阶魔方阵，阶魔方阵，b b 是是(41)的全的全 1 1 列向量。列向量。提示
36、用 rref 可观察 A 不满秩，但 b 在 A 的值空间中，这类方程用无数解。矩阵除法能正确求得这类方程的特解。逆阵法不能求得这类方程的特解。注意特解和齐次解目的A 不满秩，但 b 在 A 的值空间中，这类方程的求解过程。rref 用于方程求解。矩阵除法能正确求得这类方程的特解。逆阵法不能求得这类方程的特解。解的验证方法。齐次解的获取。全解的获得。解答（1）借助增广矩阵用指令rref 求解A=magic(4);A=magic(4);%产生产生 3 3 阶魔方阵阶魔方阵b=ones(4,1);b=ones(4,1);%全全 1 1 列向量列向量R,C=rref(A,b)R,C=rref(A,b
37、)%求解，并判断解的唯一性求解，并判断解的唯一性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关于以上结果的说明：R 阶梯阵提供的信息前 4 列是原 A 阵经消元变换后的阶梯阵；而第 5 列是原 b 向量经相同变换后的结果。-精品精品-R 的前三列为“基”，说明原A 阵秩为 3；而第4 列的前三个元素，表示原A 阵的第 4 列由其前三列线性组合而成时的加权系数，即方程的一个解。R 的第 5 列表明：b 可由原 A 阵的前三列线性表出；b 给出了方程的一
38、个解；由于原 A 阵“缺秩”，所以方程的确解不唯一。C 数组提供的信息该数组中的三个元素表示变换取原A 阵的第 1，2，3 列为基。该数组的元素总数就是“原A 阵的秩”（2）矩阵除求得的解x=Abx=AbWarning: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 求
39、解的。）（3）逆阵法的解xx=inv(A)*bxx=inv(A)*bWarning: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)b_rref=A(:,C)*R(C,5)b_d=A*xb_d=A*xb_inv=A*xxb_inv=A*xxb_rref=1 1 1 1b_d=1 1 1 1b_inv=0.7344 1.5469%验算验算 r
40、refrref 的解的解%验算矩阵除的解验算矩阵除的解%验算逆阵法的解验算逆阵法的解-精品精品-1.1719 1.8594显然，在本例中，逆阵法的解是错误的。原因是：A 不满秩，A 的逆阵在理论上不存在。这里所给出的逆阵是不可信的。（4）求齐次解xg=null(A)xg=null(A)xg=0.2236 0.6708 -0.6708 -0.2236%Ax=0%Ax=0 的齐次解的齐次解（5）方程的全解齐次解的任何倍与特解之和就构成方程的全解。下面通过一组随机系数验证。rng defaultrng default%为本书结果可被读者核对而设，并非必要。为本书结果可被读者核对而设，并非必要。f=r
41、andn(1,6)f=randn(1,6)%6%6 个随机系数个随机系数xx=repmat(x,1,6)+xg*fxx=repmat(x,1,6)+xg*f%产生产生 6 6 个不同的特解个不同的特解A*xxA*xx%所得结果的每列都应该是全所得结果的每列都应该是全 1 1，即等于，即等于 b.b.f=0.5377 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.289
42、0 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 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解答（在用除法和逆阵法求解时出现）警告信息中RCOND=1.306145e-017 是矩阵 A 的估计条
43、件倒数。该数愈接近0，A 就愈“病态”；该数接近1 时，A 就愈“良态”。该条件数由 rcond(A)给出。注意：rcond 条件倒数与 cond 条件数的算法不同。12。1111求矩阵求矩阵AxAx b b的解，的解，A A 为为 4 4 阶魔方阵，阶魔方阵，b b 3 4提示由 rref 可以看出 A 不满秩，b 不在 A 的值空间中，方程没有准确解。但可求最小二乘近似解。-精品精品-目的A 不满秩，b 不在 A 的值空间中，方程没有准确解。解答（1）借助增广矩阵用指令rref 求解A=magic(4);A=magic(4);b=(1:4);b=(1:4);R,C=rref(A,b)R,C
44、=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%产生产生 3 3 阶魔方阵阶魔方阵%求解，并判断解的唯一性求解，并判断解的唯一性（2）用伪逆求最小二乘近似解（超出范围，仅供参考。）x=pinv(A)*bx=pinv(A)*bb_pinv=A*xb_pinv=A*xx=0.0235 0.1235 0.1235 0.0235b_pinv=1.3000 2.9000 2.1000 3.7000%非准确解非准确解%验算验算解答C 表明，A 的秩为 3，A 不满秩；R 第 5 列第 4 元素非零，说明 b 不在 A 的值空间中，
45、因此方程没有准确解，但可以求最小二乘近似解。1212求求0.5t 10e0.2tsinsint 0的实数解。的实数解。提示在适当范围内，作图观察一元复杂函数的形态：观察解的存在性；解的唯一性。进而，借助图形法求近似解。匿名函数的使用方法。fzero 指令的用法。目的作图法求一元复杂函数解上的作用：观察解的存在性；解的唯一性；得近似解。匿名函数的使用方法。fzero 指令的用法。解答（1）作图观察函数并求近似解-精品精品-t=-1:0.001:5;t=-1:0.001:5;y=(t)-0.5+t-10*exp(-0.2*t).*abs(sin(sin(t);y=(t)-0.5+t-10*exp(
46、-0.2*t).*abs(sin(sin(t);plot(t,y(t)plot(t,y(t)%利用匿名函数求利用匿名函数求 y y 函数值函数值grid on,shggrid on,shgtt1,yy1=ginput(1)tt1,yy1=ginput(1)%从图形获得近似解从图形获得近似解tt1=2.7370yy1=0.0097420-2-4-6-8-10-12-1012345（2）进一步利用 fzero 求精确解t,yt=fzero(y,tt1)t,yt=fzero(y,tt1)t=2.7341yt=2.2204e-015说明假如在从图上获取数据前，先把零点附近图形放大，可以得到精度更高的近
47、似解。sin(x y)01313求解二元函数方程组求解二元函数方程组的解。的解。cos(x y)0目的solve 指令的用法。解答（1）符号法只能得到两组解S=solve(sin(x-y),cos(x+y),x,y);S=solve(sin(x-y),cos(x+y),x,y);X=S.x,Y=S.yX=S.x,Y=S.yX=1/4*pi-1/4*pi-精品精品-Y=1/4*pi-1/4*pi（2）数值法可以看到许多解，并从中找到规律x=-6:0.1:6;y=x;X,Y=meshgrid(x,y);x=-6:0.1:6;y=x;X,Y=meshgrid(x,y);Z1=sin(X-Y);Z1=
48、sin(X-Y);Z2=cos(X+Y);Z2=cos(X+Y);contour(X,Y,Z1,3)contour(X,Y,Z1,3)%在在 Z1Z1 的取值范围内取的取值范围内取 3 3 个等分点作为等位线取值。个等分点作为等位线取值。%由于由于 Z1Z1 关于关于 z1=0z1=0 对称，所以中间线，即对称，所以中间线，即 Z1=0Z1=0 的等位线。的等位线。axis squareaxis squarehold onhold oncontour(X,Y,Z2,3)contour(X,Y,Z2,3)%注意：所有绿线交点都是解注意：所有绿线交点都是解hold offhold offgrid
49、on;shggrid on;shg1414假定某窑工艺瓷器的烧制成品合格率为假定某窑工艺瓷器的烧制成品合格率为目的0.1570.157，现该窑烧制，现该窑烧制 100100件瓷器，请画出合格产品数的概率分布曲线。件瓷器，请画出合格产品数的概率分布曲线。-精品精品-指令 binopdf 的用法。绘图指令 stem 的用法。解答N=100;p=0.157;N=100;p=0.157;k=0:N;k=0:N;pdf=binopdf(k,N,p);pdf=binopdf(k,N,p);stem(k,pdf)stem(k,pdf)grid ongrid onshgshg0.12%给定二项分布的特征参数给
50、定二项分布的特征参数%定义事件定义事件 A A 发生的次数数组发生的次数数组%算出各发生次数的概率算出各发生次数的概率0.10.080.060.040.0200204060801001515试产生均值为试产生均值为 4 4，标准差为标准差为 2 2 的的(10000 1)的正态分布随机数组的正态分布随机数组 a a，分别用分别用 histhist 和和 histfithistfit 绘制该数组的频数直方图，绘制该数组的频数直方图，观察两张图形的差观察两张图形的差异。除异。除 histfithistfit 上的拟合红线外，你能使这两个指令汇出相同的频上的拟合红线外，你能使这两个指令汇出相同的频数