admin 管理员组文章数量: 887021
2023年12月23日发(作者：创建网站的基本流程的教案)
习题二
1. 如何理解“矩阵是MATLAB最基本的数据对象”？
答：因为向量可以看成是仅有一行或一列的矩阵，单个数据（标量）可以看成是仅含一个元素的矩阵，故向量和单个数据都可以作为矩阵的特例来处理。
2. 设A和B是两个同维同大小的矩阵，问：
(1) A*B和A.*B的值是否相等？
答：不相等。
(2) A./B和B.A的值是否相等？
答：相等。
(3) A/B和BA的值是否相等？
答：不相等。
(4) A/B和BA所代表的数学含义是什么？
答：A/B等效于B的逆右乘A矩阵，即A*inv(B)，而BA等效于B矩阵的逆左乘A矩阵，即inv(B)*A。
3. 写出完成下列操作的命令。
(1) 将矩阵A第2~5行中第1, 3, 5列元素赋给矩阵B。
答：B=A(2:5,1:2:5); 或B=A(2:5,[1 3 5])
(2) 删除矩阵A的第7号元素。
答：A(7)=[]
(3) 将矩阵A的每个元素值加30。
答：A=A+30;
(4) 求矩阵A的大小和维数。
答：size(A);
ndims(A);
(5) 将向量 t 的0元素用机器零来代替。
答：t(find(t==0))=eps;
(6) 将含有12个元素的向量 x 转换成矩阵。
答：reshape(x,3,4);
(7) 求一个字符串的ASCII码。
答：abs(‘123’); 或double(‘123’);
因此，矩阵是MATLAB最基本、最重要的数据对象。
(8) 求一个ASCII码所对应的字符。
答：char(49);
4. 下列命令执行后，L1、L2、L3、L4的值分别是多少？
A=1:9;B=10-A;...
L1=A==B;
L2=A<=5;
L3=A>3&A<7;
L4=find(A>3&A<7);
答：L1的值为[0, 0, 0, 0, 1, 0, 0, 0, 0]
L2的值为[1, 1, 1, 1, 1, 0, 0, 0, 0]
L3的值为[0, 0, 0, 1, 1, 1, 0, 0, 0]
L4的值为[4, 5, 6]
5. 已知
100.7780234145655
A32503269.54543.14
完成下列操作：
(1) 取出A的前3行构成矩阵B，前两列构成矩阵C，右下角32子矩阵构成矩阵D，B与C的乘积构成矩阵E。
答：B=A(1:3,:);
C=A(:,1:2);
D=A(2:4,3:4);
E=B*C;
(2) 分别求E
011111，11，01，00答：E 0111110010 00 find(A>=10&A<25)=[1; 5]。 6. 当A=[34, NaN, Inf, -Inf, -pi, eps, 0]时，分析下列函数的执行结果：all(A)，any(A)，isnan(A)，isinf(A)，isfinite(A)。 答：all(A)的值为0 any(A) 的值为1 isnan(A) 的值为[ 0, 1, 0, 0, 0, 0, 0] isinf(A) 的值为[ 0, 0, 1, 1, 0, 0, 0] isfinite(A) 的值为[1, 0, 0, 0, 1, 1, 1] 7. 用结构体矩阵来存储5名学生的基本情况数据，每名学生的数据包括学号、姓名、专业和6门课程的成绩。 答：student(1).id='0001'; student(1).name='Tom'; student(1).major='computer'; student(1).grade=[89,78,67,90,86,85]; 8. 建立单元矩阵B并回答有关问题。 B{1,1}=1; B{1,2}='Brenden'; B{2,1}=reshape(1:9,3,3); B{2,2}={12,34,2;54,21,3;4,23,67}; (1) size(B)和ndims(B)的值分别是多少？ 答：size(B) 的值为2, 2。 ndims(B) 的值为2。 (2) B(2)和B(4)的值分别是多少？ 147，B(4)= 258答：B(2)=369[12][34][2][54][21][3] [4][23][67](3) B(3)=[]和B{3}=[]执行后，B的值分别是多少？ 答：当执行B(3)=[]后， B={1, [1, 4, 7; 2, 5, 8; 3, 6, 9], {12, 34, 2; 54, 21, 3; 4, 23, 67}} 当执行B{3}=[]后， B={1,[]; [1, 4, 7; 2, 5, 8; 3, 6, 9], {12, 34, 2; 54, 21, 3; 4, 23, 67}} 习题三 1. 写出完成下列操作的命令。 (1) 建立3阶单位矩阵A。 答：A=eye(3); (2) 建立5×6随机矩阵A，其元素为[100,200]范围内的随机整数。 答：round(100+(200-100)*rand(5,6)); (3) 产生均值为1，方差为0.2的500个正态分布的随机数。 答：1+sqrt(0.2)*randn(5,100); (4) 产生和A同样大小的幺矩阵。 答：ones(size(A)); (5) 将矩阵A对角线的元素加30。 答：A+eye(size(A))*30; (6) 从矩阵A提取主对角线元素，并以这些元素构成对角阵B。 答：B=diag(diag(A)); 2. 使用函数，实现方阵左旋90o或右旋90o的功能。例如，原矩阵为A，A左旋后得到B，右旋后得到C。 14710 A2581136912101112789 B456123321654 C987121110 答： B=rot90(A); C=rot90(A,-1); 3. 建立一个方阵A，求A的逆矩阵和A的行列式的值，并验证A与A-1是互逆的。 答： A=rand(3)*10; B=inv(A); C=det(A); 先计算B*A，再计算A*B，由计算可知B*A=A*B，即A·A-1= A-1·A是互逆。 4. 求下面线性方程组的解。 4x12x2x323x1x22x310 12x3x821 答： A=[4,2,-1;3,-1,2;12,3,0]; b=[2;10;8]; x=inv(A)*b 6.0000 26.6667方程组的解为x= EMBED 4 27.333332 295. 求下列矩阵的主对角线元素、上三角阵、下三角阵、秩、范数、条件数和迹。 112514(1) A30511150 0.43432(2) B 8.9421答： (1) 取主对角线元素： diag(A); 上三角阵： triu(A)； 下三角阵： tril(A); 秩： rank(A); 范数： norm(A,1); 或 norm(A); 或 norm(A,inf); 条件数： cond(A,1); 或 cond(A,2); 或 cond(A,inf) 迹： trace(A); 6. 求矩阵A的特征值和相应的特征向量。 (2)【请参考(1)】。 10.51 A110.2520.50.25 答： [V,D]=eig(A); 习题四 1. 从键盘输入一个4位整数，按如下规则加密后输出。加密规则：每位数字都加上7，然后用和除以10的余数取代该数字；再把第一位与第三位交换，第二位与第四位交换。 答： 2. 分别用if语句和switch语句实现以下计算，其中a、b、c的值从键盘输入。 a=input('请输入4位整数：'); A=[a/1000,a/100,a/10,a]; A=fix(rem(A,10)); A=rem(A+7,10); b=A(3)*1000+A(4)*100+A(1)*10+A(2); disp(['加密后的值为：',num2str(b)]); ax2bxc, 0.5x1.5yasincbx, 1.5x3.5 clnb, 3.5x5.5x答：(1) 用if语句实现计算： end y=a*((sin(b))^c)+x; y=a*x^2+b*x+c; a=input('请输入a的值：'); b=input('请输入b的值：'); c=input('请输入c的值：'); x=input('请输入x的值：'); end end disp(['y=',num2str(y)]); (2) 用switch语句实现计算： a=input('请输入a的值：'); b=input('请输入b的值：'); c=input('请输入c的值：'); x=input('请输入x的值：'); switch fix(x/0.5) y=log(abs(b+c/x)); case {1,2} y=a*x^2+b*x+c; case num2cell(3:6) y=a*((sin(b))^c)+x; case num2cell(7:10) y=log(abs(b+c/x)); end disp(['y=',num2str(y)]); 3. 产生20个两位随机整数，输出其中小于平均值的偶数。 答： A=fix(10+89*rand(1,20)); sum=0; for i=1:20 sum=sum+A(i); end B=A(find(A<(sum/20))); C=B(find(rem(B,2)==0)); disp(C); 4. 输入20个数，求其中最大数和最小数。要求分别用循环结构和调用MATLAB的max函数、min函数来实现。 答： (1) 用循环结构实现： v_max=0; v_min=0; for i=1:20 x=input(['请输入第', num2str(i), '数：']); if x> v_max v_max=x; end; if x< v_min end disp(['最大数为：', num2str(v_max)]); disp(['最小数为：', num2str(v_min)]); (2) 用max函数、min函数实现： for i=1:5 end disp(['最大数为：', num2str(max(A))]); disp(['最小数为：', num2str(min(A))]); 5. 已知：s122223数求s的值。 答： (1) 用循环结构实现： s=0; for i=0:63 s=s+2^i; end s (2) 调用sum函数实现： s=0:63; 263，分别用循环结构和调用MATLAB的sum函 v_min=x; end; A(i)=input(['请输入第', num2str(i), '数：']); s=2.^s; sum(s) 6. 当n分别取100、1000、10000时，求下列各式的值。 1111(1) 1(1)n1(ln2) 234n111(2) 1() 357411111(3) n() 4166443224466(4) 133557 (2n)(2n)(2n1)(2n1) 2要求分别用循环结构和向量运算（使用sum或prod函数）来实现。 (1) 用循环结构实现： sum=0; for k=1:100 end sum 使用sum函数： sum=sum+(-1)^(k+1)/k; 答： x=[]; for k=1:10000 end sum(x) (2) 用循环结构实现： sum=0; for k=1:100 end sum 使用sum函数： x=[]; sum=sum+(-1)^(k+1)/(2*k-1); x=[x, (-1)^(k+1)/k]; for k=1:100 end (3) 用循环结构实现： sum=0; for k=1:100 end sum 使用sum函数实现： x=[]; for k=1:100 end sum(x) (4) 用循环结构实现： t=1; for k=1:100 end t 使用prod函数实现： x=[]; for k=1:100 end prod(x) 7. 编写一个函数文件，求小于任意自然数n的斐波那契（Fibnacci）数列各项。斐波那契数列定义如下： x=[x, ((2*k)*(2*k))/((2*k-1)*(2*k+1))]; t=t*(((2*k)*(2*k))/((2*k-1)*(2*k+1))); x=[x, 1/(4^k)]; sum=sum+1/(4^k); sum(x) x=[x, (-1)^(k+1)/(2*k-1)]; f11, n1f21, n2 fff, n2n1n2n答： function x=fibnacci(n) for i=1:n if i<=2 x(i)=1; else x(i)=x(i-1)+x(i-2); end end 8. 编写一个函数文件，用于求两个矩阵的乘积和点乘，然后在命令文件中调用该函数。 答： 函数文件myfnc.m： function [x, y]= myfnc(A, B) try x=A*B; catch end y=A.*B; 命令文件myexe.m： A=input('请输入矩阵A：'); B=input('请输入矩阵B：'); [x, y]=myfnc(A, B); if length(x)==0 end disp('矩阵A和矩阵B的点乘为：'); y display('两矩阵的维数不匹配，无法进行乘积运算！'); disp('矩阵A和矩阵B的乘积为：'); x else x=[]; 9. 先用函数的递归调用定义一个函数文件求im，然后调用该函数文件求i1nkkk。 2k1k1k110050101答： 10.写出下列程序的输出结果。 ① s=0; a=[12,13,14;15,16,17;18,19,20;21,22,23]; for k=a end s 答：执行结果为 s=108 x = 4 12 20 ② 执行后的结果为： for j=1:4 end if rem(k(j),2)~=0 end s=s+k(j); 函数文件myfnc.m： function sum=myfnc(n, m) if n<=1 end 1在命令窗口中调用myfnc.m文件，计算kk： k1k1k1k21005010sum=1; sum= myfnc (n-1, m)+n^m; else sum=myfnc(100, 1)+ myfnc(50, 2)+myfnc(10,-1) y= 2 4 6 第4章 数值运算 习题 4 及解答 1 根据题给的模拟实际测量数据的一组t和 y(t)试用数值差分diff或数值梯度gradient指令计算y(t)，然后把y(t)和y(t)曲线绘制在同一张图上，观察数值求导的后果。（模拟数据从prob_获得） 〖目的〗 强调：要非常慎用数值导数计算。 练习mat数据文件中数据的获取。 实验数据求导的后果 把两条曲线绘制在同一图上的一种方法。 〖解答〗 （1）从数据文件获得数据的指令 假如prob_文件在当前目录或搜索路径上 clear load prob_ （2）用diff求导的指令 dt=t(2)-t(1); yc=diff(y)/dt; %注意yc的长度将比y短1 plot(t,y,'b',t(2:end),yc,'r') grid on 1.510.50-0.5-1-1.5-201234567 （3）用gradent求导的指令（图形与上相似） dt=t(2)-t(1); yc=gradient(y)/dt; plot(t,y,'b',t,yc,'r') grid on 〖说明〗 不到万不得已，不要进行数值求导。 假若一定要计算数值导数，自变量增量dt 要取得比原有数据相对误差高1、2个量级以上。 求导会使数据中原有的噪声放大。 2 采用数值计算方法，画出y(x)计算y(4.5)。 x0sintdt在[0, 10]区间曲线，并t〖提示〗 指定区间内的积分函数可用cumtrapz指令给出。 y(4.5)在计算要求不太高的地方可用find指令算得。 〖目的〗 指定区间内的积分函数的数值计算法和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 21.81.61.41.210.80.60.40.200246810 3 求函数f(x)esin3x的数值积分s 0f(x)dx，并请采用符号计算 尝试复算。 〖提示〗 数值积分均可尝试。 符号积分的局限性。 〖目的〗 符号积分的局限性。 〖解答〗 dx=pi/2000; x=0:dx:pi; s=trapz(exp(sin(x).^3))*dx s = 5.1370 符号复算的尝试 syms x f=exp(sin(x)^3); ss=int(f,x,0,pi) Warning: Explicit integral could not be found. > In at 58 ss = int(exp(sin(x)^3),x = 0 .. pi) 4 用quad求取度为109。 1.75exsinxdx的数值积分，并保证积分的绝对精〖目的〗 quadl，精度可控，计算较快。 近似积分指令trapz获得高精度积分的内存和时间代价较高。 〖解答〗 %精度可控的数值积分 fx=@(x)exp(-abs(x)).*abs(sin(x)); format long sq=quadl(fx,-10*pi,1.7*pi,1e-7) sq = 1.98 %近似积分算法 x=linspace(-10*pi,1.7*pi,1e7); dx=x(2)-x(1); st=trapz(exp(-abs(x)).*abs(sin(x)))*dx st = 1.30 %符号积分算法 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.911 5 求函数f(t)(sin5t)值点。 2e0.06t1.5tcos2t1.8t0.5在区间[5,5]中的最小2〖目的〗 理解极值概念的邻域性。 如何求最小值。 学习运用作图法求极值或最小值。 感受符号法的局限性。 〖解答〗 （1）采用fminbnd找极小值点 在指令窗中多次运行以下指令，观察在不同数目子区间分割下，进行的极小值搜索。然后从一系列极小值点中，确定最小值点。 clear ft=@(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); %使极小值点做与目标值相应的重新排列 fobj,tmin （2）最后确定的最小值点 在N1,2,,10的不同分割下，经观察，最后确定出 最小值点是 -1.28498111480531 相应目标值是 -0.186 （3）采用作图法近似确定最小值点（另一方法） （A）在指令窗中运行以下指令： clear ft=@(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]=ginput(1) tmin2 = -1.285 fobj2 = -0.186 出现具有相同数值的刻度区域表明已达最小可分辨状态 （4）符号法求最小值的尝试 syms t 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); %求导函数 tmin=solve(dfdt,t) %求导函数的零点 fobj3=subs(fts,t,tmin) %得到一个具体的极值点 tmin = -.66486955e-2 fobj3 = .89951674 〖说明〗 最小值是对整个区间而言的，极小值是对邻域而言的。 在一个区间中寻找最小值点，对不同子区间分割进行多次搜索是必要的。这样可以避免把极小值点误作为最小值点。最小值点是从一系列极小值点和边界点的比较中确定的。 作图法求最小值点，很直观。假若绘图时，自变量步长取得足够小，那么所求得的最小值点有相当好的精度。 符号法在本例中，只求出一个极值点。其余很多极值点无法秋初，更不可能得到最小值。 d2y(t)6 设23dy(t)2y(t)1,y(0)1,dy(0)0，用数值法和符号法求dtdtdt y(t)t0.5。 〖目的〗 学习如何把高阶微分方程写成一阶微分方程组。 ode45解算器的导数函数如何采用匿名函数形式构成。 如何从ode45一组数值解点，求指定自变量对应的函数值。 〖解答〗 （1）改写高阶微分方程为一阶微分方程组 令y1(t)y(t),y2(t)dy(t)，于是据高阶微分方程可写出 dtdy1(t)y2(t)dt dy(t)22y1(t)3y2(t)1dt （2）运行以下指令求y(t)的数值解 format long ts=[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*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 = .7895891685 〖说明〗 第<4>条指令中的导数函数也可采用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检验）。 〖目的〗 体验矩阵值空间的基向量组的不唯一性，但它们可以互为线性表出。 利用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 59 B2 = -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.3536 B3 = 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.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 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.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 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倍数时，则矩阵的秩等于 (n/2+2)。关于魔方阵的有关历史，请见第6.1.3节。 8 已知由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 1024 -1024 2048 -6144 24572 ans = 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/V er_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.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.0000i -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.0000i er_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.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.0000i Dc = 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 Column 5 0 0 0 0 0.0144 - 0.0104i eigc = 1.0e+011 * 5.2687 5.2313 5.2313 5.1725 5.1724 （4）对A采用Jordan分解并检验 [VJ,DJ]=jordan(A); %求出准确的广义特征值，使A*VJ=VJ*D成立。 DJ AJ=VJ*DJ/VJ er_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 0 AJ = 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.4572 er_AJ = 2.0500e-011 〖说明〗 指令condeig的第3输出量eigc给出的是所谓的“矩阵特征值条件数”。当特征条件数与1eps相当时，就意味着矩阵A可能“退化”，即矩阵可能存在“代数重数”大于“几何重数”的特征值。此时，实施Jordan分解更适宜。 顺便指出：借助condeig算得的特征值条件数与cond指令算得的矩阵条件数是两个不同概念。前者描述特征值的问题，后者描述矩阵逆的问题。 本例矩阵A的特征值条件数很高，表明分解不可信。验算也表明，相对误差较大。 当对矩阵A进行Jordan分解时，可以看到，A具有5重根。当对Jordan分解进行验算时，相对误差很小。 9 求矩阵Axb的解，A为3阶魔方阵，b是(31)的全1列向量。 〖提示〗 了解magic指令 rref 用于方程求解。 矩阵除法和逆阵法解方程。 〖目的〗 满秩方阵求解的一般过程。 rref 用于方程求解。 矩阵除法和逆阵法解方程。 〖解答〗 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.0667 C = 1 2 3 x = 0.0667 0.0667 0.0667 xx = 0.0667 0.0667 0.0667 〖说明〗 rref指令通过对增广矩阵进行消去法操作完成解方程。由分解得到的3根“坐标向量”和（或）C3指示的3根基向量，可见A3满秩，因此方程解唯一。 在本例情况下，矩阵除、逆阵法、rref法所得解一致。 10 求矩阵Axb的解，A为4阶魔方阵，b是(41)的全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 = 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 0 C = 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）逆阵法的解 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 1 b_d = 1 1 1 1 b_inv = 0.7344 1.5469 1.1719 1.8594 显然，在本例中，逆阵法的解是错误的。原因是：A不满秩，A的逆阵在理论上不存在。这里所给出的逆阵是不可信的。 （4）求齐次解 xg=null(A) xg = 0.2236 0.6708 -0.6708 -0.2236 %Ax=0的齐次解 （5）方程的全解 齐次解的任何倍与特解之和就构成方程的全解。下面通过一组随机系数验证。 rng default %为本书结果可被读者核对而设，并非必要。 f=randn(1,6) %6个随机系数 xx=repmat(x,1,6)+xg*f %产生6个不同的特解 A*xx %所得结果的每列都应该是全1，即等于b. f = 0.5377 1.8339 -2.2588 0.8622 0.3188 -1.3077 xx = 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.2924 ans = 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的估计条件倒数。该数愈接近0，A就愈“病态”；该数接近1时，A就愈“良态”。该条件数由rcond(A)给出。注意：rcond条件倒数与cond条件数的算法不同。 1211 求矩阵Axb的解，A为4阶魔方阵，b3。 4〖提示〗 由rref可以看出A不满秩，b不在A的值空间中，方程没有准确解。 但可求最小二乘近似解。 〖目的〗 A不满秩，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 1 C = 1 2 3 5 （2）用伪逆求最小二乘近似解（超出范围，仅供参考。） x=pinv(A)*b b_pinv=A*x x = 0.0235 0.1235 %非准确解 %验算 0.1235 0.0235 b_pinv = 1.3000 2.9000 2.1000 3.7000 〖解答〗 C表明，A的秩为3，A不满秩；R第5列第4元素非零，说明b不在A的值空间中，因此方程没有准确解，但可以求最小二乘近似解。 12 求0.5t10e0.2tsin[sint]0的实数解。 〖提示〗 在适当范围内，作图观察一元复杂函数的形态：观察解的存在性；解的唯一性。进而，借助图形法求近似解。 匿名函数的使用方法。 fzero指令的用法。 〖目的〗 作图法求一元复杂函数解上的作用：观察解的存在性；解的唯一性；得近似解。 匿名函数的使用方法。 fzero指令的用法。 〖解答〗 （1）作图观察函数并求近似解 t=-1:0.001:5; y=@(t)-0.5+t-10*exp(-0.2*t).*abs(sin(sin(t))); plot(t,y(t)) %利用匿名函数求y函数值 grid on,shg [tt1,yy1]=ginput(1) %从图形获得近似解 tt1 = 2.7370 yy1 = 0.0097 420-2-4-6-8-10-12-1012345 （2）进一步利用fzero求精确解 [t,yt]=fzero(y,tt1) t = 2.7341 yt = 2.2204e-015 〖说明〗 假如在从图上获取数据前，先把零点附近图形放大，可以得到精度更高的近似解。 sin(xy)013 求解二元函数方程组的解。 cos(xy)0〖目的〗 solve指令的用法。 〖解答〗 （1）符号法 只能得到两组解 S=solve('sin(x-y)','cos(x+y)','x','y'); X=S.x,Y=S.y X = [ 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); Z1=sin(X-Y); Z2=cos(X+Y); contour(X,Y,Z1,3) %在Z1的取值范围内取3个等分点作为等位线取值。 %由于Z1关于z1=0对称，所以中间线，即Z1=0的等位线。 axis square hold on contour(X,Y,Z2,3) %注意：所有绿线交点都是解 hold off grid on;shg 14 假定某窑工艺瓷器的烧制成品合格率为0.157，现该窑烧制100件瓷器，请画出合格产品数的概率分布曲线。 〖目的〗 指令binopdf的用法。 绘图指令stem的用法。 〖解答〗 N=100;p=0.157; k=0:N; pdf=binopdf(k,N,p); stem(k,pdf) grid on shg %给定二项分布的特征参数 %定义事件A发生的次数数组 %算出各发生次数的概率 0.120.10.080.060.040.100 15 试产生均值为4，标准差为2的(100001)的正态分布随机数组 a ， 分别用hist和histfit绘制该数组的频数直方图，观察两张图形的差异。除histfit上的拟合红线外，你能使这两个指令汇出相同的频数直方图吗？ 〖目的〗 指令normrnd的应用，及随机发生器的初始状态控制。 指令hist 与 histfit的异同。 〖解答〗 （1）绘制频数直方图 rng default mu=4;sigma=2; %为本例结果可被读者重现而设，并非必要。 %定义均值和标准差 %a=4+2*randn(10000,1); a=normrnd(mu,sigma,10000,1); %产生正态分布随机数组a <3> subplot(2,1,1),hist(a) %绘制a的频数统计直方图 subplot(2,1,2),histfit(a) 3-44-4-202468112 （2）对两个绘图指令中的直方条的数目设置相同的值 subplot(2,1,1),hist(a,50) %绘制a的频数统计直方图 subplot(2,1,2),histfit(a,50) 8-48-4-202468112 〖说明〗 第<3>条指令可用一下指令等价代替 a=4+2*randn(10000,1); 16 从数据文件prob_得到随机数组R，下面有一段求取随机数组全部数据最大值、均值和标准差的程序。 Mx=max(max(R)),Me=mean(mean(R)),St=std(std(R)), 试问该程序所得的结果都正确吗？假如不正确，请写出正确的程 序。 〖目的〗 mat数据文件中数据的获取。 max, mean, std的正确运用。 〖解答〗 （1）获取随机数组R 先要保证数据文件prob_在当前目录上或搜索路径上，然后运行以下指令。 load prob_data416 （2）运行题给程序 Mx=max(max(R)) Me=mean(mean(R)) S=std(std(R)) Mx = 0.9997 Me = 0.5052 S = 0.0142 （3）正确程序 Mx1=max(R(:)) Me1=mean(R(:)) S1=std(R(:)) %<7> Mx1 = 0.9997 Me1 = 0.5052 S1 = 0.2919 〖说明〗 两种方法得到的最大值和均的计算结果相同。求取R全部数据标准差的正确方法见第<7>条指令。 题给方法求出的S不是R全部数据的标准差，而是对于“各列数据样本标准差”的标准差。该S很小，说明各列求得的标准差相互之间离散度很小。 17 已知有理分式R(x)N(x)，其中N(x)(3x3x)(x30.5)，D(x)D(x)(x22x2)(5x32x21)。（1）求该分式的商多项式Q(x)和余多项式r(x)。（2）用程序验算D(x)Q(x)r(x)N(x)是否成立。 〖目的〗 指令conv, deconv的应用。 计算结果的演算。 〖解答〗 （1）求商及余 clear N=conv([3 0 1 0],[1 0 0 0.5]) D=conv([1 2 -2],[5 2 0 1]) [Q,r]=deconv(N,D) N = 3.0000 0 1.0000 1.5000 0 0.5000 0 D = 5 12 -6 -3 2 -2 Q = 0.6000 -1.4400 r = -0.0000 0.0000 21.8800 -5.3400 -5.5200 4.5800 -2.8800 （2）验算 NN1=conv(Q,D); m=length(r); NN1(end-m+1:end)=NN1(end-m+1:end)+r %关键指令 err=norm(N-NN1)/norm(N) %多项式系数向量的相对范数误差 NN1 = 3.0000 0 1.0000 1.5000 0 0.5000 0 err = 0 18 现有一组实验数据x, y，试求这组数据的拟合多项式。 〖目的〗 prob_数据文件中数据的获取。 多项式拟合指令polyfit和多项式求值指令polyval的应用。 〖解答〗 load prob_data418 %要注意搜索路径 P=polyfit(x,y,5) %采用5阶多项式拟合 xx=linspace(x(1),x(end),5*length(x));%或为应用需要，或为精细绘图，增加数据点 yy=polyval(P,xx); plot(x,y,'r.',xx,yy,'-b') P = -0.0039 0.0338 -0.0227 -0.4456 0.9590 -0.0364 0.60.40.20-0.2-0.4-0.6-0.8-1-1.2-1.4-101234 〖说明〗 拟合多项式的阶数要适当，不宜过高。比如，对于本例而言，10阶以上多项式也许救过高了。 19 已知系统冲激响应为h(n)=[0.05,0.24,0.40,0.24,0.15,-0.1,0.1] ，系 统输入u(n)由指令rng('default');u=2*(randn(1,100)>0.5)-1产生，该输入信号的起始作用时刻为0。试用直杆图（提示：用stem指令）画出分别显示该系统输入、输出信号的两张子图。 〖目的〗 序列卷积的计算。 如何序列的起点的正确记述。 非数在绘图中的应用。 〖解答〗 clear h=[0.05,0.24,0.40,0.24,0.15,-0.1,0.1]; rng('default');u=2*(randn(1,100)>0.5)-1; %以下采用非平凡区间法求输出序列 n1=0;n2=length(h)-1; n3=0;n4=length(u)-1; ny1=n1+n3;ny2=n2+n4; y=conv(u,h); %以下为绘图 ny=length(y); uu=[u,NaN*ones(1,length(h)-1)]; subplot(2,1,1),stem(0:ny-1,uu,'filled') title('Input u'),axis([0,ny,-1,1]) subplot(2,1,2),stem(0:ny-1,y,'filled') title('Output y'),axis([0,ny,-1,1]) Input u10.50-0.5-1%系统冲激响应 %输入序列 %使uu与y长度相同 <11> %使横坐标刻度从0开始 Output y10.50-0.5-1 〖说明〗 第<11>条指令的作用：使uu与y长度相同。那些增添点的值取NaN，可使绘图正确反映从第100时刻起系统没有输入信号的事实。 ============================ 以下是参考题============================= 20 求方程ax2bxc0的解。其中ac1,b100000000。 〖目的〗 符号计算没有中间计算误差。 MATLAB提供指令是数值可靠的。 数值计算中要尽量避免“数值相近大数间的相减”操作。 〖解答〗 （1）符号法求解 syms z aa=sym('1');cc=aa;bb=sym('-100000000'); z=solve(aa*z^2+bb*z+cc,z) zz=vpa(z) z = 50000000+2499999999999999^(1/2) 500099999999^(1/2) zz = 99999999.999999989999999999999999 .10001e-7 （2）利用多项式求根指令求解 format long a=1;c=1;b=-1e8; y=roots([a,b,c]) %算法可靠 %但由于两个根的大小相差很大很大，所以显示结果可能引起误读 y = 1.0e+008 * 1.00 0.00 y(1) ans = 9.999999999999999e+007 y(2) %单独显示可避免误读 ans = 1.000e-008 bb24ac（3）利用二次方程求根公式x计算 2ax1=(-b+sqrt(b*b-4*a*c))/2/a x2=(-b-sqrt(b*b-4*a*c))/2/a x1 = 100000000 x2 = 7.458e-009 %存在数值相近大数相减操作 （4）在利用求根公式求得x1之后，再利用x1x2xx2=c/a/x1 xx2 = 1.000e-008 c求x2 a%采用此法算得的第二个根是比较精确的
版权声明：本文标题：matlab课后习题答案 内容由网友自发贡献，该文观点仅代表作者本人， 转载请联系作者并注明出处：http://www.freenas.com.cn/jishu/1703276136h445168.html， 本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌抄袭侵权/违法违规的内容，一经查实，本站将立刻删除。
发表评论