MATLAB教程2023a习题解答1

本文格式为Word版，下载可任意编辑——MATLAB教程2023a习题解答1

???

MATLABR2023a课后习题答案全解第一章基础准备及入门

习题1及解答

?1.数字1.5e2，1.5e3中的哪个与1500一致吗？

〖解答〗1.5e3

?2.请指出如下5个变量名中，哪些是合法的？

abcd-2xyz_33chana变量ABCDefgh

〖解答〗

2、5是合法的。

?3.在MATLAB环境中，比1大的最小数是多少？

〖解答〗1+eps

?4.设a=-8,运行以下三条指令，问运行结果一致吗？为什么？

w1=a^(2/3)w2=(a^2)^(1/3)w3=(a^(1/3))^2〖解答〗

（1）不同。具体如下

w1=a^(2/3)

%仅求出主根

1

w2=(a^2)^(1/3)w3=(a^(1/3))^2w1=w2=4.0000w3=

%求出(-8)^2的主根%求出(-8)主根后再平方

-2.0000+3.4641i

-2.0000+3.4641i

（2）复数的多方根的，下面是求取全部方根的两种方法：（A）根据复数方根定义

a=-8;n=2;m=3;

ma=abs(a);aa=angle(a);fork=1:m

%m决定循环次数%计算各根的相角%计算各根

sa(k)=(aa+2*pi*(k-1))*n/m;

end

result=(ma^(2/3)).*exp(j*sa)result=

-2.0000+3.4641i4.0000-0.0000i-2.0000-3.4641i

（B）利用多项式r?a?0求根

p=[1,0,0,-a^2];r=roots(p)r=

-2.0000+3.4641i-2.0000-3.4641i4.0000

32

?5.指令clear,clf,clc各有什么用处？

〖解答〗clear清除工作空间中所有的变量。clf清除当前图形。clc清除命令窗口中所有显示。

?6.以下两种说法对吗？（1）“MATLAB进行数值的表达精度与

其指令窗中的数据显示精度一致。〞（2）MATLAB指令窗中显示的数值有效位数不超过7位。〞

2

〖解答〗

（1）否；（2）否。

?123?456?，?7.想要在MATLAB中产生二维数组S??下面哪些指令????789??能实现目的？

（1）S=[1,2,3;4,5,6;7,8;9]

（2）S=[123;456;789]

（3）S=[1，2，3；4，5，6；7，8，9]〖解答〗

前两种输入方法可以，后一种方法不行。

%整个指令在中文状态下输入

?8.试为例1.3-5编写一个解题用的M脚本文件？

〖解答〗

直接点击新文件图标，出现M文件编辑器窗口；在该M文件编辑器中，输入例1.3-5中的全部指令；并另存为p109.m，便得到所需的脚本文件。

第2章符号运算

习题2及解答

?/1说出以下四条指令产生的结果各属于哪种数据类型，是“双精

度〞对象，还是“符号〞符号对象？

3/7+0.1;sym(3/7+0.1);sym('3/7+0.1');vpa(sym(3/7+0.1))

〖目的〗

?不能从显示形式判断数据类型，而必需依靠class指令。〖解答〗

3

c1=3/7+0.1c2=sym(3/7+0.1)c3=sym('3/7+0.1')c4=vpa(sym(3/7+0.1))Cs1=class(c1)Cs2=class(c2)Cs3=class(c3)Cs4=class(c4)c1=0.5286c2=37/70c3=

0.52857142857142857142857142857143c4=

0.52857142857142857142857142857143Cs1=doubleCs2=symCs3=symCs4=sym

?/2在不加专门指定的状况下，以下符号表达式中的哪一个变量被

认为是自由符号变量.

sym('sin(w*t)'),sym('a*exp(-X)'),sym('z*exp(j*th)')

〖目的〗

?理解自由符号变量的确认规则。〖解答〗

symvar(sym('sin(w*t)'),1)ans=w

symvar(sym('a*exp(-X)'),1)ans=a

4

symvar(sym('z*exp(j*th)'),1)ans=z

?/3求以下两个方程的解

3（1）试写出求三阶方程x?44.5?0正实根的程序。注意：只要正

实根，不要出现其他根。

（2）试求二阶方程x2?ax?a2?0在a?0时的根。

〖目的〗

?体验变量限定假设的影响〖解答〗

（1）求三阶方程x?44.5?0正实根

reset(symengine)symsxpositivesolve(x^3-44.5)ans=

(2^(2/3)*89^(1/3))/2

%确保下面操作不受前面指令运作的影响

3（2）求五阶方程x?ax?a?0的实根

symsapositivesolve(x^2-a*x+a^2)>Insolveat83ans=

[emptysym]

symsxclearsymsapositivesolve(x^2-a*x+a^2)ans=

a/2+(3^(1/2)*a*i)/2a/2-(3^(1/2)*a*i)/2

%注意：关于x的假设没有去除

22Warning:Explicitsolutioncouldnotbefound.

5

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

sinint(x)1.81.61.41.210.80.60.40.202323x456%将得到一个特别经典函数

（2）数值计算复验

tt=0:0.001:4.5;tt(1)=eps;

yn=trapz(sin(tt)./tt)*0.001yn=

1.6541

??/12在n?0的限制下，求y(n)?13?20sinnxdx的一般积分表达式，

并计算y()的32位有效数字表达。

〖目的〗

?一般符号解与高精度符号数值解。

11

〖解答〗

symsx

symsnpositivef=sin(x)^n;

yn=int(f,x,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.有序列x(k)?ak，h(k)?bk，（在此k?0，a?b），求这两个序

列的卷积y(k)??h(n)x(k?n)。

n?0k〖目的〗

?符号离散卷积直接法和变换法。〖解答〗（1）直接法

symsabknx=a^k;h=b^k;

w=symsum(subs(h,k,n)*subs(x,k,k-n),n,0,k)%据定义y1=simple(w)w=

piecewise([a=b,b^k+b^k*k],[ab,(a*a^k-b*b^k)/(a-b)])y1=

piecewise([a=b,b^k+b^k*k],[ab,(a*a^k-b*b^k)/(a-b)])

（2）变换法（复验）

symsz

X=ztrans(a^k,k,z);H=ztrans(b^k,k,z);y2=iztrans(H*X,z,k)y2=

%通过Z变换及反变换

piecewise([b0,(a*a^k)/(a-b)-(b*b^k)/(a-b)])

〖说明〗

12

?符号计算不同途径产生的结果在形式上有可能不同，而且往往无法依靠符号计算本身的

指令是它们一致。此时，必须通过手工解决。

t?0?14.设系统的冲激响应为h(t)?e?3t，求该系统在输入u(t)?cost，

作用下的输出。

〖目的〗

?符号连续函数卷积的直接法和变换法。?符号变量限定性定义的作用。?laplace,ilaplace指令的应用。〖解答〗（1）直接法

symst

h=exp(-3*t);u=cos(t);symstao;

h_tao=subs(h,t,tao);u_t_tao=subs(u,t,t-tao);hu_tao=h_tao*u_t_tao;

hut=simple(int(hu_tao,tao,0,t))%直接卷积hut=

(3*cos(t))/10-3/(10*exp(3*t))+sin(t)/10

（2）变换法（复验）

symss;

HU=laplace(h,t,s)*laplace(u,t,s);huL=simple(ilaplace(HU,s,t))huL=

%拉氏变换及反变换

(3*cos(t))/10-3/(10*exp(3*t))+sin(t)/10

?15.求f(t)?Ae??t,??0的Fourier变换。

〖目的〗

?符号变量限定性定义的作用。?fourier指令的应用。〖解答〗

symsAtw

a=sym('a','positive');f=A*exp(-a*abs(t));y=fourier(f,t,w)F=simple(y)y=

(2*A*a)/(a^2+w^2)F=

13

(2*A*a)/(a^2+w^2)

??t??A1??16.求f(t)??????0?????t??的Fourier变换，并画出A?2,??2时t??的幅频谱。

〖目的〗

?单位阶跃符号函数heaviside的应用。?subs实现多变量置换。?ezplot的使用。〖解答〗

symstAw;

tao=sym('tao','positive');

f=A*((1+t/tao)*(heaviside(t+tao)-heaviside(t))+(1-t/tao)*(heaviside(t)-heaviside(t-tao)));Fw=fourier(f,t,w);Fws=simple(Fw)

Fw2=subs(Fws,[A,tao],[2,2])ezplot(abs(Fw2))gridFws=

-(4*A*(cos((tao*w)/2)^2-1))/(tao*w^2)Fw2=

-(8*cos(w)^2-8)/(2*w^2)

14

abs(8cos(w)2-8)/(2abs(w)2)43.532.521.510.50-6-4-20w246

?17.求F(s)?〖解答〗

symsst

s?3的Laplace反变换。

s3?3s2?6s?4F=(s+3)/(s^3+3*s^2+6*s+4);f=simple(ilaplace(F,s,t))f=

(3^(1/2)*sin(3^(1/2)*t)-2*cos(3^(1/2)*t)+2)/(3*exp(t))

?18.利用符号运算证明Laplace变换的时域求导性质：

?df(t)?L??s?L?f(t)??f(0)。??dt?〖目的〗

?符号计算用于定理证明。〖解答〗

symsts;y=sym('f(t)');df=diff(y,t);

Ldy=laplace(df,t,s)

15

Ldy=

s*laplace(f(t),t,s)-f(0)

??kT?19.求f(k)?ke的Z变换表达式。

〖目的〗

?注意：变换中，被变换变量的约定。〖解答〗

symslambdakTz;f_k=k*exp(-lambda*k*T);F_z=simple(ztrans(f_k,k,z))F_z=

(z*exp(T*lambda))/(z*exp(T*lambda)-1)^2

?20.求方程x2?y2?1,xy?2的解。

〖目的〗

?solve指令中，被解方程的正确书写，输出量的正确次序。〖解答〗

eq1='x^2+y^2=1';eq2='x*y=2';

[x,y]=solve(eq1,eq2,'x','y')x=

(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)/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)(1/2-(15^(1/2)*i)/2)^(1/2)-(1/2-(15^(1/2)*i)/2)^(1/2)

?21.求图p2-1所示信号流图的系统传递函数，并对照胡寿松主编

“自动控制原理〞中的例2-21结果，进行局部性验证。

16

图p2-1

〖目的〗

?理解和把握信号流图传递函数的“代数状态方程解法〞。

?并设法用胡寿松主编的“自动控制原理〞的例2-21进行局部性验证。

〖解答〗

（1）求传递函数

symsG1G2G3G4G5G6G7H1H2H3H4H5A=[

0000-H3-H4;0G200-H2G6;

G10-H1000;00G300G7;000G400;0G5000-H5];b=[1;0;0;0;0;0];c=[000010];Y2U=c*((eye(size(A))-A)\\b);%求传递函数[NN,DD]=numden(Y2U);DD=sort(DD);

%分开出分子、分母多项式

%分母多项式排序

disp([blanks(5),'传递函数Y2U为'])pretty(NN/DD)传递函数Y2U为

(G1G4(G2G3+G5G7+G3G5G6+G2G3H5))/

(H5+G2H1+G3G4H2+G1G5H4+G5G6H1+G2H1H5+G3G4H2H5+

G1G2G3G4H3+G1G4G5G7H3-G4G5G7H1H2+G1G3G4G5G6H3

17

+

G1G2G3G4H3H5+G1G3G4G5H2H4+1)

（2）局部性验证

symsabcdefg

y2u=subs(Y2U,[G1,G2,G3,G4,G5,G6,G7,H1,H2,H3,H4,H5],[a,e,f,1,b,c,0,g,0,0,0,d]);

[nn,dd]=numden(y2u);dd=sort(dd);

disp([blanks(5),'局部性验证用的传递函数y2u'])pretty(nn/dd)

局部性验证用的传递函数y2u

a(ef+bcf+def)d+eg+bcg+deg+1

此结果与胡寿松主编的“自动控制原理〞例2-21一致。

?22.采用代数状态方程法求图p2-2所示结构框图的传递函数

Y。WY和U

图p2-2

〖目的〗

?运用“代数状态方程解法〞求输入和扰动同时存在的结构框图的传递函数。

〖解答〗

（1）理论演绎对于结构框图写出状态方程

18

?x?Ax?bU?fW??Y?cx?dU?gW此式第一个方程关于x的解可写为

（p2-1）

x?(I?A)?1bU?(I?A)?1fW

把此式代入式（p2-1）的其次个方程，加以整理后可得

（p2-2）

Y?c(I?A)?1b?dU?c(I?A)?1f?gW

据此可写出传递函数

????Y?c(I?A)?1b?dU

（p2-3）（p2-4）

Y?c(I?A)?1f?gN（2）列出“元素级〞状态方程值得提醒：在编写M码之前，最好先在草稿纸上，细心“元素级〞状态方程是避免出错的冲要措施。对此，不要掉以轻心。本例的“元素级〞状态方程如下

?x1??0?x??G?2??2?x3???0????x4??0??0?x5???Y??01?G1??x1??G1??0??x??0??0?0?2???????0000???x3???0??U??G3??W（p2-5）???????H1000??x4??0??0??H2000??0??????H2???x5???000??x?0?U?(?1)?W000?G2?G10（3）编写相应的M码

symsG1G2G3H1H2A=[

000-G1-G1;00000;

G20-G200;0H1000;0H2000];b=[G1;0;0;0;0];f=[0;0;G3;0;-H2];c=[01000];d=0;g=-1;

R=c/(eye(size(A))-A);Y2U=R*b+d;Y2W=R*f+g;

%中间变量%计算传递函数Y/U%计算传递函数Y/W%分开出分子、分母多项式

%分母多项式排序

[NU,DU]=numden(Y2U);DU=sort(DU);

disp([blanks(5),'传递函数Y2U为'])

19

pretty(NU/DU)[NW,DW]=numden(Y2W);NW=sort(NW);DW=sort(DW);

disp([blanks(5),'传递函数Y2W为'])pretty(NW/DW)传递函数Y2U为

G1G2

G1G2H1+G1G2H2+1传递函数Y2W为

G2G3+G1G2H1+1-G1G2H1+G1G2H2+1

?23.求微分方程yy?5?x4?0的通解，并绘制任意常数为1时解的图

形。

〖目的〗

?理解指令dsolve的正确使用。?对dsolve输出结果的正确理解。

?ezplot指令绘图时，如何进行线色控制。?如何覆盖那些不能反映图形窗内容的图名。〖解答〗（1）求通解

reset(symengine)clearsymsyx

y=dsolve('0.2*y*Dy+0.25*x=0','x')y=

2^(1/2)*(C3-(5*x^2)/8)^(1/2)-2^(1/2)*(C3-(5*x^2)/8)^(1/2)

（2）根据所得通解中不定常数的符号写出“对其进行数值替代的指令〞

yy=subs(y,'C3',1)yy=

2^(1/2)*(1-(5*x^2)/8)^(1/2)-2^(1/2)*(1-(5*x^2)/8)^(1/2)

%将通解中的C3用1代替

（3）观测通解中两个分解的平方是否一致yy(1)^2==yy(2)^2

ans=

20

?1??2??11.求矩阵Ax?b的解，A为4阶魔方阵，b???。

?3????4?〖提醒〗

?由rref可以看出A不满秩，b不在A的值空间中，方程没有确凿解。?但可求最小二乘近似解。

〖目的〗

?A不满秩，b不在A的值空间中，方程没有确凿解。〖解答〗（1）借助增广矩阵用指令rref求解

A=magic(4);b=(1:4)';

%产生3阶魔方阵

[R,C]=rref([A,b])R=

%求解，并判断解的唯一性

1001001030001-3000001C=

1235

（2）用伪逆求最小二乘近似解（超出范围，仅供参考。）

x=pinv(A)*bb_pinv=A*xx=0.02350.12350.12350.0235b_pinv=1.30002.90002.10003.7000

%非确凿解%验算

〖解答〗

?C说明，A的秩为3，A不满秩；R第5列第4元素非零，说明b不在A的值空间中，

因此方程没有确凿解，但可以求最小二乘近似解。

?12.求?0.5?t?10e?0.2tsin[sint]?0的实数解。

46

〖提醒〗

?在适当范围内，作图观测一元繁杂函数的形态：观测解的存在性；解的唯一性。进而，

借助图形法求近似解。?匿名函数的使用方法。?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))gridon,shg

%利用匿名函数求y函数值%从图形获得近似解

[tt1,yy1]=ginput(1)tt1=2.7370yy1=0.0097

420-2-4-6-8-10-12-1012345

（2）进一步利用fzero求确切解

[t,yt]=fzero(y,tt1)t=2.7341yt=

2.2204e-015

47

〖说明〗

?假使在从图上获取数据前，先把零点附近图形放大，可以得到精度更高的近似解。

?13.求解二元函数方程组?〖目的〗

?solve指令的用法。〖解答〗（1）符号法只能得到两组解

?sin(x?y)?0的解。

?cos(x?y)?0S=solve('sin(x-y)','cos(x+y)','x','y');X=S.x,Y=S.yX=

[1/4*pi][-1/4*pi]Y=

[1/4*pi][-1/4*pi]

（2）数值法

Z1=sin(X-Y);Z2=cos(X+Y);

可以看到大量解，并从中找到规律

x=-6:0.1:6;y=x;[X,Y]=meshgrid(x,y);

contour(X,Y,Z1,3)%在Z1的取值范围内取3个等分点作为等位线取值。axissquareholdon

contour(X,Y,Z2,3)%注意：所有绿线交点都是解holdoffgridon;shg

%由于Z1关于z1=0对称，所以中间线，即Z1=0的等位线。

48

?14假定某窑工艺瓷器的烧制成品合格率为0.157，现该窑烧制100

件瓷器，请画出合格产品数的概率分布曲线。

〖目的〗

?指令binopdf的用法。?绘图指令stem的用法。〖解答〗

N=100;p=0.157;k=0:N;

%给定二项分布的特征参数%定义事件A发生的次数数组%算出各发生次数的概率

pdf=binopdf(k,N,p);stem(k,pdf)gridonshg

49

0.120.10.080.060.040.020230406080100

?15试产生均值为4，标准差为2的(10000?1)的正态分布随机数组

a，分别用hist和histfit绘制该数组的频数直方图，观测两张图形的差异。除histfit上的拟合红线外，你能使这两个指令汇出一致的频数直方图吗？

〖目的〗

?指令normrnd的应用，及随机发生器的初始状态控制。?指令hist与histfit的异同。〖解答〗

（1）绘制频数直方图

rngdefault

%为本例结果可被读者重现而设，并非必要。

mu=4;sigma=2;

%定义均值和标准差

%a=4+2*randn(10000,1);

a=normrnd(mu,sigma,10000,1);

%产生正态分布随机数组a

subplot(2,1,1),hist(a)%绘制a的频数统计直方图subplot(2,1,2),histfit(a)

50

=1？');%该指令只能在指令窗中运行tt=linspace(-5,5,N+1);fork=1:N

[tmin(k),fobj(k)]=fminbnd(ft,tt(k),tt(k+1));end

[fobj,ii]=sort(fobj);%将目标值由小到大排列tmin=tmin(ii);%使微小值点做与目标值相应的重新排列fobj,tmin

（2）最终确定的最小值点在N?1,2,?,10的不同分割下，经观测，最终确定出最小值点是-1.28498111480531相应目标值是-0.18604801006545

（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)gridon,shg

（B）经观测后，把最小值附近邻域放到足够大，然后运行以下指令，那放大图形被推向前台，与此同时光标变为“十字线〞，利用它点击极值点可得到最小值数据

[tmin2,fobj2]=ginput(1)

234

tmin2=

-1.28500000993975fobj2=

-0.18604799369136

出现具有一致数值的刻度区域说明已达最小可分辩状态

（4）符号法求最小值的尝试

symst

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=

-.60100931947716486053884417850955e-2fobj3=

.89909908144684551670208797723124

〖说明〗

35

?最小值是对整个区间而言的，微小值是对邻域而言的。?在一个区间中寻觅最小值点，对不同子区间分割进行屡屡探寻是必要的。这样可以避免

把微小值点误作为最小值点。最小值点是从一系列微小值点和边界点的比较中确定的。?作图法求最小值点，很直观。假若绘图时，自变量步长取得足够小，那么所求得的最小

值点有相当好的精度。?符号法在本例中，只求出一个极值点。其余好多极值点无法秋初，更不可能得到最小值。

d2y(t)dy(t)dy(0)?3?2y(t)?1,y(0)?1,?0，用数值法和符号法求?6.设2dtdtdty(t)t?0.5。

〖目的〗

?学习如何把高阶微分方程写成一阶微分方程组。

?ode45解算器的导数函数如何采用匿名函数形式构成。

?如何从ode45一组数值解点，求指定自变量对应的函数值。〖解答〗

（1）改写高阶微分方程为一阶微分方程组

令y1(t)?y(t),y2(t)?dy(t)，于是据高阶微分方程可写出dtdy1(t)??y2(t)?dt?dy(t)?2??2y1(t)?3y2(t)?1?dt

（2）运行以下指令求y(t)的数值解

formatlongts=[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）符号法求解

symst;

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=

36

.78958035647060552916850705213780

〖说明〗

?第条指令中的导数函数也可采用M函数文件表达，具体如下。

functionS=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);[R,ci]=rref(A);B1=A(:,ci)

B2=orth(A)[V,D]=eig(A);B3=V(:,1:rv)B1=

64239555417474640262732343541232249151485859B2=

-0.35360.54010.3536-0.3536-0.3858-0.3536-0.3536-0.2315-0.3536-0.35360.07720.3536-0.3536-0.07720.3536-0.35360.2315-0.3536-0.35360.3858-0.3536-0.3536-0.54010.3536B3=

%采用8阶魔方阵作为试验矩阵

%直接从A中取基向量%求A值空间的正交基

%非零特征值数就是矩阵的秩

%取A的非零特征值对应的特征向量作基

rv=sum(sum(abs(D))>1000*eps);

37

0.35360.62700.39130.3536-0.4815-0.24580.3536-0.3361-0.10040.35360.1906-0.04510.35360.0451-0.19060.35360.10040.33610.35360.24580.48150.3536-0.3913-0.6270

（2）验证A的任何列可用B1线性表出

B1_A=rref([B1,A])%若B1_A矩阵的下5行全为0，B1_A=

%就说明A可以被B1的3根基向量线性表出

1001001100101001034-3-47001001-3-445-70000000000000000000000000000000000000000000000000000000

B2_A=rref([B2,A])B2_A=

Columns1through7

1.000000-91.9239-91.9239-91.9239-91.923901.0000051.8459-51.8459-51.845951.8459001.00009.8995-7.0711-4.24261.414200000000000000000000000000000000000Columns8through11

-91.9239-91.9239-91.9239-91.923951.8459-51.8459-51.845951.8459-1.41424.24267.0711-9.899500000000000000000000

38

B3_A=rref([B3,A])B3_A=

Columns1through7

1.00000091.923991.923991.923991.923901.0000042.3447-38.1021-33.859429.6168001.000012.6462-16.8889-21.131525.374100000000000000000000000000000000000Columns8through11

91.923991.923991.923991.923925.3741-21.1315-16.888912.646229.6168-33.8594-38.102142.344700000000000000000000

〖说明〗

?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=

%为紧凑地显示特征值而写

-911-2163-25270-69141-4211684-575575-11493451-138013891-38917782-2334593365

39

1024-10242048-614424572ans=

Columns1through4

-0.0181-0.0054-0.0171i-0.0054+0.0171i0.0144-0.0104iColumn5

0.0144+0.0104i

（2）验算说明相对误差较大

AE=V*D/V

er_AE=norm(A-AE,'fro')/norm(A,'fro')AE=

1.0e+004*

Columns1through4

-0.0009+0.0000i0.0011-0.0000i-0.0021+0.0000i0.0063-0.0000i

0.0070-0.0000i-0.0069+0.0000i0.0141-0.0000i-0.0421+0.0000i

-0.0575+0.0000i0.0575-0.0000i-0.1149+0.0000i0.3451-0.0000i

0.3891-0.0000i-0.3891+0.0000i0.7781-0.0000i-2.3343+0.0000i

0.1024-0.0000i-0.1024+0.0000i0.2048-0.0000i-0.6144+0.0000iColumn5

-0.0252+0.0000i0.1684-0.0000i-1.3800+0.0000i9.3359-0.0001i2.4570-0.0000ier_AE=

6.9310e-005

%相对F范数

（3）一个更严谨的特征值分解指令

[Vc,Dc,eigc]=condeig(A)Vc=

Columns1through4

-0.0000-0.0000+0.0000i-0.0000-0.0000i0.0000+0.0000i

0.02060.0207+0.0000i0.0207-0.0000i0.0207+0.0000i-0.1397-0.1397+0.0000i-0.1397-0.0000i-0.1397+0.0000i

0.95740.95740.95740.95740.25190.2519-0.0000i0.2519+0.0000i0.2519-0.0000i

%eigc中的高值时，说明相应的特征值不可信。

40

Column5

0.0000-0.0000i0.0207-0.0000i-0.1397-0.0000i0.95740.2519+0.0000iDc=

Columns1through4

-0.01810000-0.0054+0.0171i0000-0.0054-0.0171i00000.0144+0.0104i0000Column500000.0144-0.0104ieigc=1.0e+011*5.26875.23135.23135.17255.1724

（4）对A采用Jordan分解并检验

[VJ,DJ]=jordan(A);DJ

AJ=VJ*DJ/VJ

er_AJ=norm(A-AJ,'fro')/norm(A,'fro')DJ=

0100000100000100000100000AJ=

1.0e+004*

-0.00090.0011-0.00210.0063-0.02520.0070-0.00690.0141-0.04210.1684-0.05750.0575-0.11490.3451-1.38010.3891-0.38910.7782-2.33459.3365

%求出确凿的广义特征值，使A*VJ=VJ*D成立。

41

0.1024-0.10240.2048-0.61442.4572er_AJ=

2.0500e-011

〖说明〗

?指令condeig的第3输出量eigc给出的是所谓的“矩阵特征值条件数〞。当特征条件数

与1eps相当时，就意味着矩阵A可能“退化〞，即矩阵可能存在“代数重数〞大于

“几何重数〞的特征值。此时，实施Jordan分解更适合。?顺便指出：借助condeig算得的特征值条件数与cond指令算得的矩阵条件数是两个不同

概念。前者描述特征值的问题，后者描述矩阵逆的问题。

?本例矩阵A的特征值条件数很高，说明分解不可信。验算也说明，相对误差较大。?当对矩阵A进行Jordan分解时，可以看到，A具有5重根。当对Jordan分解进行验算

时，相对误差很小。

?9.求矩阵Ax?b的解，A为3阶魔方阵，b是(3?1)的全1列向量。

〖提醒〗

?了解magic指令?rref用于方程求解。

?矩阵除法和逆阵法解方程。

〖目的〗

?满秩方阵求解的一般过程。?rref用于方程求解。

?矩阵除法和逆阵法解方程。〖解答〗

A=magic(3);

%产生3阶魔方阵%（3*1）全1列向量%矩阵除解方程%逆阵法解方程

b=ones(3,1);x=A\\b

[R,C]=rref([A,b])%GaussJordan消去法解方程，同时判断解的唯一性xx=inv(A)*bR=

1.0000000.066701.000000.0667001.00000.0667C=

123x=0.06670.06670.0667xx=0.0667

42

0.06670.0667

〖说明〗

?rref指令通过对增广矩阵进行消去法操作完成解方程。由分解得到的3根“坐标向量〞

和（或）C3指示的3根基向量，可见A3满秩，因此方程解唯一。?在本例状况下，矩阵除、逆阵法、rref法所得解一致。

?10.求矩阵Ax?b的解，A为4阶魔方阵，b是(4?1)的全1列向量。

〖提醒〗

?用rref可观测A不满秩，但b在A的值空间中，这类方程用无数解。?矩阵除法能正确求得这类方程的特解。?逆阵法不能求得这类方程的特解。?注意特解和齐次解

〖目的〗

?A不满秩，但b在A的值空间中，这类方程的求解过程。?rref用于方程求解。

?矩阵除法能正确求得这类方程的特解。?逆阵法不能求得这类方程的特解。?解的验证方法。?齐次解的获取。?全解的获得。〖解答〗

（1）借助增广矩阵用指令rref求解

A=magic(4);

%产生3阶魔方阵%全1列向量

b=ones(4,1);

[R,C]=rref([A,b])R=

%求解，并判断解的唯一性

1.0000001.00000.058801.000003.00000.1176001.0000-3.0000-0.058800000C=

123

关于以上结果的说明：?R阶梯阵提供的信息

?前4列是原A阵经消元变换后的阶梯阵；而第5列是原b向量经一致变换后的结

果。

?R的前三列为“基〞，说明原A阵秩为3；而第4列的前三个元素，表示原A阵

的第4列由其前三列线性组合而成时的加权系数，即方程的一个解。

43

?R的第5列说明：b可由原A阵的前三列线性表出；b给出了方程的一个解；由于

原A阵“缺秩〞，所以方程的确解不唯一。

?C数组提供的信息

?该数组中的三个元素表示变换取原A阵的第1，2，3列为基。?该数组的元素总数就是“原A阵的秩〞

（2）矩阵除求得的解

x=A\\b

Warning:Matrixisclosetosingularorbadlyscaled.Resultsmaybeinaccurate.RCOND=1.306145e-017.x=0.05880.1176-0.05880

运行结果指示：矩阵除法给出的解与rref解一致。（实际上，MATLAB在设计“除法〞程序时，针对“b在A值空间中〞的状况，就是用rref求解的。）

（3）逆阵法的解

xx=inv(A)*b

Warning:Matrixisclosetosingularorbadlyscaled.Resultsmaybeinaccurate.RCOND=1.306145e-017.xx=0.0469