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環(huán)境中，比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決定循環(huán)次數%計算各根的相角%計算各根

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]〖解答〗

前兩種輸入方法可以，后一種方法不行。

%整個指令在中文狀態(tài)下輸入

?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)')

〖目的〗

?理解自由符號變量的確認規(guī)則?！冀獯稹?/p>

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〖目的〗

?符號離散卷積直接法和變換法?！冀獯稹剑?）直接法

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，

作用下的輸出。

〖目的〗

?符號連續(xù)函數卷積的直接法和變換法。?符號變量限定性定義的作用。?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指令的應用?！冀獯稹?/p>

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的使用?！冀獯稹?/p>

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?〖目的〗

?符號計算用于定理證明?！冀獯稹?/p>

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變換表達式。

〖目的〗

?注意：變換中，被變換變量的約定?！冀獯稹?/p>

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指令中，被解方程的正確書寫，輸出量的正確次序?！冀獯稹?/p>

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

〖目的〗

?理解和把握信號流圖傳遞函數的“代數狀態(tài)方程解法〞。

?并設法用胡壽松主編的“自動控制原理〞的例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.采用代數狀態(tài)方程法求圖p2-2所示結構框圖的傳遞函數

Y。WY和U

圖p2-2

〖目的〗

?運用“代數狀態(tài)方程解法〞求輸入和擾動同時存在的結構框圖的傳遞函數。

〖解答〗

（1）理論演繹對于結構框圖寫出狀態(tài)方程

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）列出“元素級〞狀態(tài)方程值得提醒：在編寫M碼之前，最好先在草稿紙上，細心“元素級〞狀態(tài)方程是避免出錯的沖要措施。對此，不要掉以輕心。本例的“元素級〞狀態(tài)方程如下

?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的值空間中，方程沒有確鑿解?！冀獯稹剑?）借助增廣矩陣用指令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

〖提醒〗

?在適當范圍內，作圖觀測一元繁雜函數的形態(tài)：觀測解的存在性；解的唯一性。進而，

借助圖形法求近似解。?匿名函數的使用方法。?fzero指令的用法。

〖目的〗

?作圖法求一元繁雜函數解上的作用：觀測解的存在性；解的唯一性；得近似解。?匿名函數的使用方法。?fzero指令的用法?！冀獯稹?/p>

（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指令的用法?！冀獯稹剑?）符號法只能得到兩組解

?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);

可以看到大量解，并從中找到規(guī)律

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發(fā)生的次數數組%算出各發(fā)生次數的概率

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

49

0.120.10.080.060.040.020230406080100

?15試產生均值為4，標準差為2的(10000?1)的正態(tài)分布隨機數組

a，分別用hist和histfit繪制該數組的頻數直方圖，觀測兩張圖形的差異。除histfit上的擬合紅線外，你能使這兩個指令匯出一致的頻數直方圖嗎？

〖目的〗

?指令normrnd的應用，及隨機發(fā)生器的初始狀態(tài)控制。?指令hist與histfit的異同?！冀獯稹?/p>

（1）繪制頻數直方圖

rngdefault

%為本例結果可被讀者重現而設，并非必要。

mu=4;sigma=2;

%定義均值和標準差

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

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

%產生正態(tài)分布隨機數組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）經觀測后，把最小值附近鄰域放到足夠大，然后運行以下指令，那放大圖形被推向前臺，與此同時光標變?yōu)椤笆志€〞，利用它點擊極值點可得到最小值數據

[tmin2,fobj2]=ginput(1)

234

tmin2=

-1.28500000993975fobj2=

-0.18604799369136

出現具有一致數值的刻度區(qū)域說明已達最小可分辯狀態(tài)

（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

?最小值是對整個區(qū)間而言的，微小值是對鄰域而言的。?在一個區(qū)間中尋覓最小值點，對不同子區(qū)間分割進行屢屢探尋是必要的。這樣可以避免

把微小值點誤作為最小值點。最小值點是從一系列微小值點和邊界點的比較中確定的。?作圖法求最小值點，很直觀。假若繪圖時，自變量步長取得足夠小，那么所求得的最小

值點有相當好的精度。?符號法在本例中，只求出一個極值點。其余好多極值點無法秋初，更不可能得到最小值。

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

〖目的〗

?學習如何把高階微分方程寫成一階微分方程組。

?ode45解算器的導數函數如何采用匿名函數形式構成。

?如何從ode45一組數值解點，求指定自變量對應的函數值?！冀獯稹?/p>

（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節(jié)。

?8.已知由MATLAB指令創(chuàng)立的矩陣A=gallery(5)，試對該矩陣進

行特征值分解，并通過驗算觀測發(fā)生的現象。

〖目的〗

?展示特征值分解可能存在的數值問題。?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用于方程求解。

?矩陣除法能正確求得這類方程的特解。?逆陣法不能求得這類方程的特解。?解的驗證方法。?齊次解的獲取。?全解的獲得?！冀獯稹?/p>

（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