數(shù)學(xué)實(shí)驗(yàn)與建模：第5章 MATLAB數(shù)值計(jì)算

1、第5章MATLAB數(shù)值計(jì)算目錄目錄v在科學(xué)和工程應(yīng)用中，往往要進(jìn)行大量的數(shù)學(xué)計(jì)算。這些運(yùn)算一般來(lái)說(shuō)難以用手工精確和快捷地進(jìn)行，而要借助計(jì)算機(jī)編制相應(yīng)的程序做近似計(jì)算并不斷更新和擴(kuò)充。MATLAB的數(shù)值分析功能十分強(qiáng)大中，本章主要講述MAYLAB在函數(shù)、插值和曲線似合分析、微積分和線性方程系統(tǒng)方面的應(yīng)用。5.1 特殊矩陣特殊矩陣5.2 矩陣分析矩陣分析5.3 矩陣分解與線性方程組求解矩陣分解與線性方程組求解5.4 數(shù)據(jù)處理與多項(xiàng)式計(jì)算數(shù)據(jù)處理與多項(xiàng)式計(jì)算5.5 傅立葉分析傅立葉分析5.6 數(shù)值微積分?jǐn)?shù)值微積分5.7 常微分方程的數(shù)值求解常微分方程的數(shù)值求解5.8 非線性方程的數(shù)值求解非線性方程的

2、數(shù)值求解5.9 稀疏矩陣稀疏矩陣5.1 特殊矩陣特殊矩陣5.1.1對(duì)角陣與三角陣對(duì)角陣與三角陣1. 矩陣的對(duì)角元素矩陣的對(duì)角元素(1)提取矩陣的對(duì)角線元素提取矩陣的對(duì)角線元素 設(shè)設(shè)A為為mn矩陣，矩陣，diag(A)函數(shù)用于提取矩陣函數(shù)用于提取矩陣A主對(duì)角線主對(duì)角線元素產(chǎn)生一個(gè)具有元素產(chǎn)生一個(gè)具有min(m,n)個(gè)元素的列向量。個(gè)元素的列向量。 diag(A)函數(shù)還有更進(jìn)一步的形式函數(shù)還有更進(jìn)一步的形式diag(A,k)，其功能是提，其功能是提取第取第k條對(duì)角線的元素。條對(duì)角線的元素。目錄目錄A=17,0,1,0,15;23,5,7,14,16;4,0,13,0,22;10,12,19,21,

3、3;11,18,25,2,19A = 17 0 1 0 15 23 5 7 14 16 4 0 13 0 22 10 12 19 21 3 11 18 25 2 19 diag(A)ans = 17 5 13 21 19 diag(A,3)ans = 0 16(2)構(gòu)造對(duì)角矩陣構(gòu)造對(duì)角矩陣 設(shè)設(shè)V為具有為具有m個(gè)元素的向量，個(gè)元素的向量，diag(V)將產(chǎn)生一將產(chǎn)生一個(gè)個(gè)mm對(duì)角矩陣，其主對(duì)角線元素即為向量對(duì)角矩陣，其主對(duì)角線元素即為向量V的元素。的元素。 diag(V)函數(shù)也有更進(jìn)一步的形式函數(shù)也有更進(jìn)一步的形式diag(V,k)，其功能是產(chǎn)生一個(gè)其功能是產(chǎn)生一個(gè)nn(n=m+)對(duì)角陣，其第

4、對(duì)角陣，其第k條對(duì)角線的元素即為向量條對(duì)角線的元素即為向量V的元素。的元素。 V=1 2 3 4 5; diag(V)ans = 1 0 0 0 0 0 2 0 0 0 0 0 3 0 0 0 0 0 4 0 0 0 0 0 5diag(V,2)ans = 0 0 1 0 0 0 0 0 0 0 2 0 0 0 0 0 0 0 3 0 0 0 0 0 0 0 4 0 0 0 0 0 0 0 5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 例例5.1 先建立先建立55矩陣矩陣A，然后將，然后將A的第的第1行元素乘行元素乘以以1，第，第2行乘以行乘以2，第，第5行乘以行乘以5。目錄目錄

5、ans = 17 0 1 0 15 46 10 14 28 32 12 0 39 0 66 40 48 76 84 12 55 90 125 10 95A = 17 0 1 0 15 23 5 7 14 16 4 0 13 0 22 10 12 19 21 3 11 18 25 2 19命令如下：命令如下：A=17,0,1,0,15;23,5,7,14,16;4,0,13,0,22;10,12,19,21,3;11,18,25,2,19;D=diag(1,2,3,4,5);D*A 2. 矩陣的三角陣矩陣的三角陣 (1)下三角矩陣下三角矩陣 求矩陣求矩陣A的下三角陣的的下三角陣的MATLAB函數(shù)

6、是函數(shù)是tril(A) tril(A)函數(shù)也有更進(jìn)一步的一種形式函數(shù)也有更進(jìn)一步的一種形式tril(A,k)，其功能是求矩陣其功能是求矩陣A的第的第k條對(duì)角線以下的元素。條對(duì)角線以下的元素。 (2)上三角矩陣上三角矩陣 在在MATLAB中，提取矩陣中，提取矩陣A的上三角矩陣的函的上三角矩陣的函數(shù)是數(shù)是triu(A)和和triu(A,k)，其用法與提取下三角，其用法與提取下三角矩陣的函數(shù)矩陣的函數(shù)tril(A)和和tril(A,k)完全相同。完全相同。目錄目錄 tril(A)ans = 17 0 0 0 0 23 5 0 0 0 4 0 13 0 0 10 12 19 21 0 11 18 25

7、 2 19 triu(A)ans = 17 0 1 0 15 0 5 7 14 16 0 0 13 0 22 0 0 0 21 3 0 0 0 0 19 A=17,0,1,0,15;23,5,7,14,16;4,0,13,0,22;10,12,19,21,3;11,18,25,2,19A = 17 0 1 0 15 23 5 7 14 16 4 0 13 0 22 10 12 19 21 3 11 18 25 2 19 tril(A,1)ans = 17 0 0 0 0 23 5 7 0 0 4 0 13 0 0 10 12 19 21 3 11 18 25 2 19 5.1.2 特殊矩陣的生

8、成特殊矩陣的生成 1. 魔方矩陣魔方矩陣魔方矩陣是魔方矩陣是n*n元素所構(gòu)成的方陣，其每個(gè)元素元素所構(gòu)成的方陣，其每個(gè)元素由不同的由不同的1n2的整數(shù)所組成，它的每行、每列的整數(shù)所組成，它的每行、每列以及對(duì)角線元素之和均相等，并等于以及對(duì)角線元素之和均相等，并等于n(1+n2)/2.函數(shù)格式為函數(shù)格式為magic(n) 例例5.2 將將101125等等25個(gè)數(shù)填入一個(gè)個(gè)數(shù)填入一個(gè)5行行5列的表格中，列的表格中，使其每行每列及對(duì)角線的和均為使其每行每列及對(duì)角線的和均為565。 命令如下：命令如下： B=100+magic(5) B=100+magic(5)B = 117 124 101 108

9、115 123 105 107 114 116 104 106 113 120 122 110 112 119 121 103 111 118 125 102 109 sum(B(1,:)ans = 565 sum(B(2,:)ans = 565 sum(B(:,4)ans = 565 B(1,1)+B(2,2)+B(3,3)+B(4,4)+B(5,5)ans = 5652. 范得蒙矩陣范得蒙矩陣 函數(shù)函數(shù)vander(V)生成以向量生成以向量V為基礎(chǔ)向量的范得為基礎(chǔ)向量的范得蒙矩陣。蒙矩陣。VANDER Vandermonde matrix. A = VANDER(V) returns th

10、e Vandermonde matrix whose columns are powers of the vector V, that is A(i,j) = v(i)(n-j). p=1 2 3 4 5p = 1 2 3 4 5 A = VANDER(p)A = 1 1 1 1 1 16 8 4 2 1 81 27 9 3 1 256 64 16 4 1 625 125 25 5 1 3. 希爾伯特矩陣 Hilbert矩陣的每個(gè)元素的值，由行數(shù)i和列數(shù)j決定，等于1/(i+j-1)，生成希爾伯特矩陣的函數(shù)是hilb(n)。MATLAB中，有一個(gè)專門求希爾伯特矩陣的逆的函數(shù)invhilb(n)

11、，其功能是求n階的希爾伯特矩陣的逆矩陣。 hilb(4)ans = 1.0000 0.5000 0.3333 0.2500 0.5000 0.3333 0.2500 0.2000 0.3333 0.2500 0.2000 0.1667 0.2500 0.2000 0.1667 0.1429 I=invhilb(4)I = 16 -120 240 -140 -120 1200 -2700 1680 240 -2700 6480 -4200 -140 1680 -4200 28004. 托普利茲矩陣托普利茲矩陣 生成托普利茲矩陣的函數(shù)是生成托普利茲矩陣的函數(shù)是toeplitz(x,y)，它生成一個(gè)

12、，它生成一個(gè)以以x為第為第1列，列，y為第為第1行的托普利茲矩陣。這里行的托普利茲矩陣。這里x, y均為均為向量，二者不必等長(zhǎng)。當(dāng)向量，二者不必等長(zhǎng)。當(dāng)x和和y的第一個(gè)元素不同時(shí)，系的第一個(gè)元素不同時(shí)，系統(tǒng)將給出提示信息并以統(tǒng)將給出提示信息并以x中的元素為準(zhǔn)中的元素為準(zhǔn) c=1 2 3 4 5; r=1.5 2.5 3.5 4.5 5.5; toeplitz(c,r)Warning: First element of input column does not match first element of input row. Column wins diagonal conflict.(Ty

13、pe warning off MATLAB:toeplitz:DiagonalConflict to suppress this warning.) In E:matlabanzhuangtoolboxmatlabelmattoeplitz.m at line 18ans = 1.0000 2.5000 3.5000 4.5000 5.5000 2.0000 1.0000 2.5000 3.5000 4.5000 3.0000 2.0000 1.0000 2.5000 3.5000 4.0000 3.0000 2.0000 1.0000 2.5000 5.0000 4.0000 3.0000

14、2.0000 1.00005. 友矩陣 Compan矩陣生成多項(xiàng)式系數(shù)向量P的伴隨矩陣，其中（A(1,:)=-P(2:n)/P(1)）,友矩陣的函數(shù)是：compan(P) 。P是一個(gè)多項(xiàng)式的系數(shù)向量，高次冪系數(shù)排在前，低次冪排在后。例5.3求x3-7x+6的根。 u=1 0 -7 6u = 1 0 -7 6 A=compan(u)A = 0 7 -6 1 0 0 0 1 0 eig(compan(u)ans = -3.0000 2.0000 1.0000該多項(xiàng)式的根是該多項(xiàng)式的根是-3,2,1友矩陣的特征根正好是多項(xiàng)友矩陣的特征根正好是多項(xiàng)式的根式的根6. 帕斯卡矩陣帕斯卡矩陣 Pascal矩

15、陣是一個(gè)實(shí)對(duì)稱的正定矩陣，它由矩陣是一個(gè)實(shí)對(duì)稱的正定矩陣，它由Pascal三角形組三角形組成，成， Pascal三角形是由三角形是由0到到2n-1階的二項(xiàng)式系數(shù)組成，把二項(xiàng)階的二項(xiàng)式系數(shù)組成，把二項(xiàng)式系數(shù)依次填寫在矩陣的左側(cè)對(duì)角線上，提取左側(cè)的式系數(shù)依次填寫在矩陣的左側(cè)對(duì)角線上，提取左側(cè)的n行行n 列即列即為為Pascal矩陣。矩陣。函數(shù)函數(shù)pascal(n)生成一個(gè)生成一個(gè)n階的帕斯卡矩陣。階的帕斯卡矩陣。階數(shù)二項(xiàng)式系數(shù)0111 121 2 131 3 3 141 4 6 4 151 5 10 10 5 11 1 1 1 1 2 3 4 1 3 6 101 4 16 20 1 1 15 61

16、51 5 151 61Pascal(4)例5.4求(x+y)5的展開(kāi)式。在MATLAB命令窗口，輸入命令：pascal(6)ans = 1 1 1 1 1 1 1 2 3 4 5 6 1 3 6 10 15 21 1 4 10 20 35 56 1 5 15 35 70 126 1 6 21 56 126 252其其次對(duì)角線次對(duì)角線上的元素上的元素1，5，10，10，5，1即為展開(kāi)即為展開(kāi)式的系數(shù)。式的系數(shù)。5.2 矩陣分析5.2.1 矩陣結(jié)構(gòu)變換1. 矩陣的轉(zhuǎn)置轉(zhuǎn)置運(yùn)算符是單撇號(hào)()。2. 矩陣的旋轉(zhuǎn) 矩陣的旋轉(zhuǎn)利用函數(shù)rot90(A,k)，功能是將矩陣A旋轉(zhuǎn)90的k倍，當(dāng)k為1時(shí)可省略。3

17、. 矩陣的左右翻轉(zhuǎn)對(duì)矩陣A實(shí)施左右翻轉(zhuǎn)的函數(shù)是fliplr(A)。4. 矩陣的上下翻轉(zhuǎn)對(duì)矩陣A實(shí)施上下翻轉(zhuǎn)的函數(shù)是flipud(A)。目錄目錄A=magic(3)A = 8 1 6 3 5 7 4 9 2 fliplr(A) %將矩陣的列左右翻轉(zhuǎn)將矩陣的列左右翻轉(zhuǎn)ans = 6 1 8 7 5 3 2 9 4 flipud(A)%將矩陣將矩陣A的行上、的行上、下翻轉(zhuǎn)下翻轉(zhuǎn)ans = 4 9 2 3 5 7 8 1 6 rot90(A)%將矩陣將矩陣A逆時(shí)針旋轉(zhuǎn)逆時(shí)針旋轉(zhuǎn)90度度ans = 6 7 2 1 5 9 8 3 4 5.2.2 矩陣的逆與偽逆1. 矩陣的逆 求一個(gè)矩陣的逆非常容易。求方

18、陣A的逆可調(diào)用函數(shù)inv(A)。例5.4 用求逆矩陣的方法解線性方程組。命令如下：一般情況下，用左除比求矩陣的逆的方法更有效，即x=Ab。目錄目錄 A=1,2,3;1,4,9;1,8,27; b=5,-2,6; x=inv(A)*bx = 23.0000 -14.5000 3.6667 x=Abx = 23.0000 -14.5000 3.6667 2. 矩陣的偽逆MATLAB中，求一個(gè)矩陣偽逆的函數(shù)是pinv(A)。例5.5 求A的偽逆，并將結(jié)果送B。命令如下：A=3,1,1,1;1,3,1,1;1,1,3,1;B=pinv(A)例5.6 求矩陣A的偽逆。在MATLAB命令窗口，輸入命令：A

19、=0,0,0;0,1,0;0,0,1;pinv(A)目錄目錄B = 0.3929 -0.1071 -0.1071 -0.1071 0.3929 -0.1071 -0.1071 -0.1071 0.3929 0.0357 0.0357 0.0357ans = 0 0 0 0 1 0 0 0 1 5.2.3 方陣的行列式方陣的行列式求方陣求方陣A所對(duì)應(yīng)的行列式的值的函數(shù)是所對(duì)應(yīng)的行列式的值的函數(shù)是det(A)。例例5.7用克萊姆用克萊姆(Cramer)方法求解線性方程組。方法求解線性方程組。程序如下：程序如下：D=2,2,-1,1;4,3,-1,2;8,5,-3,4;3,3,-2,2; %定義系數(shù)

20、矩陣定義系數(shù)矩陣b=4;6;12;6; %定義常數(shù)項(xiàng)向量定義常數(shù)項(xiàng)向量D1=b,D(:,2:4); %用方程組的右端向量置換用方程組的右端向量置換D的第的第1列列D2=D(:,1:1),b,D(:,3:4); %用方程組的右端向量置換用方程組的右端向量置換D的第的第2列列D3=D(:,1:2),b,D(:,4:4); %用方程組的右端向量置換用方程組的右端向量置換D的第的第3列列D4=D(:,1:3),b; %用方程組的右端向量置換用方程組的右端向量置換D的第的第4列列DD=det(D);x1=det(D1)/DD;x2=det(D2)/DD;x3=det(D3)/DD;x4=det(D4)/

21、DD;x1,x2,x3,x4目錄目錄ans = 1 1 -1 -1 5.2.4 矩陣的秩MATLAB中，求矩陣秩的函數(shù)是rank(A)。例如，求例5.7中方程組系數(shù)矩陣D的秩，命令是：說(shuō)明D是一個(gè)滿秩矩陣。目錄目錄 D=2,2,-1,1;4,3,-1,2;8,5,-3,4;3,3,-2,2;r=rank(D)r = 4 5.2.5 向量和矩陣的范數(shù)(不講)1. 計(jì)算向量3種常用范數(shù)的函數(shù)(1)norm(V)或norm(V,2) 計(jì)算向量V的2范數(shù)=sum(abs(V).2)(1/2)(2)norm(V,1)=sum(abs(V) 計(jì)算向量V的1范數(shù)(3)norm(V,inf) 計(jì)算向量V的范數(shù)

22、= max(abs(V)目錄目錄例5.8 已知V，求V的3種范數(shù)。 v1=norm(V,1) %求V的1范數(shù)sum(abs(V) v1 = 5/2 v2=norm(V) %求V的2范數(shù)sum(abs(V).2)(1/2)v2 = 3/2 v3=norm(V,inf) %求V的范數(shù)max(abs(V)v3 = 1 2. 矩陣的范數(shù)及其計(jì)算函數(shù)MATLAB中提供了求3種矩陣范數(shù)的函數(shù)，其函數(shù)調(diào)用格式與求向量的范數(shù)的函數(shù)完全相同NORM(X) is the largest singular value of X, max(svd(X).NORM(X,2) is the same as NORM(X)

23、.NORM(X,1) is the 1-norm of X, the largest column sum, = max(sum(abs(X).NORM(X,inf) is the infinity norm of X, the largest row sum, = max(sum(abs(X).A=17,0,1,0,15;23,5,7,14,16;4,0,13,0,22;10,12,19,21,3;11,18,25,2,19;a1=norm(A,1) %求A的1范數(shù)a2=norm(A) %求A的2范數(shù)ainf=norm(A,inf) %求A的范數(shù)例5.9 求矩陣A的三種范數(shù)。命令如下：a1

24、= 75 a2 = 2790/47 ainf = 75 5.2.6 矩陣的條件數(shù)和跡1. 的條件數(shù)MATLAB中，計(jì)算矩陣A的3種條件數(shù)的函數(shù)是：(1)cond(A,1) 計(jì)算A的1范數(shù)下的條件數(shù)(2)cond(A)或cond(A,2) 計(jì)算A的2范數(shù)數(shù)下的條件數(shù)(3)cond(A,inf) 計(jì)算A的 范數(shù)下的條件數(shù)目錄目錄例5.10 求矩陣X的三種條件數(shù)。命令如下：A=2,2,3;4,5,-6;7,8,9;C1=cond(A,1)C2=cond(A)C3=cond(A,inf)C1 = 1044/7 C2 = 7126/81 C3 = 144 2. 矩陣的跡MATLAB中，求矩陣的跡的函數(shù)是

25、trace(A)。例如，X=2 2 3;4 5 -6;7 8 9;trace(X)ans = 16目錄目錄5.2.7 矩陣的特征值與特征向量MATLAB中，計(jì)算矩陣A的特征值和特征向量的函數(shù)是eig(A)，常用的調(diào)用格式有3種：(1)E=eig(A) 求矩陣A的全部特征值，構(gòu)成向量E。(2)V,D=eig(A) 求矩陣A的全部特征值，構(gòu)成對(duì)角陣D，并求A的特征向量構(gòu)成V的列向量。(3)V,D=eig(A,nobalance) 與第2種格式類似，但第2種格式中先對(duì)A作相似變換后求矩陣A的特征值和特征向量，而格式3直接求矩陣A的特征值和特征向量。例5.11 用3種不同的格式求A的特征值和特征向量。

26、命令如下：A=1,2,2;1,-1,1;4,-12,1;E=eig(A)V,D=eig(A)V,D=eig(A,nobalance)目錄目錄例5.12用求特征值的方法解方程3x5-7x4+5x2+2x-18。命令如下：p=3,-7,0,5,2,-18;A=compan(p); %A的友矩陣x1=eig(A) %求A的特征值x2=roots(p) %直接求多項(xiàng)式p的零點(diǎn)兩種方法求得的方程的根是完全一致的，實(shí)際上，roots函數(shù)正是應(yīng)用求友矩陣的特征值的方法來(lái)求方程的根。x1 = 2.1837 1.0000 + 1.0000i 1.0000 - 1.0000i -0.9252 + 0.7197i

27、-0.9252 - 0.7197ix2 = 2.1837 1.0000 + 1.0000i 1.0000 - 1.0000i -0.9252 + 0.7197i -0.9252 - 0.7197i5.2.8 MATLAB在三維向量中的應(yīng)用1. 向量共線或共面的判斷例5.13 設(shè)X=(1，1，1)，Y=(-1，2，1)，Z=(2，2，2)，判斷這三個(gè)向量的共線共面問(wèn)題。命令如下：X=1,1,1;Y=-1,2,1;Z=2,2,2;XY=X;Y;YZ=Y;Z;ZX=Z;X;XYZ=X;Y;Z;rank(XY)rank(YZ)rank(ZX)rank(XYZ)目錄目錄ans = 2ans = 2ans

28、 = 1ans = 2 2. 向量方向余弦的計(jì)算例5.14設(shè)向量V=(5，-3，2)，求V的方向余弦。建立一個(gè)函數(shù)文件direct.m：function f=f(v)r=norm(v);if r=0 f=0else f=v(1)/r,v(2)/r,v(3)/r;endreturn在MATLAB命令窗口，輸入命令：v=5,-3,2;f=direct(v)目錄目錄f = 2220/2737 -1332/2737 888/2737 3. 向量的夾角例5.15 設(shè)U=(1，0，0)，V=(0，1，0)，求U,V間的夾角.命令如下：U=1,0,0;V=0,1,0;r1=norm(U);r2=norm(V

29、);UV=U*V;cosd=UV/r1/r2;D=acos(cosd)4. 兩點(diǎn)間的距離例5.16 設(shè) U=(1，0，0)，V=(0，1，0)，求U、V兩點(diǎn)間的距離。命令如下：U=1,0,0;V=0,1,0;UV=U-V;D=norm(UV)目錄目錄 5. 向量的向量積例5.17設(shè)U=(2，-3，1)，V=(3，0，4)，求UV。命令如下：U=2,-3,1;V=3,0,4;W=eye(3);A1=W(1,:);U;V;A2=W(2,:);U;V;A3=W(3,:);U;V;UV=det(A1),det(A2),det(A3)UV= -12 -5 96. 向量的混合積例5.18 設(shè)U=(0，0，

30、2)，V=(3，0，5)，W=(1，1，0)，求以這三個(gè)向量構(gòu)成的六面體的體積。命令如下：U=0,0,2;V=3,0,5;W=1,1,0;A=U;V;W;det(A)ans = 6目錄目錄7. 點(diǎn)到平面的距離例5.19求原點(diǎn)到平面X+Y+Z=1的距離。命令如下：u=0,0,0;v=1,1,1; % A=B=C=1,u1=u2=u3=0,D=-1r=abs(u*v-1)/norm(v,2)r =0.5774目錄目錄5.3 矩陣分解與線性方程組求解矩陣分解與線性方程組求解5.3.1矩陣分解1. 實(shí)對(duì)稱矩陣的QDQ分解例5.20設(shè)對(duì)稱矩陣A，對(duì)A進(jìn)行QDQ分解。命令如下：A=2,1,4,6;1,2,

31、1,5;4,1,3,4;6,5,4,2;Q,D=eig(A)Q*D*Qans = 2.0000 1.0000 4.0000 6.0000 1.0000 2.0000 1.0000 5.0000 4.0000 1.0000 3.0000 4.0000 6.0000 5.0000 4.0000 2.0000結(jié)果與A相等，說(shuō)明確實(shí)將A分解為了QDQ的乘積。目錄目錄 例5.21求下列二次型的標(biāo)準(zhǔn)形式及變換矩陣。命令如下：A=1,2,1;2,1,1;1,1,3;Q,D=eig(A)進(jìn)一步作線性變換即得關(guān)于u,v,w的標(biāo)準(zhǔn)二次型：2. 矩陣的LU分解MATLAB中，完成LU分解的函數(shù)是：(1)L,U=lu

32、(A) 將方陣A分解為交換下三角矩陣L和上三角矩陣U，使 A=LU。(2)L,U,P=lu(A) 將方陣A分解為下三角矩陣L和上三角矩陣U，使 PA=LU。LU分解常用于求行列式以及解線性方程組。目錄目錄3. 矩陣的QR分解（正交分解）對(duì)矩陣A進(jìn)行QR分解的函數(shù)是Q,R=qr(A)，根據(jù)方陣A，求一個(gè)正交矩陣Q和一個(gè)上三角矩陣R，使A=Q*R。例如，對(duì)矩陣A進(jìn)行QR分解的命令是：A=2,1,-2;1,2,1;2,5,3;Q,R=qr(A)目錄目錄 5.3.2 線性方程組求解線性方程組求解1. 線性方程組解的一般討論線性方程組解的一般討論解線性方程組的一般解線性方程組的一般函數(shù)文件函數(shù)文件如下：

33、如下：function x,y=line_solution(A,b) m,n=size(A);y=; if norm(b)0 %非齊次方程組 if rank(A)=rank(A,b) %方程組相容 if rank(A)=m %有唯一解 x=Ab; else %方程組有無(wú)窮多個(gè)解,基礎(chǔ)解系 disp(原方程組有有無(wú)窮個(gè)解，其齊次方程組的基礎(chǔ)解系為y，特解為x); y=null(A,r); x=Ab; end else %方程組不相容，給出最小二乘法解 disp(方程組的最小二乘法解是：); x=Ab; 目錄目錄end else %齊次方程組 if rank(A)=n %列滿秩 x=zero(m,

34、1) %0解 else %非0解 disp(方程組有無(wú)窮個(gè)解，基礎(chǔ)解系為x); x=null(A,r); end endreturn2. 應(yīng)用舉例例5.23求線性方程組的解。在MATLAB命令窗口，輸入命令：A=2,2,-1,1;4,3,-1,2;8,5,-3,4;3,3,-2,2;b=4,6,12,6;x,y=line_solution(A,b) %調(diào)用自定義函數(shù)目錄目錄x = 1 1 -1 -1 y = 例例5.24求下列線性方程組的解。在MATLAB命令窗口，輸入命令：A=2,7,3,1;3,5,2,2;9,4,1,7;b=6,4,2;x,y=line_solution(A,b)原方程組

35、有有無(wú)窮個(gè)解，其齊次方程組的基礎(chǔ)解系為y，特解為xWarning: Rank deficient, rank = 2 tol = 8.6112e-015. In E:matlabanzhuangworkline_solution.m at line 10 x = -2/11 10/11 0 0 y = 1/11 -9/11 -5/11 1/11 1 0 0 1 5.4 數(shù)據(jù)處理與多項(xiàng)式計(jì)算數(shù)據(jù)處理與多項(xiàng)式計(jì)算5.4.1 數(shù)據(jù)統(tǒng)計(jì)與分析1. 求矩陣最大和最小元素(1)求向量的最大最小元素y=max(X) 返回向量X的最大元素存入y。y,I=max(X) 返回向量X的最大元素存入y，最大元素的序號(hào)

36、存入I。(2)求矩陣的最大和最小元素max(A) 返回一個(gè)行向量，向量的第i個(gè)元素是A矩陣的第i列上的最大元素。Y,U=max(A) 返回兩個(gè)行向量，Y向量記錄A的每列的最大元素，U向量記錄每列最大元素的行號(hào)。max(A,dim) dim取1或2。dim取1時(shí)，該函數(shù)和max(A)完全相同。dim取2時(shí)，該函數(shù)返回一個(gè)列向量，其第i個(gè)元素是A矩陣的第i行上的最大元素。目錄目錄(3)兩個(gè)向量或矩陣對(duì)應(yīng)元素的比較U=max(A,B) A,B是兩個(gè)同型的向量或矩陣。結(jié)果U是與A,B同型的向量或矩陣，U的每個(gè)元素等于A,B對(duì)應(yīng)元素的較大者。U=max(A,n) n是一個(gè)標(biāo)量。結(jié)果U是與A同型的向量或矩

37、陣，U的每個(gè)元素等于A對(duì)應(yīng)元素和n中的較大者。min函數(shù)的用法和max完全相同。目錄目錄例5.25 求矩陣A的每行及每列的最大和最小元素，并求整個(gè)矩陣的最大和最小元。命令如下：A=13,-56,78;25,63,-235;78,25,563;1,0,-1;max(A,2) %求每行最大元素min(A,2) %求每行最小元素max(A) %求每列最大元素min(A) %求每列最小元素max(max(A) %求整個(gè)矩陣的最大元素min(min(A) %求整個(gè)矩陣的最小元素目錄目錄ans = 78 63 563 1 ans = -56 -235 25 -1 ans = 78 63 563 ans =

38、 1 -56 -235 ans = 563 ans = -235 2. 求矩陣的平均值和中值 求矩陣和向量元素的平均值的函數(shù)是mean，求中值的函數(shù)是median。它們的調(diào)用方法和max函數(shù)完全相同。3. 矩陣元素求和與求積矩陣和向量求和與求積的基本函數(shù)是sum和prod，其使用方法和max類似。目錄目錄例5.26求矩陣A的每行元素的乘積和全部元素的乘積。命令如下：A=1,2,3,4;5,6,7,8;9,10,11,12;S=prod(A,2)prod(S) %求A的全部元素的乘積目錄目錄4. 矩陣元素累加和與累乘積MATLAB中，使用cumsum和cumprod函數(shù)能方便地求得向量和矩陣元素

39、的累加和與累乘積向量，函數(shù)的用法和sum及prod相同例5.27求向量X=(1！，2！，3！，10！)。命令如下：Format long;X=cumprod(1:10)X = 1 2 6 24 120 720 5040 40320 362880 36288005. 標(biāo)準(zhǔn)方差 MATLAB中，提供了計(jì)算數(shù)據(jù)序列的標(biāo)準(zhǔn)方差的函數(shù)std。對(duì)于向量X，std(X)返回一個(gè)標(biāo)準(zhǔn)方差。對(duì)于矩陣A，std(A)返回一個(gè)行向量，它的各個(gè)元素便是矩陣A各列或各行的標(biāo)準(zhǔn)方差。std函數(shù)的一般調(diào)用格式為：std(A,FLAG,dim) 其中dim取1或2。 當(dāng)dim=1時(shí)，求各列元素的標(biāo)準(zhǔn)方差； 當(dāng)dim=2時(shí)，則

40、求各行元素的標(biāo)準(zhǔn)方差。 FLAG取0或1。目錄目錄6. 元素排序 MATLAB中對(duì)向量X是排序函數(shù)是sort(X)，函數(shù)返回一個(gè)對(duì)X中的元素按升序排列的新向量。 sort函數(shù)也可以對(duì)矩陣A的各列(或行)重新排序，其調(diào)用格式為：Y,I=sort(A,dim) 其中dim指明對(duì)A的列還是行進(jìn)行排序，若dim=1，則按列排，若dim=2，則按行排。Y是排序后的矩陣，而I記錄Y中的元素在A中位置。目錄目錄例5.28對(duì)矩陣做各種排序。命令如下：A=1,-8,5;4,12,6;13,7,-13;sort(A) %對(duì)A的每列按升序排序-sort(-A,2) %對(duì)A的每行按降序排序X,I=sort(A) %對(duì)

41、A按列排序，并將每個(gè)元素所在行號(hào)送矩陣I目錄目錄ans = 1 -8 -13 4 7 5 13 12 6ans = 5 1 -8 12 6 4 13 7 -13X = 1 -8 -13 4 7 5 13 12 6I = 1 1 3 2 3 1 3 2 25.4.2 數(shù)值插值數(shù)值插值1. 一維數(shù)值插值一維數(shù)值插值 interp1函數(shù)調(diào)用格式為：Y1=interp1(X,Y,X1,method) 函數(shù)根據(jù)X、Y的值，計(jì)算函數(shù)在X1處的值。X、Y是兩個(gè)等長(zhǎng)的已知向量，分別描述采樣點(diǎn)和樣本值，X1是一個(gè)向量或標(biāo)量，描述欲插值的點(diǎn)，Y1是一個(gè)與X1等長(zhǎng)的插值結(jié)果。 method是插值方法，允許的取值有l(wèi)

42、inear(線性插值)、nearest(最近插值)、spline(三次樣條插值)、cubic（三次多項(xiàng)式插值），缺省值是linear。目錄目錄例例5.29用不同的插值方法計(jì)算sin(x)在/2點(diǎn)的值。這是一個(gè)一維插值問(wèn)題。在MATLAB命令窗口，輸入命令：X=0:0.2:pi;Y=sin(X); %給出X、Yinterp1(X,Y,pi/2) %用缺省方法(即線性插值方法)計(jì)算sin(/2)interp1(X,Y,pi/2,nearest) %用最近方法計(jì)算sin(/2)interp1(X,Y,pi/2,linear) %用線性方法計(jì)算sin(/2)interp1(X,Y,pi/2,splin

43、e) %用三次樣條方法計(jì)算sin(/2)interp1(X,Y,pi/2,cubic) %用三次多項(xiàng)式方法計(jì)算sin(/2) MATLAB中有一個(gè)專門的三次樣條插值函數(shù)Y1=spline(X,Y,X1)，其功能及使用方法與函數(shù)Y1=interp1(X,Y,X1,spline)完全相同。目錄目錄 例5.30 已知檢測(cè)參數(shù)f隨時(shí)間t的采樣結(jié)果，用數(shù)值插值法計(jì)算t=2，7，12，17，22，17，32，37，42，47，52，57時(shí)f的值。這是一個(gè)一維數(shù)值插值問(wèn)題，命令如下：T=0:5:65;X=2:5:57;F=3.2015,2.2560,879.5,1835.9,2968.8,4136.2,52

44、37.9,6152.7,6725.3,6848.3,6403.5,6824.7,7328.5,7857.6;F1=interp1(T,F,X) %用線性方法插值F1=interp1(T,F,X,nearest) %用最近方法插值F1=interp1(T,F,X,spline) %用三次樣條方法插值 F1=interp1(T,F,X,cubic) %用三次多項(xiàng)式方法插值目錄目錄2. 二維數(shù)值插值 MATLAB中，提供了解決二維插值問(wèn)題的函數(shù)。其調(diào)用格式為：Z1=interp2(X,Y,Z,X1,Y1,method) 其中X、Y是兩個(gè)向量，分別描述兩個(gè)參數(shù)的采樣點(diǎn)，Z是與參數(shù)采樣點(diǎn)對(duì)應(yīng)的采樣變量的

45、樣本值，X1、Y1是兩個(gè)向量或標(biāo)量，描述欲插值的點(diǎn)。method的取值與一維插值函數(shù)相同。目錄目錄例5.31設(shè)Z=x2+y2，對(duì)Z函數(shù)在(0，1)(0，2)區(qū)域內(nèi)進(jìn)行插值。命令如下：x=0:0.1:10;y=0:0.2:20;X,Y=meshgrid(x,y);Z=X.2+Y.2;interp2(x,y,Z,0.5,0.5) %對(duì)函數(shù)在(0.5,0.5)點(diǎn)進(jìn)行插值interp2(x,y,Z,0.5 0.6,0.4) %對(duì)函數(shù)在(0.5,0.4)點(diǎn)和(0.6，0.4)點(diǎn)進(jìn)行插值interp2(x,y,Z,0.5 0.6,0.4 0.5) %對(duì)函數(shù)在(0.5,0.4)點(diǎn)和(0.6，0.5)點(diǎn)進(jìn)行插

46、值interp2(x,y,Z,0.5 0.6,0.4 0.5)%對(duì)函數(shù)在(0.5,0.4),(0.6,0.4),(0.5,0.5)和(0.6,0.5)點(diǎn)進(jìn)行插值目錄目錄ans = 0.5100ans = 0.4100 0.5200ans = 0.4100 0.6200ans = 0.4100 0.5200 0.5100 0.62003. 三維數(shù)值插值 對(duì)三維函數(shù)插值的函數(shù)是interp3，其使用方法和interp2相同。其調(diào)用格式為：W1=interp3(X,Y,Z,W,X1,Y1,Z1,method) 函數(shù)返回三維插值結(jié)果。其中X、Y、Z是三個(gè)向量，分別描述三個(gè)參數(shù)的采樣點(diǎn)，W是與參數(shù)采樣點(diǎn)

47、對(duì)應(yīng)的采樣變量的樣本值，X1、Y1、Z1是三個(gè)向量或標(biāo)量，描述欲插值的點(diǎn)。method是插值方法，可選，其缺省值是 line。method的取值與一、二維插值函數(shù)相同。目錄目錄5.4.3 曲線擬合 MATLAB中，提供了解決使用最小二乘法進(jìn)行曲線擬合的函數(shù)。調(diào)用格式為：P,S=polyfit(X,Y,m) 函數(shù)根據(jù)采樣點(diǎn)X和采樣點(diǎn)函數(shù)值Y，產(chǎn)生一個(gè)m次多項(xiàng)式P及其在采樣點(diǎn)的誤差向量S。 其中X、Y是兩個(gè)等長(zhǎng)的向量，P是一個(gè)長(zhǎng)度為m+1的向量。目錄目錄例5.32 用一個(gè)5次多項(xiàng)式在區(qū)間0，2內(nèi)逼近函數(shù)sin(x)。命令如下：目錄目錄 X=linspace(0,2*pi,50);Y=sin(X);

48、P,S=polyfit(X,Y,5) %得到5次多項(xiàng)式的系數(shù)和誤差plot(X,Y,k*,X,polyval(P,X),k-)P = -65/11683 287/3284 -58/147 578/2153 1528/1737 143/13993 S = R: 6x6 double df: 44 normr: 132/3913 5.4.4 多項(xiàng)式計(jì)算1. 多項(xiàng)式的建立 已知一個(gè)多項(xiàng)式的全部根X求多項(xiàng)式系數(shù)的函數(shù)是poly(X)，該函數(shù)返回以X為全部根的一個(gè)多項(xiàng)式P，當(dāng)X是一個(gè)長(zhǎng)度為m的向量時(shí)，P是一個(gè)長(zhǎng)度為m+1的向量。目錄目錄求多項(xiàng)式例已知多項(xiàng)式的根為,34,34, 3, 2, 1jj R=-

49、1 -2 -3 -4+3*j -4-3*j; %輸入多項(xiàng)式根的向量輸入多項(xiàng)式根的向量 P=poly(R)%構(gòu)造多項(xiàng)式的系數(shù)向量構(gòu)造多項(xiàng)式的系數(shù)向量P = 1 14 84 244 323 1502. 多項(xiàng)式求根多項(xiàng)式求根求多項(xiàng)式求多項(xiàng)式p(x)的根的函數(shù)是的根的函數(shù)是roots(P)，這里，這里，P是是p(x)的系數(shù)向量，該函數(shù)返回方程的系數(shù)向量，該函數(shù)返回方程p(x)=0的全部根的全部根(含含重根，復(fù)根重根，復(fù)根)。例：已知例：已知5階多項(xiàng)式的系數(shù)向量為階多項(xiàng)式的系數(shù)向量為P1828586730，求多項(xiàng)式的根求多項(xiàng)式的根 P=1 8 28 58 67 30; R=roots(P)R = -1.

50、00000000000000 + 2.00000000000000i -1.00000000000000 - 2.00000000000000i -2.99999999999998 -2.00000000000001 -1.00000000000000 3. 多項(xiàng)式求值多項(xiàng)式求值求多項(xiàng)式求多項(xiàng)式p(x)在某點(diǎn)或某些點(diǎn)的函數(shù)值的函數(shù)是在某點(diǎn)或某些點(diǎn)的函數(shù)值的函數(shù)是polyval(P,x)。若。若x為一數(shù)值，則求多項(xiàng)式在該點(diǎn)的為一數(shù)值，則求多項(xiàng)式在該點(diǎn)的值；若值；若x為向量或矩陣，則對(duì)向量或矩陣中的每個(gè)元為向量或矩陣，則對(duì)向量或矩陣中的每個(gè)元素求其多項(xiàng)式的值。素求其多項(xiàng)式的值。例例5.33 已知一

51、個(gè)多項(xiàng)式已知一個(gè)多項(xiàng)式3x5+4x3-5x2-7.2x+5，計(jì)算：，計(jì)算： (1)計(jì)算計(jì)算f(x)=0 的全部根。的全部根。(2)由方程由方程f(x)=0的根構(gòu)造一個(gè)多項(xiàng)式的根構(gòu)造一個(gè)多項(xiàng)式g(x)，并與，并與f(x)進(jìn)行對(duì)比。進(jìn)行對(duì)比。(3)計(jì)算計(jì)算f(5)、f(7.8)、f(9.6)、f(12.3)的值。的值。命令如下：命令如下：P=3,0,4,-5,-7.2,5;X=roots(P) %求方程求方程f(x)=0的根的根G=poly(X) %求多項(xiàng)式求多項(xiàng)式g(x)X0=5,7.8,9.6,12.3;f=polyval(P,X0) %求多項(xiàng)式求多項(xiàng)式f(x)在給定點(diǎn)的值在給定點(diǎn)的值 多項(xiàng)式

52、求值還有一個(gè)函數(shù)是多項(xiàng)式求值還有一個(gè)函數(shù)是polyvalm，其調(diào)用格式與，其調(diào)用格式與polyval相相同，但含義不同。同，但含義不同。polyvalm函數(shù)要求函數(shù)要求x為方陣，它以方陣為自變量求為方陣，它以方陣為自變量求多項(xiàng)式的值。多項(xiàng)式的值。目錄目錄X = -0.30455728045903 + 1.62172723875407i -0.30455728045903 - 1.62172723875407i -1.00664150535070 1.01901848386455 0.59673758240420 G = 1.00000000000000 0.00000000000000 1.3

53、3333333333333 -1.66666666666667 -2.39999999999999 1.66666666666666f = 1.0e+005 * 0.09719000000000 0.88158079040000 2.47625833280000 8.51195163290000 4. 多項(xiàng)式的四則運(yùn)算(1)多項(xiàng)式的加減法(2)多項(xiàng)式的乘法函數(shù)conv(P1,P2)用于求多項(xiàng)式P1和P2的乘積。(3)多項(xiàng)式的除法函數(shù)Q,r=deconv(P1,P2)用于對(duì)多項(xiàng)式P1和P2作除法運(yùn)算。其中Q返回多項(xiàng)式P1除以P2的商式，r返回P1除以P2的余式。這里，Q和r仍是多項(xiàng)式系數(shù)向量。d

54、econv是conv的逆函數(shù)，即有P1=conv(P2,Q)+r。目錄目錄例5.34設(shè)有兩個(gè)多項(xiàng)式，計(jì)算：(1)求f(x)+g(x)、f(x)-g(x)。(2)求f(x)g(x)、f(x)/g(x)。目錄目錄373)(657253)(22345xxxgxxxxxxf在在MATLAB命令窗口，輸入命令命令窗口，輸入命令：f=3,-5,2,-7,5,6;g=3,5,-3;g1=0,0,0,g;f+g1 %求f(x)+g(x)f-g1 %求f(x)-g(x)conv(f,g) %求f(x)*g(x)Q,r=deconv(f,g) %求f(x)/g(x)，商式送Q，余式送r。5. 多項(xiàng)式的導(dǎo)函數(shù)對(duì)多項(xiàng)

55、式求導(dǎo)數(shù)的函數(shù)是：p=polyder(P) 求多項(xiàng)式P的導(dǎo)函數(shù)p=polyder(P,Q) 求P*Q的導(dǎo)函數(shù)p,q=polyder(P,Q) 求P/Q的導(dǎo)函數(shù)，導(dǎo)函數(shù)的分子存入p，分母存入q。目錄目錄例5.35求有理分式 的導(dǎo)數(shù)。命令如下：目錄目錄100765105853236910245xxxxxxxxxP=3,5,0,-8,1,-5;Q=10,5,0,0,6,0,0,7,-1,0,-100;p,q=polyder(P,Q)5.4.5 函數(shù)的最大值與最小值MATLAB中用于求最小值的函數(shù)是：fmin(f,a,b) 求單變量函數(shù)f(x)在區(qū)間(a,b)上的最小值點(diǎn)。fmins(F,X0) 求多

56、變量函數(shù)F(x)在估計(jì)值X0附近的最小值點(diǎn)。MATLAB沒(méi)有專門提供求函數(shù)最大值點(diǎn)的函數(shù)，但只要注意到-f(x)在區(qū)間(a,b)上的最小值點(diǎn)就是f(x)在(a,b)的最大值點(diǎn)，所以fmin(-f,a,b)返回函數(shù)f(x)在區(qū)間(a,b)上的最大值。目錄目錄例5.36 求函數(shù)f(x)在區(qū)間(-10,-1)和(1,10)上的最小值點(diǎn)。首先建立函數(shù)文件fx.m：function f=f(x)f=x-1/x+5;return再在MATLAB命令窗口，輸入命令：fmin(fx,-10,-1) %求函數(shù)在區(qū)間(-10，-1)內(nèi)的最小值點(diǎn)fmin(f,1,10) %求函數(shù)在區(qū)間(1，10)內(nèi)的最小值點(diǎn)。注意

57、函數(shù)名f不用加目錄目錄例5.37 設(shè)有函數(shù)f(x,y,z)，求函數(shù)f在(0.5,0.5,0.5)附近的最小值。建立函數(shù)文件fxyz.m：function f=f(u)x=u(1);y=u(2);z=u(3);f=x+y.2./x/4+z.2./y+2./z;return在MALAB命令窗口，輸入命令：U=fmins(fxyz,0.5,0.5,0.5) %求函數(shù)的最小值點(diǎn)fxyz(U) %求函數(shù)的最小值目錄目錄5.5 傅立葉分析(不講)MATLAB中，提供了對(duì)向量(或直接對(duì)矩陣的行或列)進(jìn)行離散傅立葉變換的函數(shù)，其調(diào)用格式是：Y=fft(X,n,dim)(1)當(dāng)X是一個(gè)向量時(shí)，返回對(duì)X的離散傅立

58、葉變換。(2)當(dāng)X是一個(gè)矩陣時(shí)，返回一個(gè)矩陣并送Y，其列(行)是對(duì)X的列(行)的離散傅立葉變換。目錄目錄例5.38 求X=(1，0，-3，5，2)的離散傅立葉逆變換。在MATLAB命令窗口，輸入命令：X=1,0,-3,5,2;Y=fft(X) %對(duì)X進(jìn)行變換3. 離散傅立葉變換的逆變換MATLAB中，對(duì)向量(或直接對(duì)矩陣的行或列)進(jìn)行離散傅立葉逆變換的函數(shù)的調(diào)用方法是：Y=ifft(X,n,dim)函數(shù)對(duì)X進(jìn)行離散傅立葉逆變換。其中X、n、dim的意義及用法和離散傅立葉變換函數(shù)fft完全相同。目錄目錄例5.39 對(duì)矩陣A的列向量、行向量分別進(jìn)行離散傅立葉變換、并對(duì)變換結(jié)果進(jìn)行逆變換。命令如下：

59、A=3,2,1,1;-5,1,0,1;3,2,1,5;fftA=fft(A) %求A的列向量的傅立葉變換fftA2=fft(A,4,2) %求A的行向量的傅立葉變換ifft(fftA) %對(duì)矩陣fftA的列向量進(jìn)行傅立葉逆變換，結(jié)果應(yīng)等于Aifft(fftA2,4,2) %對(duì)矩陣fftA2的行向量進(jìn)行傅立葉逆變換，其結(jié)果應(yīng)等于A目錄目錄5.6 數(shù)值微積分5.6.1 數(shù)值微分MATLAB中，沒(méi)有直接提供求數(shù)值導(dǎo)數(shù)的函數(shù)，只有計(jì)算向前差分的函數(shù)。DX=diff(X) 計(jì)算向量X的向前差分，DX(i)=X(i+1)-X(i)，0in。DX=diff(X,n) 計(jì)算X的n階向前差分，diff(X,2)

60、=diff(diff(X)。DX=diff(A,n,dim) 計(jì)算矩陣A的n階差分，dim=1時(shí)(缺省狀態(tài))，按列計(jì)算差分，dim=2，按行計(jì)算差分。目錄目錄例5.40 求向量sin(X)的13階差分。設(shè)X由0，2間均勻分布的10個(gè)點(diǎn)組成。命令如下：X=linspace(0,2*pi,10);Y=sin(X);DY=diff(Y); %計(jì)算Y的一階差分D2Y=diff(Y,2); %計(jì)算Y的二階差分，也可用命令diff(DY)計(jì)算D3Y=diff(Y,3); %計(jì)算Y的三階差分，也可用diff(D2Y)或diff(DY,2)目錄目錄例5.41 用不同的方法求函數(shù)f(x)的數(shù)值導(dǎo)數(shù)，并在同一個(gè)坐