MATLAB--R2012a课后普习题答案全解

●MATLAB R2012a
●课后习题答案全解
第一章基础准备及入门
习题1及解答
⏹ 1.数字1.5e2，1.5e3 中的哪个与1500相同吗？
〖解答〗
1.5e3
⏹ 2.请指出如下5个变量名中，哪些是合法的？
abcd-2 xyz_3 3chan a变量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) %仅求出主根
w2=(a^2)^(1/3) %求出(-8)^2的主根 w3=(a^(1/3))^2
%求出(-8)主根后再平方
w1 = -2.0000 + 3.4641i
w2 = 4.0000 w3 =
-2.0000 + 3.4641i
（2）复数的多方根的，下面是求取全部方根的两种方法： （A ）根据复数方根定义
a=-8;n=2;m=3;
ma=abs(a);aa=angle(a); for k=1:m
%m 决定循环次数 sa(k)=(aa+2*pi*(k-1))*n/m;
%计算各根的相角 end
result=(ma^(2/3)).*exp(j*sa) %计算各根
result =
-2.0000 + 3.4641i 4.0000 - 0.0000i -2.0000 - 3.4641i
（B ）利用多项式02
3=-a r 求根
p=[1,0,0,-a^2]; r=roots(p) r =
-2.0000 + 3.4641i -2.0000 - 3.4641i 4.0000
⏹ 5.指令clear, clf, clc 各有什么用处？
〖解答〗 clear 清除工作空间中所有的变量。

clf 清除当前图形。

clc 清除命令窗口中所有显示。

⏹ 6.以下两种说法对吗？（1）“MATLAB 进行数值的表达精度与
其指令窗中的数据显示精度相同。

” （2）MATLAB 指令窗中显示的数值有效位数不超过7位。

”
〖解答〗
（1）否；（2）否。

⏹ 7.想要在MATLAB 中产生二维数组⎥⎥⎥⎦
⎤⎢⎢⎢⎣⎡=987654321S ，下面哪些指令能实现目的？
（1） S=[1,2,3;4,5,6;7,8;9]
（2） S=[1 2 3;4 5 6;7 8 9]
（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 指令。

〖解答〗
c1=3/7+0.1
c2=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.5286
c2 =
37/70
c3 =
0.52857142857142857142857142857143
c4 =
0.52857142857142857142857142857143
Cs1 =
double
Cs2 =
sym
Cs3 =
sym
Cs4 =
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
symvar(sym('z*exp(j*th)'),1) ans = z
⏹ 3求以下两个方程的解
（1）试写出求三阶方程05.443
=-x 正实根的程序。

注意：只要正
实根，不要出现其他根。

（2）试求二阶方程022=+-a ax x 在0>a 时的根。

〖目的〗
● 体验变量限定假设的影响 〖解答〗
（1）求三阶方程05.443
=-x 正实根
reset(symengine)
%确保下面操作不受前面指令运作的影响
syms x positive solve(x^3-44.5) ans =
(2^(2/3)*89^(1/3))/2
（2）求五阶方程02
2=+-a ax x 的实根
syms a positive
%注意：关于x 的假设没有去除
solve(x^2-a*x+a^2) Warning: Explicit solution could not be found.
> In solve at 83 ans =
[ empty sym ]
syms x clear syms a positive solve(x^2-a*x+a^2) 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))等来观察。

〖目的〗
●理解准确符号数值的创建法。

●高精度误差的观察。

〖解答〗
（1）x=7/3
x=7/3;a=x,b=sym(x),c=sym(x,'d'),d=sym('7/3'),
a =
2.3333
b =
7/3
c =
2.3333333333333334813630699500209
d =
7/3
v1=vpa(abs(a-d)),v2=vpa(abs(b-d)),v3=vpa(abs(c-d))
v1 =
0.0
v2 =
0.0
v3 =
0.00000000000000014802973661668756666666667788716
（2）x=pi/3
x=pi/3;a=x,b=sym(x),c=sym(x,'d'),d=sym('pi/3'),
a =
1.0472
b =
pi/3
c =
1.047197551196597631317786181171
d =
pi/3
v1=vpa(abs(a-d)),v2=vpa(abs(b-d)),v3=vpa(abs(c-d)) v1 = 0.0 v2 = 0.0 v3 =
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.5310
b =
1275352044764433/281474976710656 c =
4.5309606547207899041040946030989 d =
pi*3^(1/3)
v1=vpa(abs(a-d)),v2=vpa(abs(b-d)),v3=vpa(abs(c-d)) v1 =
0.00000000000000026601114166290944374842393221638 v2 =
0.00000000000000026601114166290944374842393221638 v3 =
0.0000000000000002660111416629094726767991785515
⏹ 5求符号矩阵⎥⎥⎥⎦
⎤⎢⎢⎢⎣⎡=3332
31
2322
21
131211
a a a a a a a a a A 的行列式值和逆，所得结果应采用“子表达式置换”简洁化。

〖目的〗
● 理解subexpr 指令。

〖解答〗
A=sym('[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, a33] DA =
a11*a22*a33 - a11*a23*a32 - a12*a21*a33 + a12*a23*a31 + a13*a21*a32 - a13*a22*a31 IAs =
[ d*(a22*a33 - a23*a32), -d*(a12*a33 - a13*a32), d*(a12*a23 - a13*a22)] [ -d*(a21*a33 - a23*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求∑∞
=0k k x 的符号解，并进而用该符号解求∑∞
=-0)31(k k ，∑∞
=0)1(k k
π
，
∑∞
=0
3
k k
的准确值。

〖目的〗
● symsum, subs 的应用。

●
从实例中，感受指令所给出的关于∑∞
=0k k x 符号解的含义。

〖解答〗
syms x k f=x^(k);
Z1=symsum(f,k,0,inf) Z1 =
piecewise([1 <= x, Inf], [abs(x) < 1, -1/(x - 1)]) %piecewise 分段函数
subs(Z1,x,{sym('-1/3'),sym('1/pi'),sym('3')}) ans =
[ 3/4, -1/(1/pi - 1), Inf]
⏹ 7对于0>x ，求1
20
11122+∞
=∑⎪
⎭⎫
⎝⎛+-+k k x x k 。

（提示：理论结果为x ln ）
〖目的〗
● 符号变量的限定性定义的作用。

〖解答〗
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,inf)) %结果与理论值lnx 相符！ s =
piecewise([abs(x - 1) < x + 1, log(x)]) %abs 绝对值
〖注意〗
● 解答中，条件abs(x - 1) < x + 1意味着：
⏹ 约束一：x-1<x+1 ⇒ 2>0 ⇒
此式总成立，说明“无约束”。

⏹ 情况二：-(x-1)<x+1 ⇒ x>0 ⇒ 此为“约束”，满足题意。

⏹ 8（1）通过符号计算求t t y sin )(=的导数
dt
dy。

（2）然后根据此结果，求-
=0t dt dy
和2
π=t dt dy 。

〖目的〗
● diff, limit 指令的应用。

● 如何理解运行结果。

〖解答〗
syms t
y=abs(sin(t))
d=diff(y) %求dy/dt
d0_=limit(d,t,0,'left') %求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_ = -1 dpi_2 = 0
⏹ 9求出dx x e x sin 7.110⎰--π
π的具有64位有效数字的积分值。

〖目的〗
● 符号积分的解析解和符号数值解。

● 符号计算和数值计算的相互校验。

〖解答〗
（1）符号积分
syms x clear syms x
y=exp(-abs(x))*abs(sin(x))
si=vpa(int(y,-10*pi,1.7*pi),64) %vpa 指定精确位数 y =
abs(sin(x))/exp(abs(x)) si =
1.087849499412904913166671875948174520895458535212845987519414166
（2）数值计算复验
xx=-10*pi:pi/100:1.7*pi;
sn=trapz(exp(-abs(xx)).*abs(sin(xx)))*pi/100 %trapz 梯形法求积分 sn =
1.0877
⏹ 10计算二重积分⎰⎰
+2
1
1
222
)(x dydx y x 。

〖目的〗
● 变上限二重积分的符号计算法。

〖解答〗
syms x y f=x^2+y^2;
r=int(int(f,y,1,x^2),x,1,2) r =
1006/105
⏹ 11在]2,0[π区间，画出dt t
t
x y x
⎰
=
sin )(曲线，并计算)5.4(y 。

〖目的〗
● 在符号计算中，经常遇到计算结果是特殊经典函数的情况。

● 如何应用subs 获得超过16位有效数字的符号数值结果。

● 初步尝试ezplot 指令的简便。

〖解答〗
（1）符号计算
syms t x;
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
1
2
34
5
6
00.20.40.60.811.21.41.61.8x
sinint(x)
（2）数值计算复验
tt=0:0.001:4.5; tt(1)=eps;
yn=trapz(sin(tt)./tt)*0.001 yn =
1.6541
⏹ 12在0>n 的限制下，求xdx n y n ⎰
=
20
sin )(π
的一般积分表达式，
并计算)3
1(y 的32位有效数字表达。

〖目的〗
● 一般符号解与高精度符号数值解。

〖解答〗
syms x
syms n positive
f=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)/2
y3s =
1.2935547796148952674767575125656
y3d =
1.2935547796148951782413405453553
⏹13.有序列k a
k
x=
)
(，k b
k
h=
)
(，（在此0
≥
k，b
a≠），求这两个序
列的卷积∑
=-
=
k
n
n
k
x
n
h
k
y
)
(
)
(
)
(。

〖目的〗
●符号离散卷积直接法和变换法。

〖解答〗
（1）直接法
syms a b k n
x=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], [a <> b, (a*a^k - b*b^k)/(a - b)]) y1 =
piecewise([a = b, b^k + b^k*k], [a <> b, (a*a^k - b*b^k)/(a - b)])
（2）变换法（复验）
syms z
X=ztrans(a^k,k,z);
H=ztrans(b^k,k,z);
y2=iztrans(H*X,z,k) %通过Z变换及反变换
y2 =
piecewise([b <> 0, (a*a^k)/(a - b) - (b*b^k)/(a - b)])
〖说明〗
● 符号计算不同途径产生的结果在形式上有可能不同，而且往往无法依靠符号计算本身的
指令是它们一致。

此时，必须通过手工解决。

⏹ 14.设系统的冲激响应为t e t h 3)(-=，求该系统在输入t t u cos )(=，
0≥t 作用下的输出。

〖目的〗
● 符号连续函数卷积的直接法和变换法。

● 符号变量限定性定义的作用。

● laplace, ilaplace 指令的应用。

〖解答〗 （1）直接法
syms t
h=exp(-3*t);u=cos(t); syms tao;
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）变换法（复验）
syms s;
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.求0,)(>=-ααt Ae t f 的Fourier 变换。

〖目的〗
● 符号变量限定性定义的作用。

● fourier 指令的应用。

〖解答〗
syms A t w
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 =
(2*A*a)/(a^2 + w^2)
⏹ 16.求⎪⎩
⎪⎨⎧>≤⎪⎪⎭
⎫
⎝⎛-
=τ
τ
τt t t A t f 01)(的Fourier 变换，并画出2,2==τA 时的幅频谱。

〖目的〗
● 单位阶跃符号函数heaviside 的应用。

● subs 实现多变量置换。

● ezplot 的使用。

〖解答〗
syms t A w;
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)) grid Fws =
-(4*A*(cos((tao*w)/2)^2 - 1))/(tao*w^2) Fw2 =
-(8*cos(w)^2 - 8)/(2*w^2)
-6
-4
-2
02
4
6
00.511.522.533.54w
abs(8 cos(w)2 - 8)/(2 abs(w)2)
⏹ 17.求4
633
)(23++++=
s s s s s F 的Laplace 反变换。

〖解答〗
syms s t
F=(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 变换的时域求导性质：
[])0()()(f t f L s dt t df L -⋅=⎥⎦
⎤⎢⎣⎡。

〖目的〗
● 符号计算用于定理证明。

〖解答〗
syms t s; y=sym('f(t)'); df=diff(y,t);
Ldy=laplace(df,t,s)
Ldy =
s*laplace(f(t), t, s) - f(0)
⏹ 19.求T
k ke k f )(λ-=的Z 变换表达式。

〖目的〗
● 注意：变换中，被变换变量的约定。

〖解答〗
syms lambda k T z; 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.求方程2,122==+xy y x 的解。

〖目的〗
● 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)/2 y =
(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结果，进行局部性验证。

图p2-1
〖目的〗
●理解和掌握信号流图传递函数的“代数状态方程解法”。

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

〖解答〗
（1）求传递函数
syms G1 G2 G3 G4 G5 G6 G7 H1 H2 H3 H4 H5
A=[ 0 0 0 0 -H3 -H4;
G1 0 -H1 0 0 0;
0 G2 0 0 -H2 G6;
0 0 G3 0 0 G7;
0 0 0 G4 0 0;
0 G5 0 0 0 -H5];
b=[ 1; 0; 0; 0; 0; 0];
c=[ 0 0 0 0 1 0];
Y2U=c*((eye(size(A))-A)\b); %求传递函数
[NN,DD]=numden(Y2U); %分离出分子、分母多项式
DD=sort(DD); %分母多项式排序
disp([blanks(5),'传递函数 Y2U 为'])
pretty(NN/DD)
传递函数 Y2U 为
(G1 G4 (G2 G3 + G5 G7 + G3 G5 G6 + G2 G3 H5)) /
(H5 + G2 H1 + G3 G4 H2 + G1 G5 H4 + G5 G6 H1 + G2 H1 H5 + G3 G4 H2 H5 +
G1 G2 G3 G4 H3 + G1 G4 G5 G7 H3 - G4 G5 G7 H1 H2 + G1 G3 G4 G5 G6 H3
+
G1 G2 G3 G4 H3 H5 + G1 G3 G4 G5 H2 H4 + 1)
（2）局部性验证
syms a b c d e f g
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 (e f +
b
c f +
d
e f)
---------------------------
d +
e g + b c g + d e g + 1
此结果与胡寿松主编的“自动控制原理”例2-21一致。

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

W
图p2-2
〖目的〗
●运用“代数状态方程解法”求输入和扰动同时存在的结构框图的传递函数。

〖解答〗
（1）理论演绎
对于结构框图写出状态方程
⎩⎨
⎧++=++=W
U Y W
U g d cx f b Ax x （p2-1）
此式第一个方程关于x 的解可写为
W U f A)(I b A)(I x 11---+-=
（p2-2）
把此式代入式（p2-1）的第二个方程，加以整理后可得
[][]
W g U d Y +-++-=--f A)c(I b A)c(I 11
据此可写出传递函数
d U
Y
+-=-b A)c(I 1 （p2-3） g N
Y
+-=-f A)c(I 1
（p2-4）
（2）列出“元素级”状态方程 值得提醒：在编写M 码之前，最好先在草稿纸上，仔细“元素级”状态方程是避免出错的冲要措施。

对此，不要掉以轻心。

本例的“元素级”状态方程如下
[]W U Y W
H G U G x x x x x H H G G G G x x x x x ⋅-+⋅+⋅=⋅⎥⎥
⎥
⎥⎥
⎥⎦⎤
⎢⎢⎢⎢⎢⎢⎣⎡-+⋅⎥⎥⎥⎥⎥⎥⎦⎤⎢
⎢⎢⎢⎢⎢⎣⎡+⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⋅⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡---=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡)1(000010 0000000000000000000000
00023154321212
21154321x （p2-5） （3）编写相应的M 码
syms G1 G2 G3 H1 H2 A=[ 0 0 0 -G1 -G1; G2 0 -G2 0 0;
0 0 0 0 0;
0 H1 0 0 0; 0 H2 0 0 0]; b=[ G1; 0; 0; 0; 0]; f=[ 0; 0; G3; 0; -H2]; c=[ 0 1 0 0 0]; d=0; g=-1;
R=c/(eye(size(A))-A);
%中间变量 Y2U=R*b+d; %计算传递函数 Y/U Y2W=R*f+g;
%计算传递函数 Y/W [NU,DU]=numden(Y2U); %分离出分子、分母多项式 DU=sort(DU);
%分母多项式排序
disp([blanks(5),'传递函数 Y2U 为'])
pretty(NU/DU) [NW,DW]=numden(Y2W); NW=sort(NW); DW=sort(DW);
disp([blanks(5),'传递函数 Y2W 为']) pretty(NW/DW) 传递函数 Y2U 为
G1 G2
----------------------- G1 G2 H1 + G1 G2 H2 + 1 传递函数 Y2W 为
G2 G3 + G1 G2 H1 + 1 - ----------------------- G1 G2 H1 + G1 G2 H2 + 1
⏹ 23.求微分方程045=+'x y y 的通解，并绘制任意常数为1时解的图
形。

〖目的〗
● 理解指令dsolve 的正确使用。

● 对dsolve 输出结果的正确理解。

● ezplot 指令绘图时，如何进行线色控制。

● 如何覆盖那些不能反映图形窗内容的图名。

〖解答〗 （1）求通解
reset(symengine) clear syms y x
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) %将通解中的C3用1代替
yy =
2^(1/2)*(1 - (5*x^2)/8)^(1/2) -2^(1/2)*(1 - (5*x^2)/8)^(1/2)
（3）观察通解中两个分解的平方是否相同 yy(1)^2==yy(2)^2
ans =
1
（4）于是可考虑函数的平方关系
syms Y
fxy=Y^2-yy(1)^2 fxy =
Y^2 + (5*x^2)/4 - 2
（5）根据平方关系式画完整曲线
clf
ezplot(fxy,[-2,2,-2,2]) axis square grid on
Y
x
Y 2 + (5 x 2)/4 - 2 = 0
-2
-1.5-1-0.5
00.51 1.52
-2-1.5-1-0.500.511.52
（6）假如直接用“分解”画曲线，那么将是不完整的
ezplot(yy(1)),hold on cc=get(gca,'Children'); set(cc,'Color','r')
ezplot(yy(2)),axis([-2 2 -2 2]) legend('y(1)','y(2)'),hold off; title(' ')
%覆盖不完全的图名
grid
axis square
-2
-1.5-1-0.5
00.51 1.52
-2-1.5-1-0.500.511.52x
y(1)y(2)
⏹ 24.求一阶微分方程2)0(,2=+=x bt at x
的解。

〖目的〗
● 初值微分方程的符号解。

● pretty 指令的使用。

〖解答〗
x=dsolve('Dx=a*t^2+b*t','x(0)=2','t') pretty(x)
%比较易读的表达形式
x = (t^2*(3*b + 2*a*t))/6 + 2
2
t (3 b + 2 a t) ---------------- + 2 6
⏹ 25.求边值问题
1)0(,0)0(,34,43==+-=+=g f g f dx
dg g f dx df 的解。

（注意：相应的数值解法比较复杂）。

〖目的〗
●边值微分方程的符号解。

〖解答〗
[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)
第3章数值数组及其运算
习题3及解答
⏹ 5.要求在闭区间]
2,0[ 上产生具有10个等距采样点的一维数组。

试用两种不同的指令实现。

〖目的〗
●数值计算中产生自变量采样点的两个常用指令的异同。

〖解答〗
%方法一
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.2832
t2 =
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.2832
⏹ 6.由指令rng('default'),A=rand(3,5)生成二维数组A，试求该数组
中所有大于0.5的元素的位置，分别求出它们的“全下标”和“单下标”。

〖目的〗
●数组下标的不同描述：全下标和单下标。

●sub2ind, int2str, disp的使用。

●随机发生器的状态控制：保证随机数的可复现性。

〖解答〗
rng('default')
A=rand(3,5)
[ri,cj]=find(A>0.5);
id=sub2ind(size(A),ri,cj);
ri=ri';cj=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 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
⏹7.采用默认全局随机流，写出产生长度为1000的“等概率双位（即
取-1，+1）取值的随机码”程序指令，并给出-1码的数目。

〖目的〗
●两种基本随机发生器的使用。

●关系运算产生逻辑数组——可用于数组的元素的标识和寻访。

●逻辑数组的应用。

● 如何判断两个整数数组是否相等。

〖解答〗
（1）运用均匀随机数解题法——解法1
rng default %为以下结果重现而设；产生默认随机流。

详见第4.3.2节
A=rand(1,1000); a=2*(A>0.5)-1; Na=sum(a==-1) Na = 512
（2）运用正态随机数解题法——解法2
randn('state',123) B=randn(1,1000); b=2*(B>0)-1; Nb=sum(b==-1) Nb = 462
（3）直接发生法——解法3
c=randsrc(1,1000,[-1,1]); Nc=sum(c==-1) Nc = 482
⏹ 2.已知矩阵⎥⎦
⎤
⎢
⎣⎡=4321A ，运行指令B1=A.^(0.5), B2=A^(0.5), 可以观察到不同运算方法所得结果不同。

（1）请分别写出根据B1, B2恢复原矩阵A 的程序。

（2）用指令检验所得的两个恢复矩阵是否相等。

〖目的〗
● 数组运算和矩阵运算的不同。

● 如何判断两个双精度数组是否相等。

● norm 指令的应用。

〖解答〗
A=[1,2;3,4]; B1=A.^0.5 B2=A^0.5 A1=B1.*B1; A2=B2*B2;
norm(A1-A2,'fro') % 求误差矩阵的F-范数，当接近eps 量级时，就认为实际相等 B1 =
1.0000 1.4142 1.7321
2.0000 B2 =
0.5537 + 0.4644i 0.8070 - 0.2124i 1.2104 - 0.3186i 1.7641 + 0.1458i ans =
8.4961e-016
⏹ 4.在时间区间 [0,10]中，绘制
t e y t
2cos 15.0--=曲线。

要求分别采取“标量循环运算法”和“数组运算法”编写两段程序绘图。

〖目的〗
● 加强理解数组运算的机理和应用。

● 初步使用subplot, plot, xlabel, ylabel 等指令绘图。

〖解答〗
%标量循环运算法
t=linspace(0,10,200); N=length(t); y1=zeros(size(t)); for k=1:N
y1(k)=1-exp(-0.5*t(k))*cos(2*t(k)); end
subplot(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
510
00.5
11.5
t
y 1
510
00.5
1
1.5
t
y 2
⏹ 8.先运行clear,format long,rng('default'),A=rand(3,3)，然后根据
A写出两个矩阵：一个对角阵B，其相应元素由A的对角元素构成；另一个矩阵C，其对角元素全为0，而其余元素与对应的A 阵元素相同。

〖目的〗
●常用指令diag的使用场合。

〖解答〗
clear,
format long
rng('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 0.097540404999410 0.957506835434298
B =
0.814723686393179 0 0
0 0.632359246225410 0
0 0 0.957506835434298
C =
0 0.913375856139019 0.278498218867048
0.905791937075619 0 0.546881519204984
0.126986816293506 0.097540404999410 0
⏹9.先运行指令x=-3*pi:pi/15:3*pi; y=x; [X,Y]=meshgrid(x,y);
warning off; Z=sin(X).*sin(Y)./X./Y; 产生矩阵Z。

（1）请问矩阵Z中有多少个“非数”数据？（2）用指令surf(X,Y,Z); shading interp观察所绘的图形。

（3）请写出绘制相应的“无裂缝”图形的全部指令。

〖目的〗
●初步感受三维曲面的绘制方法。

●非数NaN的产生，非数的检测，和对图形的影响。

●sum的应用。

●eps如何克服“被零除”的尴尬。

〖解答〗
x=-3*pi:pi/15:3*pi; y=x;
[X,Y]=meshgrid(x,y); warning off
Z=sin(X).*sin(Y)./X./Y; NumOfNaN=sum(sum(isnan(Z))) %计算“非数”数目
subplot(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
10.下面有一段程序，企图用来解决如下计算任务：有矩阵
⎥⎥⎥⎥
⎦
⎤⎢⎢⎢
⎢⎣⎡++++=k k k k k k k k 10229221911 A ，当k 依次取10, 9, 8, 7, 6, 5, 4, 3, 2, 1时，计算矩阵k A “各列元素的和”，并把此求和结果存放为矩阵Sa
的第k 行。

例如3=k 时，A 阵为⎥⎥⎥⎦
⎤
⎢⎢⎢⎣⎡306329522841 ，
此时它各列元素 的和是一个)101(⨯行数组[]87156 ，并把它保存为Sa 的第3行。

问题：该段程序的计算结果对吗？假如计算结果不正确，请指出
错误发生的根源，并改正之。

〖目的〗
● 正确理解sum 的工作机理。

● reshape 的应用。

〖解答〗
（1）企图用以下程序完成题目要求。

for k=10:-1:1 A=reshape(1:10*k,k,10);
Sa(k,:)=sum(A);
end Sa 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 122 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时，]10,,2,1[1 A ，Sa 的第
1行应该与1A 相同。

上述程序的错误是对sum 理解不正确。

sum 对二维数组，求和按列施行；而对一维数组，不管行数组或列数组，总是求那数组所有元素的和。

正确的程序应该写成
for k=10:-1:1 A=reshape(1:10*k,k,10); Sa(k,:)=sum(A); if k==1
Sa(k,:)=A;
end
end Sa 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 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
第4章数值运算
习题 4 及解答
⏹ 1.根据题给的模拟实际测量数据的一组t和)(t y试用数值差分diff
或数值梯度gradient指令计算)(t
y'曲线绘制在
y'，然后把)(t y和)(t 同一张图上，观察数值求导的后果。

（模拟数据从prob_data401.mat获得）
〖目的〗
●强调：要非常慎用数值导数计算。

●练习mat数据文件中数据的获取。

●实验数据求导的后果
●把两条曲线绘制在同一图上的一种方法。

〖解答〗
（1）从数据文件获得数据的指令
假如prob_data401.mat文件在当前目录或搜索路径上
clear
load prob_data401.mat
（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 01234567
-2-1.5
-1
-0.5
0.5
11.5
（3）用gradent 求导的指令（图形与上相似）
dt=t(2)-t(1);
yc=gradient(y)/dt;
plot(t,y,'b',t,yc,'r')
grid on
〖说明〗
● 不到万不得已，不要进行数值求导。

● 假若一定要计算数值导数，自变量增量dt 要取得比原有数据相对误差高1、2个量级以
上。

● 求导会使数据中原有的噪声放大。

⏹ 2.采用数值计算方法，画出dt t
t x y x ⎰=
0sin )(在]10 ,0[区间曲线，并计算)5.4(y 。

〖提示〗
● 指定区间内的积分函数可用cumtrapz 指令给出。

● )5.4(y 在计算要求不太高的地方可用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 0246810
00.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
⏹ 3.求函数x e
x f 3sin )(=的数值积分⎰=π
0 )(dx x f s ，并请采用符号计算尝试复算。

〖提示〗
● 数值积分均可尝试。

● 符号积分的局限性。

〖目的〗
● 符号积分的局限性。

〖解答〗
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 sym.int at 58
ss =
int(exp(sin(x)^3),x = 0 .. pi)
⏹ 4.用quad 求取dx x e x sin 7.15⎰--π
π的数值积分，并保证积分的绝对精
度为910-。

〖目的〗
● quadl ，精度可控，计算较快。

● 近似积分指令trapz 获得高精度积分的内存和时间代价较高。

〖解答〗
%精度可控的数值积分
fx=@(x)exp(-abs(x)).*abs(sin(x));
format long
sq=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
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
⏹ 5.求函数5.08.12cos 5.1)5(sin )(2
06.02++-=t t t e t t f t 在区间]5,5[-中的最
小值点。

〖目的〗
● 理解极值概念的邻域性。

● 如何求最小值。

● 学习运用作图法求极值或最小值。

● 感受符号法的局限性。

〖解答〗
（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 ？');%该指令只能在指令窗中运行 t t=linspace(-5,5,N+1);
f or k=1:N
[tmin(k),fobj(k)]=fminbnd(ft,tt(k),tt(k+1));
e nd
[fobj,ii]=sort(fobj); %将目标值由小到大排列
t min=tmin(ii); %使极小值点做与目标值相应的重新排列 fobj,tmin
（2）最后确定的最小值点
在10,,2,1 =N 的不同分割下，经观察，最后确定出
最小值点是 -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)
grid on,shg
（B ）经观察后，把最小值附近邻域放到足够大，然后运行以下指令，那放大图形被推向前台，与此同时光标变为“十字线”，利用它点击极值点可得到最小值数据
[tmin2,fobj2]=ginput(1)
tmin2 =
-1.28500000993975
fobj2 =
-0.18604799369136
出现具有相同数值的刻度区域表明已达最小可分辨状态
（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 =
-.60100931947716486053884417850955e-2
fobj3 =
.89909908144684551670208797723124
〖说明〗
● 最小值是对整个区间而言的，极小值是对邻域而言的。

● 在一个区间中寻找最小值点，对不同子区间分割进行多次搜索是必要的。

这样可以避免
把极小值点误作为最小值点。

最小值点是从一系列极小值点和边界点的比较中确定的。

● 作图法求最小值点，很直观。

假若绘图时，自变量步长取得足够小，那么所求得的最小
值点有相当好的精度。

● 符号法在本例中，只求出一个极值点。

其余很多极值点无法秋初，更不可能得到最小值。

⏹ 6.设0)0(,1)0(,1)(2)(3)(22===+-dt
dy y t y dt t dy dt t y d ，用数值法和符号法求5.0)(=t t y 。

〖目的〗
● 学习如何把高阶微分方程写成一阶微分方程组。

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

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

〖解答〗
（1）改写高阶微分方程为一阶微分方程组
令dt
t dy t y t y t y )()(),()(21==，于是据高阶微分方程可写出 ⎪⎩⎪⎨⎧++-==1)(3)(2)()()(21221t y t y dt
t dy t y dt t dy
（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 =
.78958035647060552916850705213780
〖说明〗
●第<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。