第4講 MATLAB數(shù)值計(jì)算

第4講MATLAB的數(shù)值計(jì)算4.1特殊矩陣

4.1.1對(duì)角陣與三角陣(diag、tril、triu)

4.1.2特殊矩陣(*)4.2矩陣分析

4.2.1矩陣結(jié)構(gòu)變換('、rot90、fliplr、flipud)

4.2.2矩陣的逆與偽逆(inv、pinv)

4.2.3方陣的行列式(det)

4.2.4矩陣的秩(rank)

4.2.5矩陣的范數(shù)(norm)

4.2.6矩陣的特征值與特征向量(eig)4.3矩陣分解與線(xiàn)性方程組求解4.3.1矩陣分解(chol、lu)4.3.2線(xiàn)性方程組求解4.4多項(xiàng)式計(jì)算多項(xiàng)式的表示、根本運(yùn)算〔加、減、乘、除、求根、求值與數(shù)組運(yùn)算、求導(dǎo)〕、特征多項(xiàng)式與矩陣運(yùn)算4.5數(shù)據(jù)處理數(shù)值插值、曲線(xiàn)擬合4.6數(shù)據(jù)統(tǒng)計(jì)與分析

最大和最小值、平均值和中值、求和與求積、累加和與累乘積、標(biāo)準(zhǔn)方差、元素排序4.7泛函指令

4.4.1函數(shù)在Matlab中的表示

4.4.2函數(shù)的極值點(diǎn)

4.4.3數(shù)值微積分

4.4.4常微分方程的數(shù)值求解4.1特殊矩陣4.1.1對(duì)角陣與三角陣1.矩陣的對(duì)角元素(1)提取矩陣的對(duì)角線(xiàn)元素設(shè)A為m×n矩陣，diag(A)函數(shù)用于提取矩陣A主對(duì)角線(xiàn)元素產(chǎn)生一個(gè)具有min(m,n)個(gè)元素的列向量。diag(A)函數(shù)還有更進(jìn)一步的形式diag(A,k)，其功能是提取第k條對(duì)角線(xiàn)的元素。(2)構(gòu)造對(duì)角矩陣設(shè)V為具有m個(gè)元素的向量，diag(V)將產(chǎn)生一個(gè)m×m對(duì)角矩陣，其主對(duì)角線(xiàn)元素即為向量V的元素。diag(V)函數(shù)也有更進(jìn)一步的形式diag(V,k)，其功能是產(chǎn)生一個(gè)n×n(n=m+k)對(duì)角陣，其第k條對(duì)角線(xiàn)的元素即為向量V的元素。k=0k<0k>0m<nm>n

例4.1先建立5×5矩陣A，然后將A的第1行元素乘以1，第2行乘以2，…，第5行乘以5。命令如下：A=[17,0,1,0,1523,5,7,14,164,0,13,0,2210,12,19,21,311,18,25,2,19];D=diag([1,2,3,4,5]);D*A問(wèn)：如何將A的第1列元素乘以1，第2列乘以2，…，第5列乘以5。

2.矩陣的三角陣(1)下三角矩陣〔tril〕求矩陣A的下三角陣的MATLAB函數(shù)是tril(A)。tril(A)函數(shù)也有更進(jìn)一步的一種形式tril(A,k)，其功能是求矩陣A的第k條對(duì)角線(xiàn)以下的元素。tril(A)tril(A,1);tril(A,-1)(2)上三角矩陣〔triu〕在MATLAB中，提取矩陣A的上三角矩陣的函數(shù)是triu(A)和triu(A,k)，其用法與提取下三角矩陣的函數(shù)tril(A)和tril(A,k)完全相同。

4.1.2特殊矩陣的生成(*)1.魔方矩陣--magic(n)其功能是生成一個(gè)n階魔方陣。例4.2生成一個(gè)5行5列的魔方陣〔其每行每列及對(duì)角線(xiàn)的和均為65〕。命令如下：B=magic(5)2.范得蒙矩陣--vander(V)生成以向量V為根底向量的范得蒙矩陣,其中：A(i,j)=v(i)^(n-j),wheren=length(v)

3.希爾伯特矩陣--hilb(n)。其中：H(i,j)=1/(i+j-1).MATLAB中，有一個(gè)專(zhuān)門(mén)求希爾伯特矩陣的逆的函數(shù)invhilb(n)，其功能是求n階的希爾伯特矩陣的逆矩陣。4.帕斯卡矩陣--pascal(n)

例4.3求(x+y)5的展開(kāi)式。在MATLAB命令窗口，輸入命令：pascal(6)ans=111111123456136101521141020355615153570126162156126252其次對(duì)角線(xiàn)上的元素1，5，10，10，5，1即為展開(kāi)式的系數(shù)。4.2矩陣分析4.2.1矩陣結(jié)構(gòu)變換1.矩陣的轉(zhuǎn)置轉(zhuǎn)置運(yùn)算符是單撇號(hào)(')。2.矩陣的旋轉(zhuǎn)矩陣的旋轉(zhuǎn)利用函數(shù)rot90(A,k)，功能是將矩陣A逆時(shí)針旋轉(zhuǎn)90o的k倍，當(dāng)k為1時(shí)可省略。3.矩陣的左右翻轉(zhuǎn)對(duì)矩陣A實(shí)施左右翻轉(zhuǎn)的函數(shù)是fliplr(A)。4.矩陣的上下翻轉(zhuǎn)對(duì)矩陣A實(shí)施上下翻轉(zhuǎn)的函數(shù)是flipud(A)。

4.2.2矩陣的逆與偽逆1.矩陣的逆求一個(gè)矩陣的逆非常容易。求方陣A的逆可調(diào)用函數(shù)inv(A)。A=[123;221;343];Y=inv(A)Y=A^(-1)例4.4用求逆矩陣的方法解線(xiàn)性方程組。命令如下：A=[1,2,3;1,4,9;1,8,27];b=[5,-2,6]';x=inv(A)*b一般情況下，用左除比求矩陣的逆的方法更有效，即x=A\b。

2.矩陣的偽逆〔用于非方陣或奇異方陣求偽逆〕MATLAB中，求一個(gè)矩陣偽逆的函數(shù)是pinv(A)。例4.5求A的偽逆，并將結(jié)果送B。命令如下：A=[3,1,1,1;1,3,1,1;1,1,3,1];B=pinv(A)例4.6求矩陣A的偽逆。在MATLAB命令窗口，輸入命令：A=[0,0,0;0,1,0;0,0,1];pinv(A)

4.2.3方陣的行列式求方陣A所對(duì)應(yīng)的行列式的值的函數(shù)是det(A)。例4.7用克萊姆(Cramer)方法求解線(xiàn)性方程組。程序如下：D=[2,2,-1,1;4,3,-1,2;8,5,-3,4;3,3,-2,2];%定義系數(shù)矩陣b=[4;6;12;6];%定義常數(shù)項(xiàng)向量D1=[b,D(:,2:4)];%用方程組的右端向量置換D的第1列D2=[D(:,1:1),b,D(:,3:4)];D3=[D(:,1:2),b,D(:,4:4)];D4=[D(:,1:3),b];DD=det(D);%前提是系數(shù)矩陣的行列式不為0x1=det(D1)/DD;x2=det(D2)/DD;x3=det(D3)/DD;x4=det(D4)/DD;[x1,x2,x3,x4]

4.2.4矩陣的秩MATLAB中，求矩陣秩的函數(shù)是rank(A)。例如，求例4.7中方程組系數(shù)矩陣D的秩，命令是：D=[2,2,-1,14,3,-1,28,5,-3,43,3,-2,2];r=rank(D)說(shuō)明D是一個(gè)滿(mǎn)秩矩陣。注：rank(A)的算法是以奇異值分解為根底的。

4.2.5矩陣的范數(shù)

矩陣范數(shù)是對(duì)矩陣元素的數(shù)量級(jí)大小的一種度量，通過(guò)這些度量，可以對(duì)線(xiàn)性方程組解的誤差進(jìn)行估計(jì)。1.計(jì)算矩陣3種常用范數(shù)的函數(shù)(1)norm(A)或norm(A,2)2—范數(shù)，即求矩陣的最大奇異值，=max(svd(A))。

(2)norm(A,1)1—范數(shù),列范數(shù)，即求列元素和的最大值，=max(sum(abs(A)))(3)norm(A,inf)∞—范數(shù)，行范數(shù)，即行元素和的最大值，=max(sum(abs(A')))

例4.8求矩陣A的三種范數(shù)。命令如下：A=[17,0,1,0,1523,5,7,14,164,0,13,0,2210,12,19,21,311,18,25,2,19];a1=norm(A,1)%求A的1—范數(shù)a2=norm(A)%求A的2—范數(shù)ainf=norm(A,inf)%求A的∞—范數(shù)4.2.6矩陣〔方陣〕的特征值與特征向量特征值和特征向量是線(xiàn)性代數(shù)中非常重要的概念，在求解數(shù)學(xué)問(wèn)題和在實(shí)際的工程應(yīng)用中占有非常重要的地位。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的列向量，使AV=VD成立，如果Ｖ為非奇異，那么有如下特征值分解：Ａ＝ＶＤＶ-1(3)[V,D]=eig(A,'nobalance')%當(dāng)矩陣A中有與截?cái)嗾`差數(shù)量級(jí)相差不遠(yuǎn)的值時(shí)，該指令可能更精確。'nobalance'起誤差調(diào)節(jié)作用。例4.9用3種不同的格式求A的特征值和特征向量。命令如下：A=[1,2,2;1,-1,1;4,-12,1];E=eig(A)[V,D]=eig(A)[V,D]=eig(A,'nobalance')4.3矩陣分解與線(xiàn)性方程組求解矩陣分解1、Cholesky分解(適用于對(duì)稱(chēng)正定矩陣的分解)格式R=chol(X)%如果X為n階對(duì)稱(chēng)正定矩陣，那么存在一個(gè)實(shí)的非奇異上三角陣R，滿(mǎn)足R'*R=X；假設(shè)X非正定，那么產(chǎn)生錯(cuò)誤信息。[R,p]=chol(X)%不產(chǎn)生任何錯(cuò)誤信息，假設(shè)X為正定陣，那么p=0，R與上相同；假設(shè)X非正定，那么p為正整數(shù)，R是一個(gè)q=p-1階的上三角陣，使得R'*R=X(1:q,1:q)。例4.10

X=pascal(4)%產(chǎn)生4階pascal矩陣[R,p]=chol(X)a=rand(4)b=a'+a%產(chǎn)生4階實(shí)對(duì)稱(chēng)矩陣，但不一定正定[r,p]=chol(b)r'*r

2.矩陣的LU分解矩陣的三角分解又稱(chēng)LU分解，它的目的是將一個(gè)矩陣分解成一個(gè)下三角矩陣L和一個(gè)上三角矩陣U的乘積，即A=LU。格式[L,U]=lu(X)%U為上三角陣，L為下三角陣或其變換形式，滿(mǎn)足LU=X。[L,U,P]=lu(X)%U為上三角陣，L為下三角陣，P為單位矩陣的行變換矩陣，滿(mǎn)足LU=PX。例4.11

A=[123;456;789];[L,U]=lu(A)[L,U,P]=lu(A)3、利用矩陣的LU和cholesky分解求線(xiàn)性方程組的解〔1〕LU分解：LU分解又稱(chēng)Gauss消去分解，可把任意方陣分解為下三角矩陣和上三角矩陣的乘積。即A=LU，L為下三角陣，U為上三角陣。那么：A*X=b變成L*U*X=b，所以X=U\(L\b)這樣可以大大提高運(yùn)算速度。(實(shí)際上，除法求線(xiàn)性方程組的解：X=A\b就是用LU分解實(shí)現(xiàn)的。)例4.12求方程組的解。A=[42-1;3-12;1130];B=[2108]';D=det(A)[L,U]=lu(A)X=U\(L\B)說(shuō)明結(jié)果中的警告是由于系數(shù)行列式為零產(chǎn)生的。可以通過(guò)A*X驗(yàn)證其正確性。

〔2〕Cholesky分解假設(shè)A為對(duì)稱(chēng)正定矩陣，那么Cholesky分解可將矩陣A分解成上三角矩陣和其轉(zhuǎn)置的乘積，即：其中R為上三角陣。方程A*X=b變成R’*R*X=b所以X=R\(R’\b)

4.3.2線(xiàn)性方程組求解線(xiàn)性方程組解的一般討論

將線(xiàn)性方程的求解分為兩類(lèi)：一類(lèi)是方程組求唯一解或求特解，另一類(lèi)是方程組求無(wú)窮解即通解?？梢酝ㄟ^(guò)系數(shù)矩陣的秩〔用rank函數(shù)求矩陣的秩〕來(lái)判斷：假設(shè)系數(shù)矩陣的秩r=n〔恰定方程，n為方程組中未知變量的個(gè)數(shù)〕，那么有唯一解；假設(shè)系數(shù)矩陣的秩r<n〔欠定方程〕，那么可能有無(wú)窮解；而對(duì)于超定方程，解常常是矛盾的，一般只能得到其最小二乘解；線(xiàn)性方程組的無(wú)窮解=對(duì)應(yīng)齊次方程組的通解+非齊次方程組的一個(gè)特解；其特解的求法屬于解的第一類(lèi)問(wèn)題，通解局部屬第二類(lèi)問(wèn)題。２、求線(xiàn)性方程組的唯一解或特解〔第一類(lèi)問(wèn)題〕〔１〕利用矩陣除法求線(xiàn)性方程組的特解〔或一個(gè)解〕方程：AX=b解法：X=A\b例4.13求方程組的解：解：A=[5600015600015600015600015];B=[10001]';R_A=rank(A)%求秩X=A\B%求解

3、解線(xiàn)性方程組的一般函數(shù)文件如下：function[x,y]=line_solution(A,b)[m,n]=size(A);y=[];ifnorm(b)>0%非齊次方程組ifrank(A)==rank([A,b])%方程組相容ifrank(A)==n%有唯一解(恰定〕x=A\b;else%方程組有無(wú)窮多個(gè)解,根底解系(欠定)disp('原方程組有有無(wú)窮個(gè)解，其齊次方程組的根底解系為y，特解為x');y=null(A,'r');x=A\b;endelse%方程組不相容，給出最小二乘法解(超定)disp('方程組的最小二乘法解是：');x=A\b;endelse%齊次方程組ifrank(A)>=n%列滿(mǎn)秩x=zero(m,1)%0解else%非0解disp('方程組有無(wú)窮個(gè)解，根底解系為x');x=null(A,'r');endendreturn2.應(yīng)用舉例例4.14求線(xiàn)性方程組的解。在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ù)例4.15求以下線(xiàn)性方程組的解。在MATLAB命令窗口，輸入命令：A=[2,7,3,1;3,5,2,2;9,4,1,7];b=[6,4,2]';[x,y]=line_solution(A,b)4.4多項(xiàng)式計(jì)算4.4.1(一元)多項(xiàng)式的表示在Matlab中，任意多項(xiàng)式都可以用它按降冪次序排列后對(duì)應(yīng)的系數(shù)行向量表示。例4.16用系數(shù)行向量表示以下多項(xiàng)式3x5-7x4+5x2+2x-18命令如下：p=[3,-7,0,5,2,-18];

4.4.2多項(xiàng)式的根本運(yùn)算1.多項(xiàng)式的四那么運(yùn)算(1)多項(xiàng)式的加減法(Matlab未提供現(xiàn)成的函數(shù))如：P1=5X3,P2=X5-2X4+5X,求P1+P2(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ù)向量。deconv是conv的逆函數(shù)，即有P1=conv(P2,Q)+r。例4.17設(shè)有兩個(gè)多項(xiàng)式，計(jì)算：f(x)=3x5-5x4+2x3-7x2+5x+6g(x)=3x2+5x-3(1)求f(x)+g(x)、f(x)-g(x)。(2)求f(x)·g(x)、f(x)/g(x)。在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，余式送rconv(Q,g)+r%驗(yàn)證注：參照本例，可編寫(xiě)一函數(shù)，實(shí)現(xiàn)多項(xiàng)式的加減運(yùn)算。(問(wèn)題1)

例4.18展開(kāi)多項(xiàng)式w=conv([1,2,2],conv([1,4],[1,1]))P=poly2str(w,'s')%將w表示成字符串形式的多項(xiàng)式，s為變量名問(wèn)題2：求的商及余多項(xiàng)式

2.多項(xiàng)式求根求多項(xiàng)式p(x)的根的函數(shù)是roots(P)，這里，P是p(x)的系數(shù)向量，該函數(shù)返回方程p(x)=0的全部根(含重根，復(fù)根)。一個(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的向量。3.多項(xiàng)式求值與數(shù)組運(yùn)算求多項(xiàng)式p(x)在某點(diǎn)或某些點(diǎn)的函數(shù)值的函數(shù)是polyval(P,x)。假設(shè)x為一數(shù)值，那么求多項(xiàng)式在該點(diǎn)的值；假設(shè)x為向量或矩陣，那么對(duì)向量或矩陣中的每個(gè)元素求其多項(xiàng)式的值(數(shù)組運(yùn)算)。例4.19一個(gè)多項(xiàng)式，計(jì)算：f(x)=3x5+4x3-5x2–4.2x+5(1)計(jì)算f(x)=0的全部根。(2)由方程f(x)=0的根構(gòu)造一個(gè)多項(xiàng)式g(x)，并與f(x)進(jìn)行比照。(3)計(jì)算f(5)、f(4.8)、f(9.6)、f(12.3)的值(數(shù)組運(yùn)算)。命令如下：P=[3,0,4,-5,-4.2,5];X=roots(P)%求方程f(x)=0的根G=poly(X)%求多項(xiàng)式g(x),并與原f(x)進(jìn)行比較X0=[5,4.8,9.6,12.3];f=polyval(P,X0)%求多項(xiàng)式f(x)在給定點(diǎn)的值4.多項(xiàng)式的求導(dǎo)(略)對(duì)多項(xiàng)式求導(dǎo)數(shù)的函數(shù)是：p=polyder(P)求多項(xiàng)式P的導(dǎo)函數(shù)p=polyder(P,Q)求P*Q的導(dǎo)函數(shù),即：p=polyder(conv(P,Q))[p,q]=polyder(P,Q)求P/Q的導(dǎo)函數(shù)，導(dǎo)函數(shù)的分子存入p，分母存入q。即：q=conv(Q,Q)p=conv(polyder(P),Q)-conv(polyder(Q),P)例4.20求如下有理分式的導(dǎo)數(shù):

3x5+5x4-8x2+1x-510x5+5x4+6x2+7x-100命令如下：P=[3,5,0,-8,1,-5];Q=[10,5,0,6,7,-100];[p,q]=polyder(P,Q)conv(Q,Q)%驗(yàn)證qconv(polyder(P),Q)-conv(polyder(Q),P)%驗(yàn)證p4.4.3特征多項(xiàng)式與多項(xiàng)式矩陣運(yùn)算給定一矩陣A，求其特征多項(xiàng)式：P=poly(A)%P為A的特征多項(xiàng)式的系數(shù)行向量例4.21求矩陣A的特征多項(xiàng)式A=[111213;141516;171819];PA=poly(A)%PA為A的特征多項(xiàng)式roots(PA)%PA的根為A的特征值eig(A)%驗(yàn)證函數(shù)插值與逼近的任務(wù)是用簡(jiǎn)單函數(shù)〔主要是多項(xiàng)式和分段多項(xiàng)式〕為離散數(shù)據(jù)建立連續(xù)模型，為各種非有理函數(shù)提供好的逼近。4.5.1數(shù)值插值一維數(shù)值插值對(duì)一維數(shù)據(jù)通過(guò)插值求得離散點(diǎn)之間的數(shù)據(jù)量。這里的一維數(shù)據(jù)包括自變量的一維數(shù)據(jù)和因變量的一維數(shù)據(jù)：4.5數(shù)據(jù)處理y0y1y2yn……自變量的一維數(shù)據(jù)因變量的一維數(shù)據(jù)插值多項(xiàng)式插值曲線(xiàn)yixi一維插值函數(shù)調(diào)用格式為：Yi=interp1(X,Y,Xi,'method')函數(shù)根據(jù)X、Y的值，計(jì)算函數(shù)在Xi處的值。X、Y是兩個(gè)等長(zhǎng)的向量，分別描述采樣點(diǎn)和樣本值，其中X中的數(shù)據(jù)必須是單調(diào)的；Xi是一個(gè)向量或標(biāo)量，描述欲插值的點(diǎn)，Yi是一個(gè)與Xi等長(zhǎng)的插值結(jié)果。method是插值方法，允許的取值有‘linear’(線(xiàn)性插值)、‘nearest’(最近插值)、‘spline’(三次樣條插值)、‘cubic’〔三次多項(xiàng)式插值〕，缺省值是‘linear’。如果X中的數(shù)據(jù)是等間隔的，在每種method之前加*〔如‘*linear’〕，可使執(zhí)行速度加快。例4.22用不同的插值方法計(jì)算sin(x)在π/2點(diǎn)的值。這是一個(gè)一維插值問(wèn)題。在MATLAB命令窗口，輸入命令：X=0:0.2:pi;Y=sin(X);%給出X、Yinterp1(X,Y,pi/2)%用缺省方法(即線(xiàn)性插值方法)計(jì)算sin(π/2)interp1(X,Y,pi/2,'nearest')%用最近方法計(jì)算sin(π/2)interp1(X,Y,pi/2,'linear')%用線(xiàn)性方法計(jì)算sin(π/2)interp1(X,Y,pi/2,'spline')%用三次樣條方法計(jì)算sin(π/2)interp1(X,Y,pi/2,'cubic')%用三次多項(xiàng)式方法計(jì)算sin(π/2)注：MATLAB中有一個(gè)專(zhuān)門(mén)的三次樣條插值函數(shù)Yi=spline(X,Y,Xi)，其功能及使用方法與函數(shù)Yi=interp1(X,Y,Xi,'spline')完全相同。

例4.23檢測(cè)參數(shù)F隨時(shí)間T的采樣結(jié)果，用數(shù)值插值法計(jì)算X=2，7，12，17，22，17，32，37，42，47，52，57時(shí)F的值。這是一個(gè)一維數(shù)值插值問(wèn)題，命令如下：T=0:5:65;F=[3.2021,2.2560,879.5,1835.9,2968.8,4136.2,5234.9,6152.7,...6725.3,6848.3,6403.5,6824.7,7328.5,7854.6];X=2:5:57;F1=interp1(T,F,X);%用線(xiàn)性方法插值X0=0:0.1:65;F0=interp1(T,F,X0);%顯示采樣點(diǎn)、插值點(diǎn)并畫(huà)出插值曲線(xiàn)plot(T,F,'.r','MarkerSize',20);holdonplot(X0,F0,X,F1,'.b','MarkerSize',20)holdoff%用最近方法插值X=2:5:57;F1=interp1(T,F,X,'nearest');X0=0:0.1:65;F0=interp1(T,F,X0,'nearest');%顯示采樣點(diǎn)、插值點(diǎn)并畫(huà)出插值曲線(xiàn)plot(T,F,'.r','MarkerSize',20);holdonplot(X0,F0,X,F1,'.b','MarkerSize',20)holdoff%用三次樣條方法插值X=2:5:57;F1=interp1(T,F,X,'spline');X0=0:0.1:65;F0=interp1(T,F,X0,'spline');%顯示采樣點(diǎn)、插值點(diǎn)并畫(huà)出插值曲線(xiàn)plot(T,F,'.r','MarkerSize',20);holdonplot(X0,F0,X,F1,'.b','MarkerSize',20)holdoff%用三次多項(xiàng)式方法插值X=2:5:57;F1=interp1(T,F,X,'cubic');X0=0:0.1:65;F0=interp1(T,F,X0,'cubic');%顯示采樣點(diǎn)、插值點(diǎn)并畫(huà)出插值曲線(xiàn)plot(T,F,'.r','MarkerSize',20);holdonplot(X0,F0,X,F1,'.b','MarkerSize',20)holdoff問(wèn)題3:用隨機(jī)函數(shù)生成n個(gè)數(shù)據(jù)點(diǎn)，對(duì)這n個(gè)數(shù)據(jù)點(diǎn)做3次樣條插值并繪出插值曲線(xiàn)?！踩鏽=10〕2.二維數(shù)值插值MATLAB中，提供了解決二維插值問(wèn)題的函數(shù)。其調(diào)用格式為：Z1=interp2(X,Y,Z,X1,Y1,'method')其中X、Y、Z都是二維數(shù)組，Ｘ、Ｙ分別描述采樣點(diǎn)的兩個(gè)自變量的樣本(坐標(biāo))值，Z那么是采樣點(diǎn)的因變量的樣本(坐標(biāo))值，X1、Y1用于描述欲插值的點(diǎn)的自變量的坐標(biāo)值，分別可以是標(biāo)量、向量或二維數(shù)組。method的取值與一維插值函數(shù)相同。一維插值是在平面直角坐標(biāo)系〔二維平面〕中進(jìn)行的，插值的結(jié)果都在插值曲線(xiàn)上；而二維插值那么是在三維空間中進(jìn)行的，插值的結(jié)果都在插值曲面上。例4.24利用二維插值函數(shù)對(duì)peaks函數(shù)進(jìn)行插值。命令如下：[X,Y]=meshgrid(-3:.5:3);%產(chǎn)生分別描述采樣點(diǎn)的自變量Ｘ、Ｙ的坐標(biāo)值的兩個(gè)二維數(shù)組，其中X=Y'Z=peaks(X,Y);[X1,Y1]=meshgrid(-3:.125:3);Z1=interp2(X,Y,Z,X1,Y1,'spline');%對(duì)函數(shù)進(jìn)行插值mesh(X,Y,Z);holdon;mesh(X1,Y1,Z1+30);%下面局部為原函數(shù)圖形，上面局部為插值曲面holdoff4.5.2曲線(xiàn)擬合MATLAB中，提供了使用最小二乘法進(jìn)行曲線(xiàn)擬合的函數(shù)。調(diào)用格式為：[P,S]=polyfit(X,Y,m)函數(shù)根據(jù)采樣點(diǎn)的自變量向量X和采樣點(diǎn)的因變量向量Y，產(chǎn)生一個(gè)m次多項(xiàng)式P及其在采樣點(diǎn)的誤差向量S。其中X、Y是兩個(gè)等長(zhǎng)的向量，P是一個(gè)長(zhǎng)度為m+1的向量〔m次多項(xiàng)式的系數(shù)向量〕。例4.25用一個(gè)5次多項(xiàng)式在區(qū)間[0，2π]內(nèi)逼近函數(shù)sin(x)。命令如下：X=linspace(0,2*pi,50);Y=sin(X);[P,S]=polyfit(X,Y,5)%得到5次多項(xiàng)式的系數(shù)和誤差plot(X,Y,‘k*’,X,polyval(P,X),‘k-’)%用*號(hào)標(biāo)出樣本點(diǎn)holdonplot(X,polyval(P,X),'.r','MarkerSize',10)holdoff4.6數(shù)據(jù)統(tǒng)計(jì)與分析1.求矩陣最大和最小元素

max、min2.求矩陣的平均值和中值

mean、median3.矩陣元素求和與求積

sum和prod4.矩陣元素累加和與累乘積MATLAB中，使用cumsum和cumprod函數(shù)能方便地求得向量和矩陣元素的累加和與累乘積向量，函數(shù)的用法和sum及prod相同，區(qū)別僅在于累加、累乘的中間結(jié)果均保存。例4.26求向量X=(1，2，4，7，15)的累加和與累乘積向量。命令如下：X=[124715];cumsum(X)cumprod(X)5.標(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í)，那么求各行元素的標(biāo)準(zhǔn)方差。FLAG取0或1，分別對(duì)應(yīng)于標(biāo)準(zhǔn)方差的兩種定義：〔help〕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)行排序，假設(shè)dim=1，那么按列排，假設(shè)dim=2，那么按行排。Y是排序后的矩陣，而I記錄Y中的元素在A中位置。例4.27對(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ì)A按列排序，并將每個(gè)元素所在行號(hào)送矩陣I4.7泛函指令Matlab提供了一些可對(duì)函數(shù)(指數(shù)學(xué)函數(shù))進(jìn)行操作的函數(shù)〔指Matlab函數(shù)〕，稱(chēng)之為泛函指令〔或泛函函數(shù)〕。使用泛函函數(shù)的步驟：１、數(shù)學(xué)函數(shù)的定義：用字符串、內(nèi)聯(lián)函數(shù)或定義函數(shù)文件的方式，將指定的數(shù)學(xué)函數(shù)通過(guò)Matlab表達(dá)式表示出來(lái)；２、泛函函數(shù)的調(diào)用：用以上定義的數(shù)學(xué)函數(shù)作為參數(shù)，調(diào)用所需的泛函函數(shù)。

4.7.1數(shù)學(xué)函數(shù)在Matlab中的表示１、用字符串表示函數(shù)：調(diào)用方法：直接調(diào)用泛函函數(shù)；注：這種方法比較過(guò)時(shí)，現(xiàn)多已不用。２、用內(nèi)聯(lián)函數(shù)表示函數(shù)：調(diào)用方法：通過(guò)內(nèi)聯(lián)函數(shù)變量名調(diào)用泛函函數(shù)；３、用函數(shù)文件表示函數(shù)：調(diào)用方法：通過(guò)函數(shù)文件名字符串或函數(shù)句柄調(diào)用泛函函數(shù)；4.7.2函數(shù)的極值點(diǎn)MATLAB中用于求極小值的函數(shù)是：1、求單變量函數(shù)fun在區(qū)間(a,b)上的極小值點(diǎn)：x=fminbnd(fun,a,b)%x為極小值點(diǎn)的橫坐標(biāo)或[x,y]=fminbnd(fun,a,b)%x，y分別為極小值點(diǎn)的橫、縱坐標(biāo)2、求多變量函數(shù)F(x)在自變量向量X0附近的極小值點(diǎn)：X=fminsearch(F,X0)%X為極小值點(diǎn)的自變量坐標(biāo)向量或[X,y]=fminsearch(F,X0)%X為極小值點(diǎn)的自變量坐標(biāo)向量，y為極小值點(diǎn)的因變量坐標(biāo)注：MATLAB沒(méi)有專(zhuān)門(mén)提供求函數(shù)極大值點(diǎn)的函數(shù)，但只要注意到-f(x)在區(qū)間(a,b)上的極小值點(diǎn)與f(x)在(a,b)的極大值點(diǎn)的橫坐標(biāo)相同、縱坐標(biāo)異號(hào)，所以可以用fminbnd(-f,a,b)返回函數(shù)f(x)在區(qū)間(a,b)上的極大值。例4.28求函數(shù)f(x)=x-－+5在區(qū)間(-10，1)和(1，10)上的極小值點(diǎn)。1、使用字符串函數(shù)的解法：f1='x-1/x+5';fminbnd(f1,-10,1)%求函數(shù)在區(qū)間(-10,1)內(nèi)的極小值點(diǎn)。2、使用內(nèi)聯(lián)函數(shù)的解法：f1=inline('x-1/x+5');fminbnd(f1,-10,1)%通過(guò)內(nèi)聯(lián)函數(shù)變量名調(diào)用，函數(shù)名f1不用加'1x3、使用函數(shù)文件的解法：編輯如下函數(shù)文件，并存儲(chǔ)為f2.m:functionf=f2(x)f=x-1/x+5;再在MATLAB命令窗口，輸入命令：fminbnd('f2',1,10)%通過(guò)函數(shù)文件名字符串調(diào)用，求函數(shù)在區(qū)間(1,10)內(nèi)的極小值點(diǎn)或f=@f2;fminbnd(f,1,10)%通過(guò)函數(shù)句柄調(diào)用例4.29設(shè)有函數(shù)f(x,y,z)=x+——+—+—，求函數(shù)f在(0.5,0.5,0.5)附近的極小值。建立函數(shù)文件fxyz.m：functionf=fxyz(u)x=u(1);y=u(2);z=u(3);f=x+y.^2./x/4+z.^2./y+2./z;在MALAB命令窗口，輸入命令：U=fminsearch('fxyz',[0.5,0.5,0.5])%求函數(shù)的極小值點(diǎn)fxyz(U)%求函數(shù)的極小值或fun=@fxyz;[U,f]=fminsearch(fun,[0.5,0.5,0.5])y2x/4z2y2z4.7.3數(shù)值微積分1、數(shù)值微分MATLAB中，沒(méi)有直接提供求數(shù)值導(dǎo)數(shù)的函數(shù)，只有計(jì)算向前差分的函數(shù)。Dy=diff(y)計(jì)算向量y的向前差分,Dy(i)=y(i+1)-y(i),0<i<n。Dy=diff(y,n)計(jì)算y的n階向前差分，diff(y,2)=diff(diff(y))。Dy=diff(A,n,dim)計(jì)算矩陣A的n階差分，dim=1時(shí)(缺省狀態(tài)),按列計(jì)算差分;dim=2，按行計(jì)算差分。注：這局部?jī)?nèi)容不屬于“泛函函數(shù)〞，是出于內(nèi)容的相關(guān)才安排在一起。另外，diff函數(shù)除了在數(shù)值計(jì)算中用于計(jì)算差分，還可在符號(hào)計(jì)算中用于求符號(hào)導(dǎo)數(shù)自變量x的向量函數(shù)y的值的向量函數(shù)y的導(dǎo)函數(shù)的值的向量即“數(shù)值導(dǎo)數(shù)〞，由導(dǎo)數(shù)的定義，可以通過(guò)差分來(lái)計(jì)算y=f(x):函數(shù)的符號(hào)表達(dá)式y(tǒng)=f'(x):導(dǎo)函數(shù)的符號(hào)表達(dá)式，即“符號(hào)導(dǎo)數(shù)〞例4.30求向量sin(X)的1~3階差分。設(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)例4.31用不同的方法求如下函數(shù)的數(shù)值導(dǎo)數(shù)，并在同一個(gè)坐標(biāo)系中做出f'(x)的圖象。程序如下：f=inline('sqrt(x.^3+2*x.^2-x+12)+(x+5).^(1/6)+5*x+2');%用內(nèi)聯(lián)函數(shù)給出函數(shù)表達(dá)式，注意變量的點(diǎn)運(yùn)算的運(yùn)用x=-3:0.01:3;%定義自變量向量%用擬合多項(xiàng)式的方法求數(shù)值導(dǎo)數(shù)，用于擬合的多項(xiàng)式的次數(shù)越小，誤差越大p=polyfit(x,f(x),5);%用5次多項(xiàng)式p擬合f(x)dp=polyder(p);%對(duì)擬合多項(xiàng)式p求導(dǎo)數(shù)dpdpx=polyval(dp,x);%求dp在假設(shè)點(diǎn)的函數(shù)值%用導(dǎo)數(shù)的定義，即通過(guò)求因變量和自變量的差分來(lái)求導(dǎo)數(shù)dx=diff(f([x,3.01]))/0.01;%直接對(duì)f(x)求數(shù)值導(dǎo)數(shù)(x0,y0)(x1,y1)⊿y=y1-y0因變量的差分dx=⊿x=x1-x0自變量的差分dyy'|x=x0=－≈－－dydx⊿y⊿x%直接用內(nèi)聯(lián)函數(shù)給出導(dǎo)函數(shù)的表達(dá)式，將自變量向量代入以求導(dǎo)g=inline('(3*x.^2+4*x-1)./sqrt(x.^3+2*x.^2-x+12)/2+1/6./(x+5).^(5/6)+5');gx=g(x);%求函數(shù)f的導(dǎo)函數(shù)g在采樣點(diǎn)的導(dǎo)數(shù)plot(x,dpx,x,dx,'g.',x,gx,'r-');%作圖2、數(shù)值積分(1)被積函數(shù)是一個(gè)解析式quad(f,a,b,tol