1、MATLAB程序設計基礎MATLAB的數值計算Matlab的數據類型o 變量 變量不需要事先聲明，也不需要指定變量類型，它會自動根據所賦予變量的值或對變量的操作來確定變量的類型；賦值過程中，如果變量已存在，則用新值代替舊值，以新的類型代替舊的類型。 變量的命名規則：l 變量名區分大小寫；l 變量名長度不超過31位，第31位之后的字符被忽略；l 變量名以字母開頭，變量名中可以包含字母、數字、下劃線，但不能使用標點。 變量一般為局部變量，即僅在其調用的M文件內部有效；若要定義全局變量，須在變量前加關鍵字global。o 常量 matlab中預定義的一些特殊的量。 i,j 虛數單位Realmin 最

2、小的正浮點數， pi 圓周率Realmax 最大的浮點數， eps 浮點運算的相對精度Inf 無窮大 NaN not a number ，不定值 例如：pians = 3.14161102221023252101/0Warning: Divide by zero.ans = Inf0/0Warning: Divide by zero.ans = NaNo 定義變量時應避免與常量名相同，如果改變了某個常量的值，可以用clear命令來恢復。pi=1pi = 1clear pipians = 3.1416o 數字變量l 數字變量的運算258*369ans = 95202x=258*369x = 95

3、2021233ans = 1860867sqrt(ans)ans = 1.3641e+003l 數字的輸入輸出格式 缺省為實數保留小數點后4位浮點數表示。 其輸入格式與C語言一致： 如：9 -73 0.1999 1.475e6 輸出格式由format命令控制，只是影響屏幕顯示效果，不影響內部存儲和計算。format long;pians = 3.14159265358979format long e;pians = 3.141592653589793e+000format long g;pians = 3.14159265358979o 字符串l 1、字符串的約定字符串用單引號輸入或賦值；字符

4、串的每個字符都是都是字符數組的一個元素；字符串和字符數組基本上等價。s=symbolics =symbolic size(s)ans = 1 8 s(3)ans =ml 字符串的轉換 double 字符串轉換為數值代碼 num2str 數字轉換為字符串 int2str 整數轉換為字符串 mat2str 矩陣轉換為字符串 str2num 轉換字符串為數字 double(s)ans = 49 50 51 50 51 52l 字符串操作 strcat strcmp strvcat strncmp findstr upper lower blanks deblank l 執行字符串 t=1/(a*b-

5、1);a=2;b=3;c=eval(t)c = 0.2000o 結構型變量 由函數struct定義，以指針操作符“.”連接結構型變量名與屬性名。 結構型變量名結構型變量名struct（元素名（元素名1，元素值，元素值1，元素名，元素名2，元素值，元素值2，） c=struct(c1,1,c2,1 2 3 4,c3,abcd)c = c1: 1 c2: 1 2 3 4 c3: abcd c.c2ans = 1 2 3 4 c.c3ans =abcdabcdo 單元型變量 單元型變量為任意類型的多維數組，其定義需用大括號，元素間用逗號隔開。 a=1,2;3,4a = 1 2 3 4 b=1:4,a

6、,abcdb = 1x4 double 2x2 double abcd?cellplot(b) 單元型變量元素的引用采用大括號為下標標識，用小括號只顯示該元素的壓縮形式。 b2ans = 1 2 3 4 b(2)ans = 2x2 doubleo 向量 向量元素用“”括起來，元素間用空格、逗號或分號分隔； 注意：空格和逗號分隔成行向量，分號分割成列向量。注意：空格和逗號分隔成行向量，分號分割成列向量。l冒號表達式生成向量 基本格式：xx1:step:x2 xx1:x2 a=1:2:12a = 1 3 5 7 9 11 a=12:-2:1a = 12 10 8 6 4 2 a=1:6a = 1

7、2 3 4 5 6l線性等分向量生成 y=linspace(x1,x2) 生成100維行向量 y=linspace(x1,x2,n) 生成n維行向量 a=linspace(1,100,6)a = 1.0000 20.8000 40.6000 60.4000 80.2000 100.0000l對數等分向量生成 y=logspace(x1,x2)生成50維對數等分向量，y(1)=10 x1 y(50)=10 x2 y=logspace(x1,x2,n)生成n維對數等分向量y(1)=10 x1 y(n)=10 x2 a=logspace(0,5,6)a = 1 10 100 1000 10000 1

8、00000o 向量的基本運算l 與數運算a = 1.0000 20.8000 40.6000 60.4000 80.2000 100.0000 a-1ans = 0 19.8000 39.6000 59.4000 79.2000 99.0000 a*2ans = 2.0000 41.6000 81.2000 120.8000 160.4000 200.0000l 點積計算 指兩個向量在其中一個向量方向上的投影的乘積。 dot(a,b) a,b必須同維。 a=1 2 3; b=3,4,5; dot(a,b)ans = 26 sum(a.*b)ans = 26l 叉積 表示過兩相交向量的交點的垂直

9、于兩向量所在平面的向量。 cross（a，b） a，b必須為三維向量。l 混合積 c=cross(a,b)c = -2 4 -2 dot(a,cross(b,c)ans = 24o 矩陣大型矩陣通借助M文件來輸入。矩陣987654321A的輸入： A = 1,2,3; 4,5,6; 7,8,9A = 1 2 3 4 5 6 7 8 9 a=1 2 34 5 67 8 9a = 1 2 3 4 5 6 7 8 9x=rand(1,5)%產生的均布隨機數組產生的均布隨機數組 x =0.9501 0.2311 0.6068 0.4860 0.8913 x(3)%尋訪數組尋訪數組x的第三個元素。的第三

10、個元素。 ans =0.6068 x(1 2 5)%尋訪數組尋訪數組x的第一、二、五個元素組成的子數組。的第一、二、五個元素組成的子數組。 ans =0.9501 0.2311 0.8913 x(1:3)%尋訪前三個元素組成的子數組尋訪前三個元素組成的子數組 ans = 0.9501 0.2311 0.6068 x(3:end)%尋訪除前尋訪除前2個元素外的全部其他元素。個元素外的全部其他元素。end是最后一是最后一個元素的下標。個元素的下標。 ans = 0.6068 0.4860 0.8913 常用的特殊矩陣單位矩陣：eye(m,n); eye(m)零 矩 陣：zeros(m,n); ze

11、ros(m)一 矩 陣：ones(m,n); ones(m)對角矩陣：對角元素向量 V=a1,a2,an A=diag(V)隨機矩陣：rand(m,n)產生一個mn的均勻分別的隨機矩陣 eye(2,3)ans= 1 0 0 0 1 0 zeros(2,3)ans= 0 0 0 0 0 0 ones(2,3)ans= 1 1 1 1 1 1 V=5 7 2; A=diag(V)A= 5 0 0 0 7 0 0 0 2 eye(2)ans= 1 0 0 1 zeros(2)ans= 0 0 0 0 ones(2)ans= 1 1 1 1 如果已知A為方陣，則V=diag(A)可以提取A的對角元素構

12、成向量V。其他特殊矩陣 compan 友矩陣函數 magic 魔方矩陣 hankel Hankel矩陣 rosser 對稱特征值測試矩陣 hilb Hilbert矩陣 pascal Pascal矩陣 invhilb 反Hilbert矩陣 vander 范德蒙矩陣 o 矩陣的基本運算l 加減運算 要求兩矩陣必須同階。 a=1 2 3;2 3 4; 3 4 5; b=1 1 1;2 2 2;3 3 3; c=a+bc = 2 3 4 4 5 6 6 7 8l 乘法 要求a為ij階，b為jk階時，ab才能相乘。 e=b,5 5 5e = 1 1 1 5 2 2 2 5 3 3 3 5 f=a*ef

13、= 14 14 14 30 20 20 20 45 26 26 26 60l 除法左除“”: 相當于Ax=B的解，x=A-1B。右除“/”：相當于xA=B的解，x=BA-1 A-1B=(BA-1)。通常，右除稍快一些，而左除可以避免奇異性。對于AxB，其中A為（nm）階矩陣：n=m且非奇異時，方程為恰定方程；nm方程為超定方程；n A=1 2 3;4 5 6;7 8 0;1 3 5; B=1 3 5;2 4 6; A/Bans = 0 0.5000 -3.0000 3.5000 -12.0000 10.2500 1.0000 0.0000(BA)ans = 0 0.5000 -3.0000 3

14、.5000 -12.0000 10.2500 1.0000 0.0000 rand(seed,12); a=rand(100)+1.e8; x=ones(100,1); b=a*x; cond(a)ans = 5.0482e+011 tic;x1=b/a;t1=toct1 = 0.0773 er1=norm(x-x1)er1 = 139.8328 tic;x1=ab;t1=toct1 = 0.0014 er2=norm(x-x1)er2 = 1.1311e-004l 矩陣與常數的運算 常數與此矩陣的各元素之間進行運算。 注意：進行數除時，常數通常只能做除數。l 矩陣的逆運算 函數 inv A=

15、2 1 -3 -1;3 1 0 7;-1 2 4 -2;1 0 -1 5; inv(A)ans = -0.0471 0.5882 -0.2706 -0.9412 0.3882 -0.3529 0.4824 0.7647 -0.2235 0.2941 -0.0353 -0.4706 -0.0353 -0.0588 0.0471 0.2941l 矩陣的行列式運算 函數 det A=2 1 -3 -1;3 1 0 7;-1 2 4 -2;1 0 -1 5; a1=det(A)a1 = -85 a2=det(inv(A)a2 = -0.0118 a1*a2ans = 1l 矩陣的冪運算 與數字的冪運算

16、形式相同，用“”算符。l 矩陣的指數運算 常用函數 expm expm1 expm2 expm3l 矩陣的對數運算 函數 logml 矩陣的開方運算 函數 sqrtm b=magic(3)b = 8 1 6 3 5 7 4 9 2 sqrtm(b)ans = 2.7065 + 0.0601i 0.0185 + 0.5347i 1.1480 - 0.5948i 0.4703 + 0.0829i 2.0288 + 0.7378i 1.3739 - 0.8207i 0.6962 - 0.1430i 1.8257 - 1.2725i 1.3511 + 1.4155i b0.5ans = 2.7065

17、+ 0.0601i 0.0185 + 0.5347i 1.1480 - 0.5948i 0.4703 + 0.0829i 2.0288 + 0.7378i 1.3739 - 0.8207i 0.6962 - 0.1430i 1.8257 - 1.2725i 1.3511 + 1.4155io 矩陣的基本函數運算l 特征值函數 函數 x,y=eig(A) 可以給出特征值和特征向量的值 x為特征向量矩陣，y為特征值矩陣。 A=7 3 -2;3 4 -1;-2 -1 3;x,y=eig(A)x = 0.5774 0.0988 -0.8105 -0.5774 -0.6525 -0.4908 0.577

18、4 -0.7513 0.3197y = 2.0000 0 0 0 2.3944 0 0 0 9.6056l 奇異值函數 函數 svd svds 矩陣非零特征值的平方根稱作矩陣的奇異值。 矩陣的奇異值大小通常決定矩陣的性態，如果矩陣的奇異值變化特別大，則矩陣中某個參數的微小變化將嚴重影響到原矩陣的參數，這樣的矩陣又稱為病態矩陣，也稱奇異矩陣。例：A=1,1;5*eps,0;0,5*eps; sqrt(eig(A*A)l 矩陣翻轉 函數 fliplr flipud rot90a = 7 3 -2 3 4 -1 -2 -1 3 fliplr(a)ans = -2 3 7 -1 4 3 3 -1 -2

19、 flipud(a)ans = -2 -1 3 3 4 -1 7 3 -2 rot90(a)ans = -2 -1 3 3 4 -1 7 3 -2l 條件數函數 判斷矩陣的“病態”程度。 函數cond 計算矩陣的條件數的值 condest計算矩陣的1范數條件數的估計值 rcond 計算矩陣的條件數的倒數值 矩陣的條件數condA=最大奇異值/最小奇異值，條件數越大，則對參數變化越敏感。l范數函數 函數 norm（X，P） P1 1范數 P2 2范數 Pinf 無窮范數 Pfro F范數 norm（X）norm（X，2）njijniTniijnjaAAAsAaA11max2111|max|max

20、l 秩函數 函數 ranke = 1 1 1 5 2 2 2 2 3 3 3 5 rank(e)ans = 2l 跡函數 矩陣所有對角線上元素的和稱為矩陣的跡。 函數 tracel 正交空間函數 函數 orth 用來求矩陣的一組正交基。l 偽逆函數 函數 pinv 求解“病態”問題時，避免產生偽解。 a=magic(4)a = 16 2 3 13 5 11 10 8 9 7 6 12 4 14 15 1 b=a*1 1 1 1; inv(a)*bWarning: Matrix is close to singular or badly scaled. Results may be inaccu

21、rate. RCOND = 1.567374e-017.ans = 0 8 0 0 pinv(a)*bans = 1.0000 1.0000 1.0000 1.0000o 稀疏矩陣 用元素的行列號和元素值來存儲一個元素。 a=speye(100)a = (1,1) 1 (2,2) 1 (3,3) 1 (4,4) 1 (99,99) 1 (100,100) 1 b=eye(100); whos a b Name Size Bytes Class a 100 x100 1604 sparse array b 100 x100 80000 double arrayGrand total is 101

22、00 elements using 81604 bytesl 常規矩陣轉換為稀疏矩陣 sparse(A) a=1 2 3;4 5 6;7 8 0; b=sparse(a)b = (1,1) 1 (2,1) 4 (3,1) 7 (1,2) 2 (2,2) 5 (3,2) 8 (1,3) 3 (2,3) 6 whos a b Name Size Bytes Class a 3x3 72 double array b 3x3 112 sparse arrayGrand total is 17 elements using 184 bytes sparse(i,j,s,n,m,nzmax) spars

23、e(i,j,s,n,m) sparse(i,j,s,n) n,m為生成稀疏矩陣的行列數，i,j,s為子矩陣，nzmax最多的非零元素數l 稀疏矩陣轉換為常規矩陣 full（）l 通用函數形式 通用函數調用格式 funm（A，funname） funname包括sin sinh asin asinh cos cosh acos acosh tan exp log log2 pow2 sqrt abs funm(a,sqrt)ans = 3.7584 - 0.2071i -0.2271 + 0.4886i 0.3887 + 0.7700i 1.9110 - 1.0514i 0.2745 - 0.0

24、130i 2.3243 + 0.0306i 2.0076 + 0.0483i 1.2246 - 0.0659i 1.3918 - 0.2331i 1.5060 + 0.5498i 1.4884 + 0.8666i 1.4447 - 1.1833i 0.4063 + 0.4533i 2.2277 - 1.0691i 1.9463 - 1.6848i 1.2506 + 2.3006i sqrtm(a)ans = 3.7584 - 0.2071i -0.2271 + 0.4886i 0.3887 + 0.7700i 1.9110 - 1.0514i 0.2745 - 0.0130i 2.3243 +

25、 0.0306i 2.0076 + 0.0483i 1.2246 - 0.0659i 1.3918 - 0.2331i 1.5060 + 0.5498i 1.4884 + 0.8666i 1.4447 - 1.1833i 0.4063 + 0.4533i 2.2277 - 1.0691i 1.9463 - 1.6848i 1.2506 + 2.3006i sqrt(a)ans = 4.0000 1.4142 1.7321 3.6056 2.2361 3.3166 3.1623 2.8284 3.0000 2.6458 2.4495 3.4641 2.0000 3.7417 3.8730 1.0

26、000p 矩陣分解函數l 特征值分解V,D=eig(X) 矩陣的特征值分解：XVVDV,D=eig(X,nobalance) 關閉平衡算法的求解方法（平衡算法對于某些問題可以得到更高的精度）。V,D=eig(A,B) 廣義特征值分解：AVBVD a=-149 -50 -154;537 180 546;-27 -9 -25;v,d=eig(a)v = 0.3162 0.4041 0.1391 -0.9487 -0.9091 -0.9740 0.0000 -0.1010 0.1789d = 1.0000 0 0 0 2.0000 0 0 0 3.0000 b=2 10 2;10 5 -8;2 -8

27、 11;v,d=eig(a,b)v = 0.8211 -0.3138 -0.0191 -0.3452 0.9495 -0.9441 -0.4546 -0.0044 0.3290d = 12.9030 0 0 0 -0.0045 0 0 0 0.0706l 奇異值分解 U,S,V=svd(X) 其中XUSV a=1;1;U,S,V=svd(a)U = 0.7071 -0.7071 0.7071 0.7071S = 1.4142 0V = 1 l LU分解 L,U=lu(A) 又稱三角分解，目的是分解成一個下三角陣L和一個上三角陣U的乘積，即ALU a=1 2 3;2 4 1;4 6 7;l,u=

28、lu(a)l = 0.2500 0.5000 1.0000 0.5000 1.0000 0 1.0000 0 0u = 4.0000 6.0000 7.0000 0 1.0000 -2.5000 0 0 2.5000注意：L實際上是一個“心理上”的下三角矩陣，它事實上是一個置換矩陣P的逆矩陣與一個真正下三角矩陣L1（其對角線元素為1）的乘積。l,u,p=lu(a)l = 1.0000 0 0 0.5000 1.0000 0 0.2500 0.5000 1.0000u = 4.0000 6.0000 7.0000 0 1.0000 -2.5000 0 0 2.5000p = 0 0 1 0 1

29、0 1 0 0 inv(p)*l*uans = 1 2 3 2 4 1 4 6 7p為置換矩陣，此時滿足AP-1LUl Chol分解 如果A為n階對稱正定矩陣，則存在一個非奇異下三角實矩陣L，使得ALLT,當限定L的對角元素為正時，這種分解是唯一的。 a=4 -1 1;-1 4.25 2.75;1 2.75 3.5; chol(a)ans = 2.0000 -0.5000 0.5000 0 2.0000 1.5000 0 0 1.0000l 正交分解 AQR a=1 1 1;2 -1 -1;2 -4 5;q,r=qr(a)q = -0.3333 -0.6667 -0.6667 -0.6667

30、-0.3333 0.6667 -0.6667 0.6667 -0.3333r = -3 3 -3 0 -3 3 0 0 -3將矩陣A做正交化分 解 ， 使 得Q*R=A，其中Q為正交矩陣（其范 數 為 1 ， 指 令norm(Q)=1)，R為對角化的上三角矩陣。矩陣的特殊操作p 變維lreshape(X,M,N) X變為MN維lreshape(X,M,N,P,)X變為MNP 或reshape(X,M N P )l“：” 操作符 a=1:12; b=reshape(a,2,6)b = 1 3 5 7 9 11 2 4 6 8 10 12 c=zeros(4,3); c(:)=a(:)c = 1

31、5 9 2 6 10 3 7 11 4 8 12p 矩陣抽取l 對角元素抽取 diag（X,k）抽取X的第k條對角線元素，k=0為主對角線，上對角線為正值，下對角線為負值。 diag（X） 抽取主對角線元素 a=pascal(4)a = 1 1 1 1 1 2 3 4 1 3 6 10 1 4 10 20 v=diag(a)v = 1 2 6 20 v=diag(a,2)v = 1 4l 三角陣的抽取 tril(X)提取X的主下三角部分 tril(X,k) 提取X的第k條對角線下面的元素 triu(X)提取X的主上三角部分 triu(X,k) 提取X的第k條對角線上面的元素 al=tril(a

32、,-1)al = 0 0 0 0 1 0 0 0 1 3 0 0 1 4 10 0 au=triu(a,-1)au = 1 1 1 1 1 2 3 4 0 3 6 10 0 0 10 20數組及其運算 數組的建立、存儲完全同于矩陣，由于計算的不同，把相同型矩陣之間的運算稱為數組運算。p 基本數組運算l 數組的四則運算 普通運算同矩陣的運算，另有點運算“.*”、”./”、”.”，即兩數組對應元素之間的運算。 a1=1 2 3;2 3 4;3 4 5; b1=1 1 1;2 2 2;3 3 3; a1./b1ans = 1.0000 2.0000 3.0000 1.0000 1.5000 2.0000 1.0000 1.3333 1.6667 a1.b1ans = 1.0000 0.5000 0.3333 1.0000 0.6667 0.5000 1.0000 0.7500 0.6000l 數組