matlab课后习题答案(1-9章)

1、1 数字1.5e2，1.5e3 中的哪个与1500相同吗？1.5e32 请指出如下5个变量名中，哪些是合法的？abcd-2 xyz_3 3chan a变量 ABCDefgh 2、5是合法的。3 在MATLAB环境中，比1大的最小数是多少？ 1+eps4 设 a = -8 , 运行以下三条指令，问运行结果相同吗？为什么？w1=a(2/3)w2=(a2)(1/3) w3=(a(1/3)2w1 = -2.0000 + 3.4641i ;w2 = 4.0000 ;w3 =-2.0000 + 3.4641i 5 指令clear, clf, clc各有什么用处？clear 清除工作空间中所有的变量。clf

2、 清除当前图形。clc 清除命令窗口中所有显示。第二章1 说出以下四条指令产生的结果各属于哪种数据类型，是“双精度”对象，还是“符号”符号对象？ 3/7+0.1双; sym(3/7+0.1)符; sym('3/7+0.1') 符; vpa(sym(3/7+0.1) 符;2 在不加专门指定的情况下，以下符号表达式中的哪一个变量被认为是自由符号变量.sym('sin(w*t)'),sym('a*exp(-X)'),sym('z*exp(j*th)')symvar(sym('sin(w*t)'),1) w a z3 （1

3、）试写出求三阶方程正实根的程序。注意：只要正实根，不要出现其他根。（2）试求二阶方程在时的根。（1）reset(symengine)syms x positivesolve(x3-44.5) ans =(2(2/3)*89(1/3)/2 （2）求五阶方程的实根syms a positive%注意：关于x的假设没有去除solve(x2-a*x+a2) Warning: Explicit solution could not be found.> In solve at 83ans = empty sym syms x clearsyms a positivesolve(x2-a*x+a2)

4、 ans = a/2 + (3(1/2)*a*i)/2 a/2 - (3(1/2)*a*i)/2 4 观察一个数（在此用记述）在以下四条不同指令作用下的异同。a =, b = sym( ), c = sym( ,'d ' ), d = sym( ' ' )在此， 分别代表具体数值 7/3 , pi/3 , pi*3(1/3) ；而异同通过vpa(abs(a-d) , vpa(abs(b-d) , vpa(abs(c-d)等来观察。l 理解准确符号数值的创建法。l 高精度误差的观察。（1）x=7/3x=7/3;a=x,b=sym(x),c=sym(x,'d

5、'),d=sym('7/3'), a = 2.3333b =7/3c =2.3333333333333334813630699500209d =7/3 v1=vpa(abs(a-d),v2=vpa(abs(b-d),v3=vpa(abs(c-d) v1 =0.0v2 =0.0v3 =0.00000000000000014802973661668756666666667788716 （2）x=pi/3x=pi/3;a=x,b=sym(x),c=sym(x,'d'),d=sym('pi/3'), a = 1.0472b =pi/3c =1.0

6、47197551196597631317786181171d =pi/3 v1=vpa(abs(a-d),v2=vpa(abs(b-d),v3=vpa(abs(c-d) v1 =0.0v2 =0.0v3 =0.00000000000000011483642827992216762806615818554 （3）x=pi*3(1/3)x=pi*3(1/3);a=x,b=sym(x),c=sym(x,'d'),d=sym('pi*3(1/3)') a = 4.5310b =1275352044764433/281474976710656c =4.5309606547

7、207899041040946030989d =pi*3(1/3) v1=vpa(abs(a-d),v2=vpa(abs(b-d),v3=vpa(abs(c-d) v1 =0.00000000000000026601114166290944374842393221638v2 =0.00000000000000026601114166290944374842393221638v3 =0.0000000000000002660111416629094726767991785515 5 求符号矩阵的行列式值和逆，所得结果应采用“子表达式置换”简洁化。l 理解subexpr指令。A=sym('

8、a11 a12 a13;a21 a22 a23;a31 a32 a33')DA=det(A)IA=inv(A);IAs,d=subexpr(IA,d) A = a11, a12, a13 a21, a22, a23 a31, a32, a33DA =a11*a22*a33 - a11*a23*a32 - a12*a21*a33 + a12*a23*a31 + a13*a21*a32 - a13*a22*a31IAs = d*(a22*a33 - a23*a32), -d*(a12*a33 - a13*a32), d*(a12*a23 - a13*a22) -d*(a21*a33 - a

9、23*a31), d*(a11*a33 - a13*a31), -d*(a11*a23 - a13*a21) d*(a21*a32 - a22*a31), -d*(a11*a32 - a12*a31), d*(a11*a22 - a12*a21)d =1/(a11*a22*a33 - a11*a23*a32 - a12*a21*a33 + a12*a23*a31 + a13*a21*a32 - a13*a22*a31) 6 求的符号解，并进而用该符号解求，的准确值。l symsum, subs的应用。l 从实例中，感受指令所给出的关于符号解的含义。syms x kf=x(k);Z1=symsu

10、m(f,k,0,inf)Z1 =piecewise(1 <= x, Inf, abs(x) < 1, -1/(x - 1) subs(Z1,x,sym('-1/3'),sym('1/pi'),sym('3') ans = 3/4, -1/(1/pi - 1), Inf 7 对于，求。（提示：理论结果为）l 符号变量的限定性定义的作用。syms k;x=sym('x','positive');f_k=2/(2*k+1)*(x-1)/(x+1)(2*k+1);s=simple(symsum(f_k,k,0,

11、inf) %结果与理论值lnx相符！ s =piecewise(abs(x - 1) < x + 1, log(x) 注意l 解答中，条件abs(x - 1) < x + 1意味着：n 约束一：x-1<x+1 2>0 此式总成立，说明“无约束”。n 情况二：-(x-1)<x+1x>0此为“约束”，满足题意。8 （1）通过符号计算求的导数。（2）然后根据此结果，求和。l diff, limit指令的应用。l 如何理解运行结果。syms ty=abs(sin(t)d=diff(y) %求dy/dtd0_=limit(d,t,0,'left') %

12、求dy/dt|t=0-dpi_2=limit(d,t,pi/2) %求dy/dt|t=pi/2 y =abs(sin(t)d =sign(sin(t)*cos(t)d0_ =-1dpi_2 =0 9 求出的具有64位有效数字的积分值。l 符号积分的解析解和符号数值解。l 符号计算和数值计算的相互校验。（1）符号积分syms x clearsyms xy=exp(-abs(x)*abs(sin(x)si=vpa(int(y,-10*pi,1.7*pi),64) y =abs(sin(x)/exp(abs(x)si =1.08784949941290491316667187594817452089

13、5458535212845987519414166 （2）数值计算复验xx=-10*pi:pi/100:1.7*pi;sn=trapz(exp(-abs(xx).*abs(sin(xx)*pi/100 sn = 1.0877 10 计算二重积分。l 变上限二重积分的符号计算法。syms x yf=x2+y2;r=int(int(f,y,1,x2),x,1,2) r =1006/105 11 在区间，画出曲线，并计算。l 在符号计算中，经常遇到计算结果是特殊经典函数的情况。l 如何应用subs获得超过16位有效数字的符号数值结果。l 初步尝试ezplot指令的简便。（1）符号计算syms t x

14、;f=sin(t)/t;y=int(f,t,0,x)% 将得到一个特殊经典函数y5=subs(y,x,sym('4.5')ezplot(y,0,2*pi) y =sinint(x)y5 =1.6541404143792439835039224868515（2）数值计算复验tt=0:0.001:4.5;tt(1)=eps;yn=trapz(sin(tt)./tt)*0.001 yn = 1.6541 12 在的限制下，求的一般积分表达式，并计算的32位有效数字表达。l 一般符号解与高精度符号数值解。syms xsyms n positivef=sin(x)n;yn=int(f,x

15、,0,pi/2) y3s=vpa(subs(yn,n,sym('1/3')y3d=vpa(subs(yn,n,1/3) yn =beta(1/2, n/2 + 1/2)/2y3s =1.2935547796148952674767575125656y3d =1.2935547796148951782413405453553 13 求方程的解。l solve指令中，被解方程的正确书写，输出量的正确次序。eq1='x2+y2=1'eq2='x*y=2'x,y=solve(eq1,eq2,'x','y') x = (1/

16、2 + (15(1/2)*i)/2)(1/2)/2 - (1/2 + (15(1/2)*i)/2)(3/2)/2 - (1/2 + (15(1/2)*i)/2)(1/2)/2 + (1/2 + (15(1/2)*i)/2)(3/2)/2 (1/2 - (15(1/2)*i)/2)(1/2)/2 - (1/2 - (15(1/2)*i)/2)(3/2)/2 - (1/2 - (15(1/2)*i)/2)(1/2)/2 + (1/2 - (15(1/2)*i)/2)(3/2)/2y = (1/2 + (15(1/2)*i)/2)(1/2) -(1/2 + (15(1/2)*i)/2)(1/2) (

17、1/2 - (15(1/2)*i)/2)(1/2) -(1/2 - (15(1/2)*i)/2)(1/2) 14 求微分方程的通解，并绘制任意常数为1时解的图形。l 理解指令dsolve的正确使用。l 对dsolve输出结果的正确理解。l ezplot指令绘图时，如何进行线色控制。l 如何覆盖那些不能反映图形窗内容的图名。（1）求通解reset(symengine)clearsyms y xy=dsolve('0.2*y*Dy+0.25*x=0','x') y = 2(1/2)*(C3 - (5*x2)/8)(1/2) -2(1/2)*(C3 - (5*x2)/

18、8)(1/2) （2）根据所得通解中不定常数的符号写出“对其进行数值替代的指令”yy=subs(y,'C3',1) %将通解中的C3用1代替 yy = 2(1/2)*(1 - (5*x2)/8)(1/2) -2(1/2)*(1 - (5*x2)/8)(1/2) （3）观察通解中两个分解的平方是否相同yy(1)2=yy(2)2 ans = 1 （4）于是可考虑函数的平方关系syms Yfxy=Y2-yy(1)2 fxy =Y2 + (5*x2)/4 - 2 （5）根据平方关系式画完整曲线clfezplot(fxy,-2,2,-2,2)axis squaregrid on （6）假

19、如直接用“分解”画曲线，那么将是不完整的 ezplot(yy(1),hold oncc=get(gca,'Children');set(cc,'Color','r')ezplot(yy(2),axis(-2 2 -2 2)legend('y(1)','y(2)'),hold off;title(' ')%覆盖不完全的图名gridaxis square 15 求一阶微分方程的解。l 初值微分方程的符号解。l pretty指令的使用。x=dsolve('Dx=a*t2+b*t','

20、;x(0)=2','t')pretty(x)%比较易读的表达形式 x =(t2*(3*b + 2*a*t)/6 + 2 2 t (3 b + 2 a t) - + 2 6 16 求边值问题的解。（注意：相应的数值解法比较复杂）。l 边值微分方程的符号解。f,g=dsolve('Df=3*f+4*g','Dg=-4*f+3*g','f(0)=0,g(0)=1') f =sin(4*t)*exp(3*t)g =cos(4*t)*exp(3*t) （1） 数值数组及其运算习题3及解答6 要求在闭区间上产生具有10个等距采样点的一

21、维数组。试用两种不同的指令实现。第 1 章 数值计算中产生自变量采样点的两个常用指令的异同。%方法一 t1=linspace(0,2*pi,10)%方法二t2=0:2*pi/9:2*pi %要注意采样间距的选择，如这里的2*pi/9. t1 = Columns 1 through 7 0 0.6981 1.3963 2.0944 2.7925 3.4907 4.1888 Columns 8 through 10 4.8869 5.5851 6.2832t2 = Columns 1 through 7 0 0.6981 1.3963 2.0944 2.7925 3.4907 4.1888 Col

22、umns 8 through 10 4.8869 5.5851 6.2832 1 由指令rng('default'),A=rand(3,5)生成二维数组A，试求该数组中所有大于0.5的元素的位置，分别求出它们的“全下标”和“单下标”。第 1 章 数组下标的不同描述：全下标和单下标。第 1 章 sub2ind, int2str, disp的使用。第 1 章 随机发生器的状态控制：保证随机数的可复现性。rng('default')A=rand(3,5)ri,cj=find(A>0.5);id=sub2ind(size(A),ri,cj);ri=ri'c

23、j=cj'disp(' ')disp('大于0.5的元素的全下标')disp('行号 ',int2str(ri)disp('列号 ',int2str(cj)disp(' ')disp('大于0.5的元素的单下标')disp(id') A = 0.8147 0.9134 0.2785 0.9649 0.9572 0.9058 0.6324 0.5469 0.1576 0.4854 0.1270 0.0975 0.9575 0.9706 0.8003 大于0.5的元素的全下标行号 1 2

24、 1 2 2 3 1 3 1 3列号 1 1 2 2 3 3 4 4 5 5 大于0.5的元素的单下标 1 2 4 5 8 9 10 12 13 15 1 已知矩阵，运行指令B1=A.(0.5), B2=A(0.5), 可以观察到不同运算方法所得结果不同。（1）请分别写出根据B1, B2恢复原矩阵A的程序。（2）用指令检验所得的两个恢复矩阵是否相等。第 1 章 数组运算和矩阵运算的不同。第 1 章 如何判断两个双精度数组是否相等。第 1 章 norm指令的应用。A=1,2;3,4;B1=A.0.5B2=A0.5A1=B1.*B1;A2=B2*B2;norm(A1-A2,'fro'

25、;)% 求误差矩阵的F-范数，当接近eps量级时，就认为实际相等B1 = 1.0000 1.4142 1.7321 2.0000B2 = 0.5537 + 0.4644i 0.8070 - 0.2124i 1.2104 - 0.3186i 1.7641 + 0.1458ians = 8.4961e-016 1 在时间区间 0,10中，绘制曲线。要求分别采取“标量循环运算法”和“数组运算法”编写两段程序绘图。 第 1 章 加强理解数组运算的机理和应用。第 1 章 初步使用subplot, plot, xlabel, ylabel等指令绘图。%标量循环运算法t=linspace(0,10,200)

26、;N=length(t);y1=zeros(size(t);for k=1:Ny1(k)=1-exp(-0.5*t(k)*cos(2*t(k);endsubplot(1,2,1),plot(t,y1),xlabel('t'),ylabel('y1'),grid on%数组运算法y2=1-exp(-0.5*t).*cos(2*t);subplot(1,2,2),plot(t,y2),xlabel('t'),ylabel('y2'),grid on 1 先运行clear,format long,rng('default'

27、;),A=rand(3,3)，然后根据A写出两个矩阵：一个对角阵B，其相应元素由A的对角元素构成；另一个矩阵C，其对角元素全为0，而其余元素与对应的A阵元素相同。第 1 章 常用指令diag的使用场合。clear,format longrng('default')A=rand(3,3)B=diag(diag(A)C=A-B A = 0.814723686393179 0.913375856139019 0.278498218867048 0.905791937075619 0.632359246225410 0.546881519204984 0.126986816293506

28、 0.097540404999410 0.957506835434298B = 0.814723686393179 0 0 0 0.632359246225410 0 0 0 0.957506835434298C = 0 0.913375856139019 0.278498218867048 0.905791937075619 0 0.546881519204984 0.126986816293506 0.097540404999410 0 1 先运行指令x=-3*pi:pi/15:3*pi; y=x; X,Y=meshgrid(x,y); warning off; Z=sin(X).*sin

29、(Y)./X./Y; 产生矩阵Z。（1）请问矩阵Z中有多少个“非数”数据？（2）用指令surf(X,Y,Z); shading interp观察所绘的图形。（3）请写出绘制相应的“无裂缝”图形的全部指令。第 1 章 初步感受三维曲面的绘制方法。第 1 章 非数NaN的产生，非数的检测，和对图形的影响。第 1 章 sum的应用。第 1 章 eps如何克服“被零除”的尴尬。x=-3*pi:pi/15:3*pi;y=x;X,Y=meshgrid(x,y);warning offZ=sin(X).*sin(Y)./X./Y;NumOfNaN=sum(sum(isnan(Z)%计算“非数”数目subpl

30、ot(1,2,1),surf(X,Y,Z),shading interp,title('有缝图')%产生无缝图XX=X+(X=0)*eps;YY=Y+(Y=0)*eps;ZZ=sin(XX).*sin(YY)./XX./YY;subplot(1,2,2),surf(XX,YY,ZZ),shading interp,title('无缝图') NumOfNaN = 181 1 下面有一段程序，企图用来解决如下计算任务：有矩阵，当依次取10, 9, 8, 7, 6, 5, 4, 3, 2, 1时，计算矩阵“各列元素的和”，并把此求和结果存放为矩阵Sa的第k行。例如时，

31、A阵为，此时它各列元素 的和是一个行数组，并把它保存为Sa的第3行。问题：该段程序的计算结果对吗？假如计算结果不正确，请指出错误发生的根源，并改正之。第 1 章 正确理解sum的工作机理。第 1 章 reshape的应用。（1）企图用以下程序完成题目要求。for k=10:-1:1A=reshape(1:10*k,k,10);Sa(k,:)=sum(A);endSa Sa = 55 55 55 55 55 55 55 55 55 55 3 7 11 15 19 23 27 31 35 39 6 15 24 33 42 51 60 69 78 87 10 26 42 58 74 90 106 1

32、22 138 154 15 40 65 90 115 140 165 190 215 240 21 57 93 129 165 201 237 273 309 345 28 77 126 175 224 273 322 371 420 469 36 100 164 228 292 356 420 484 548 612 45 126 207 288 369 450 531 612 693 774 55 155 255 355 455 555 655 755 855 955 （2）正确性分析除k=1外，计算所得Sa所有行的结果都正确。但k=1时，Sa的第1行应该与相同。上述程序的错误是对sum理

33、解不正确。sum对二维数组，求和按列施行；而对一维数组，不管行数组或列数组，总是求那数组所有元素的和。正确的程序应该写成for k=10:-1:1A=reshape(1:10*k,k,10);Sa(k,:)=sum(A);if k=1Sa(k,:)=A;endendSa Sa = 1 2 3 4 5 6 7 8 9 10 3 7 11 15 19 23 27 31 35 39 6 15 24 33 42 51 60 69 78 87 10 26 42 58 74 90 106 122 138 154 15 40 65 90 115 140 165 190 215 240 21 57 93 12

34、9 165 201 237 273 309 345 28 77 126 175 224 273 322 371 420 469 36 100 164 228 292 356 420 484 548 612 45 126 207 288 369 450 531 612 693 774 55 155 255 355 455 555 655 755 855 955 第 1 章 数值运算习题 4 及解答n 根据题给的模拟实际测量数据的一组和 试用数值差分diff或数值梯度gradient指令计算，然后把和曲线绘制在同一张图上，观察数值求导的后果。（模拟数据从prob_data401.mat获得）l 强

35、调：要非常慎用数值导数计算。l 练习mat数据文件中数据的获取。l 实验数据求导的后果l 把两条曲线绘制在同一图上的一种方法。（1）从数据文件获得数据的指令假如prob_data401.mat文件在当前目录或搜索路径上clearload prob_data401.mat （2）用diff求导的指令dt=t(2)-t(1);yc=diff(y)/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;plo

36、t(t,y,'b',t,yc,'r')grid on 说明l 不到万不得已，不要进行数值求导。l 假若一定要计算数值导数，自变量增量dt 要取得比原有数据相对误差高1、2个量级以上。l 求导会使数据中原有的噪声放大。n 采用数值计算方法，画出在区间曲线，并计算。提示l 指定区间内的积分函数可用cumtrapz指令给出。l 在计算要求不太高的地方可用find指令算得。l 指定区间内的积分函数的数值计算法和cumtrapz指令。l find指令的应用。dt=1e-4;t=0:dt:10;t=t+(t=0)*eps;f=sin(t)./t;s=cumtrapz(f)*

37、dt;plot(t,s,'LineWidth',3)ii=find(t=4.5);s45=s(ii) s45 = 1.6541 n 求函数的数值积分，并请采用符号计算尝试复算。提示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=int(f,x,0,pi) Warning: Explicit integral could not be found.> In at 58

38、ss =int(exp(sin(x)3),x = 0 . pi) n 用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 %近似积分算法x=linspace(-10*pi,1.7*pi,1e7);dx=x(2)-x(1);st=trapz(exp(-abs(x).*abs(sin(x)*dx

39、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 n 求函数在区间中的最小值点。l 理解极值概念的邻域性。l 如何求最小值。l 学习运用作图法求极值或最小值。l 感受符号法的局限性。（1）采用fminbnd找极小值点在指令窗中多次运行以下指令，观察在不同数目子区间分割下，进行的极小值搜索。然后从一系列极小值点中，确定最小值点。clearft=(t)sin(5*t)

40、.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);%将目标值由小到大排列tmin=tmin(ii);%使极小值点做与目标值相应的重新排列fobj,tmin（2）最后确定的最小值

41、点在的不同分割下，经观察，最后确定出最小值点是 -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）经观察后，把最小值附近邻域放到足够大，然后运行以下指令，那放大图形被推向前台，与此同时光标变为“十字线”，利用它点击极值点可得到最小值数据tmin2,fobj2=gin

42、put(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 =-.60100931947716486053884417850955e-2fobj3 =.89909908

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

44、值解点，求指定自变量对应的函数值。（1）改写高阶微分方程为一阶微分方程组令，于是据高阶微分方程可写出（2）运行以下指令求y(t)的数值解format longts=0,1;y0=1;0;dydt=(t,y)y(2);-2*y(1)+3*y(2)+1;%<4>%匿名函数写成的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

45、*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 第<4>条指令中的导数函数也可采用M函数文件表达，具体如下。function S=prob_DyDt(t,y)S=y(2);-2*y(1)+3*y(2)+1; n 求矩阵的解，A为3阶魔方阵，b是的全1列向量。提示l 了解magic指令l rref 用于方程求解。l 矩

46、阵除法和逆阵法解方程。l 满秩方阵求解的一般过程。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指令通过对增广矩阵

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

48、助增广矩阵用指令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 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.