第6講：MATLAB數(shù)據(jù)分析與應(yīng)用

第6講：MATLAB數(shù)據(jù)分析與應(yīng)用6.1多項(xiàng)式運(yùn)算6.2線性方程組的求解6.3非線性方程組的數(shù)值求解6.4函數(shù)的零點(diǎn)和極值6.5基本統(tǒng)計(jì)處理6.6數(shù)據(jù)插值和曲線擬合16.1多項(xiàng)式運(yùn)算matlab語(yǔ)言把多項(xiàng)式表達(dá)成一個(gè)行向量，該向量中的元素是按多項(xiàng)式降冪排列的。f(x)=anxn+an-1xn-1+……+a0

可用行向量p=[anan-1……a1a0]表示。1.poly(x)——x是矩陣或向量，產(chǎn)生首元素是1的n+1維的特征多項(xiàng)式的系數(shù)向量。2例1:

a=[12.1229-5.7345-0.3884];p=poly(a)

p=1.00-6.00-72.00-27.00例2:

a=[123;456;780];p=poly(a)

p=1.00-6.00-72.00-27.00友矩陣：poly(compan(p))3p是多項(xiàng)式p(x)=x3-6x2-72x-27的matlab描述方法，可用：p1=poly2str(p,'x')—函數(shù)，顯示多項(xiàng)式的數(shù)學(xué)形式：

p1=x^3-6x^2-72x-27matlab指令為

p=[1,-6,-72,-27];

p1=poly2str(p,'x')4例3:

a=[123;456;780];p=poly(a)

p=1.00-6.00-72.00-27.00r=roots(p)

r=12.12-5.73——顯然r是矩陣a的特征值-0.392.roots(A)——A是多項(xiàng)式的系數(shù)，求多項(xiàng)式的根

5當(dāng)然我們可用poly令其返回多項(xiàng)式形式p2=poly(r)

p2=1.00-6.00-72.00-27.00

matlab規(guī)定多項(xiàng)式系數(shù)向量用行向量表示，一組根用列向量表示。6Y=polyval(P,x)：代數(shù)多項(xiàng)式求值，若x為一數(shù)值，則求在該點(diǎn)的值；若x為向量或矩陣，則對(duì)向量或矩陣中的每個(gè)元素求其值。polyvalm(P,A)：矩陣多項(xiàng)式求值，A為方陣，它以方陣為自變量求多項(xiàng)式的值。設(shè)A為方陣，P代表多項(xiàng)式x3-5x2+8，則后者含義是：A^3-5*A^2+8*eye(size(A))，而前者是：A.^3-5*A.^2+8*ones(size(A))。實(shí)驗(yàn)：p=[1,-5,0,8];A=compan(p);polyvalm(p,A)3.求多項(xiàng)式的值7例4:

a(x)=x2+2x+3;b(x)=4x2+5x+6;c=(x2+2x+3)(4x2+5x+6)a=[123];b=[456];c=conv(a,b)=conv([123],[456])

c=4.0013.0028.0027.0018.00p=poly2str(c,'x')

p=4x^4+13x^3+28x^2+27x+184.conv(),convs()——多項(xiàng)式乘運(yùn)算8a=[1,2,3];c=[4.00,13.00,28.00,27.00,18.00]d=deconv(c,a)

d=4.005.006.00[d,r]=deconv(c,a)余數(shù)c除a后的整數(shù)5.deconv()——多項(xiàng)式除運(yùn)算9命令格式：

polyder(p):求p的微分polyder(a,b):求多項(xiàng)式a,b乘積的微分[p,q]=polyder(a,b):求多項(xiàng)式a,b商的微分例5：a=[12345];poly2str(a,'x')

ans=x^4+2x^3+3x^2+4x+5b=polyder(a)

b=4664poly2str(b,'x')

ans=4x^3+6x^2+6x+46.polyder函數(shù)——多項(xiàng)式微分10

1.利用左除運(yùn)算符的直接解法對(duì)于方程ax=b，a為an×m矩陣，有三種情況：當(dāng)n=m時(shí)，此方程成為“恰定”方程當(dāng)n>m時(shí)，此方程成為“超定”方程當(dāng)n<m時(shí)，此方程成為“欠定”方程matlab定義的左除運(yùn)算可以很方便地解上述三種方程。6.2線性方程組的求解11方程ax=b(a為非奇異)，x=a-1b矩陣逆兩種解:x=inv(a)

b—求逆運(yùn)算解方程x=a\b—左除運(yùn)算解方程（等價(jià)）例6:x1+2x2=82x1+3x2=13a=[12;23];b=[8;13];

x=inv(a)*b

x=a\bx=x=2.002.003.003.00

1）“恰定”方程組12方程ax=b,m<n時(shí)，不存在唯一解。方程解(a'a)x=a'b,x=(a'a)-1a'b——求逆法x=a\b——最小二乘法

(等價(jià))例7:

x1+2x2=12x1+3x2=23x1+4x2=3，a=[12;23;34];b=[1;2;3];

x=a\b

x=inv(a'

a)

a'

b

x=x=1.001.0000.00

2）“超定”方程組13方程ax=b,m>n時(shí)，即不定情況，有無(wú)窮多個(gè)解存在。matlab可求出兩個(gè)解：1）具有最多零元素的解——除法求解2）具有最小長(zhǎng)度或范數(shù)的解——偽逆pinv例8:x1+2x2+3x3=12x1+3x2+4x3=2，a=[123;234];b=[1;2];

x=a\b

x=pinv(a)

b

x=1.00x=0.8300.330-0.173）“欠定”方程組14

矩陣分解是指根據(jù)一定的原理用某種算法將一個(gè)矩陣分解成若干個(gè)矩陣的乘積。

(1)LU分解：將非奇異的矩陣分解為一個(gè)下三角矩陣L和一個(gè)上三角矩陣U的乘積，即A=L*U。Ax=b可變換為x=U\(L\b)或x=U\(L\Pb)。

分解格式(矩陣A必須是方陣)：[L,U]=lu(A)：產(chǎn)生一個(gè)上三角陣U和一個(gè)下三角陣L。[L,U,P]=lu(A)：產(chǎn)生一個(gè)上三角陣U和一個(gè)下三角陣L以及一個(gè)置換矩陣P2.利用矩陣分解的直接解法15命令如下：A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4];b=[13,-9,6,0]';[L,U]=lu(A);x=U\(L\b)或采用LU分解的第2種格式，命令如下：[L,U,P]=lu(A);x=U\(L\P*b)例9

用LU分解求解例7-1中的線性方程組。16

(2)QR分解：把矩陣A分解為一個(gè)正交矩陣Q和一個(gè)上三角矩陣R的乘積，即A=Q*R。Ax=b可變換為x=R\(Q\b)或x=E(R\(Q\b))。

分解格式：[Q,R]=qr(A)：產(chǎn)生一個(gè)正交矩陣Q和一個(gè)上三角矩陣R。[Q,R,E]=qr(A)：產(chǎn)生一個(gè)正交矩陣Q、一個(gè)上三角矩陣R以及一個(gè)置換矩陣E。QR分解17(3)Cholesky分解：把對(duì)稱正定的矩陣A分解成一個(gè)上三角矩陣R和其轉(zhuǎn)置矩陣R'的乘積，即A=R'R。Ax=b可變成x=R\(R'\b)。分解格式：R=chol(A)：產(chǎn)生一個(gè)上三角陣R。[R,p]=chol(A)：A對(duì)稱正定時(shí)，p=0；否則p為一個(gè)正整數(shù)。如果A滿秩，則R為一個(gè)滿足R‘R=A(1:p-1,1:p-1)的階為p-1的上三角陣。Cholesky分解18@函數(shù)Matlab6.0以后的版本中可以用@函數(shù)進(jìn)行函數(shù)調(diào)用，其調(diào)用格式：handle=@function這類似于C++語(yǔ)言中的引用。實(shí)例：利用句柄傳遞數(shù)據(jù)。fhandle=@humps;%為hump函數(shù)創(chuàng)建一個(gè)函數(shù)句柄，并將它指定為fhandle變量x=fminbnd(@humps,0,1)%傳遞humps句柄給fminbnd函數(shù)，傳遞所有變量匿名函數(shù)19對(duì)于非線性方程組F(X)=0，用fsolve函數(shù)求其數(shù)值解。fsolve函數(shù)的調(diào)用格式為：X=fsolve('fun',X0)

X=fsolve('fun',X0,options)其中X為返回的解，fun是用于定義需求解的非線性方程組的函數(shù)文件名，X0是初值，options為選項(xiàng)設(shè)定，默認(rèn)缺省。如果想改變其中某個(gè)選項(xiàng)，則可以調(diào)用optimset()函數(shù)來(lái)完成。6.3非線性方程組的求解20例10：求解方程在[-1,10]內(nèi)的解。

首先作圖如下：

fplot(‘[5*x^2*sin(x)-exp(-x),0]’,[0,10])發(fā)現(xiàn)在[0,10]區(qū)間中有4個(gè)解，分別在0,3,6,9附近，所以用命令：

fun=inline(‘5*x.^2.*sin(x)-exp(-x)’);fsolve(fun,[0369],1e-6)

ans=0.50183.14076.28329.4248注意：用inline建立的函數(shù)fun調(diào)用時(shí)不加單引號(hào)；而用M函數(shù)文件建立的函數(shù)則要加單引號(hào)或用@函數(shù)。21

首先編制函數(shù)文件fc.m如下：functiony=fc(x)y=[x(1)-0.7*sin(x(1))-0.2*cos(x(2));x(2)-0.7*cos(x(1))+0.2*sin(x(2))];

在MATLAB命令窗口輸入：x0=[0.50.5];fsolve(‘fc’,x0)或fsolve(@fc,x0)ans=0.52650.5079例11：求非線性方程組在(0.5,0.5)附近的解。22對(duì)于函數(shù)f(x)的零點(diǎn)，用fzero函數(shù)求解?；靖袷綖椋簒=fzero('fun',X0)[x,fval]=fzero(fun,X0,options)求函數(shù)零點(diǎn)x和該點(diǎn)的函數(shù)值fval，fun可以是字符串、內(nèi)聯(lián)函數(shù)或M函數(shù)文件名的函數(shù)句柄。X0是初值，options為選項(xiàng)設(shè)定，默認(rèn)缺省。6.4函數(shù)的零點(diǎn)23例12：求的零點(diǎn)。首先作圖確定零點(diǎn)范圍：x=-5:0.1:5;f=x.^3-2*x-5;plot(x,f,x,0);發(fā)現(xiàn)零點(diǎn)大致在2附近，編寫M文件，函數(shù)名zero1.m

functionf=zeros(x)f=x.^3-2*x-5;計(jì)算2附近的零點(diǎn)和該點(diǎn)的函數(shù)值[xz,yz]=fzero(@zero1,2)

xz=yz=2.0946-8.8818e-01624對(duì)于函數(shù)f(x)的極值，用fmidbnd函數(shù)求解極小值。基本格式為：

x=fmidbnd('fun',x1,x2)[x,fval]=fmidbnd(fun,x1,x2,options)求函數(shù)fun在(x1,x2)間的極小點(diǎn)x及該點(diǎn)函數(shù)值fval，用法同fzero。6.5函數(shù)的極值25例13：求的極小值。首先作圖確定極小值范圍：x=-5:0.1:5;f=(x-3).^2-1;plot(x,f);發(fā)現(xiàn)在(2,4)間有極小值，編M文件，函數(shù)名為m1.m

functionf=m1(x)

f=(x-3).^2-1;計(jì)算(2,4)間的極小點(diǎn)和該點(diǎn)的函數(shù)值[xm,ym]=fminbnd(‘m1’,2,4)xm=ym=3-1261、查取最大值——max函數(shù)[Y,I]=max(X)：返回矩陣X的各列中的最大元素值及其該元素的位置并賦予行向量Y與I；當(dāng)X為向量時(shí)，則Y與I為單變量。[Y,I]=max(X,[],DIM)：按數(shù)組X的第DIM維的方向查取其最大的元素值及其該元素的位置賦予向量Y與I。6.6基本統(tǒng)計(jì)處理27U=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同型的向量或矩陣，U的每個(gè)元素等于A對(duì)應(yīng)元素和n中的較大者。28例14

查找下面數(shù)列x的最大值。

x=[359618]%產(chǎn)生數(shù)列xx=359618

y=max(x)%查出數(shù)列x中的最大值賦予yy=9

[y,l]=max(x)%查出數(shù)列x中的最大值及其該元素的位置賦予y,ly=9l=329例15

分別查找下面3×4的二維數(shù)組x中各列和各行元素中的最大值。

x=[1842;9625;3671]%產(chǎn)生二維數(shù)組xx=184296253671

y=max(x)%查出二維數(shù)組x中各列元素的最大值產(chǎn)生賦予行向量yy=987530

[y,l]=max(x)%查出二維數(shù)組x中各列元素的最大值，元素的行下標(biāo)賦予y,ly=9875l=2132[y,l]=max(x,[],1)%本命令的結(jié)果同上

y=9875l=2132

[y,l]=max(x,[],2)%DIM=2，在各行中查找y=8l=29173

31例16

試取下面兩個(gè)2×3的二維數(shù)組x、y所有同一位置上的元素值大者構(gòu)成一個(gè)新矩陣p。

x=[456;148];%產(chǎn)生二維數(shù)組x

y=[175;457];%產(chǎn)生二維數(shù)組y

p=max(x,y)%在x,y同一位置上的兩個(gè)元素中查找出最大值，賦予與x,y同樣大小的二維數(shù)組pp=476458322、查取最小值MIN函數(shù)用來(lái)查取數(shù)據(jù)序列的最小值。3、求中值Y=median(X)：將median(X)返回矩陣X各列元素的中值；若X為向量，則Y為單變量。Y=median(X,DIM)：按數(shù)組X的第DIM維方向的元素求其中值賦予向量Y。若DIM=1，為按列操作；若DIM=2，為按行操作。若X為二維數(shù)組，則Y為一個(gè)向量；若X為一維數(shù)組，則Y為單變量。33例17

對(duì)二維數(shù)組x試從不同維方向求出其中值。x=[1842;9625;3671];y0=median(x)%按列操作y0=3642y1=median(x,1)%按列操作,結(jié)果為行向量y1=3642y2=median(x,2)%按行操作,結(jié)果為列向量y2=3.00005.50004.5000344、求和Y=sum(X)：將sum(X)返回矩陣X各列元素之和賦予行向量Y；若X為向量，則Y為單變量。Y=sum(X,DIM)：按數(shù)組X的第DIM維的方向的元素求其和賦予Y。若DIM=1，為按列操作；若DIM=2，為按行操作。若X為二維數(shù)組，Y為一個(gè)向量；若X為一維數(shù)組，則Y為單變量。355、求平均值Y=mean(X)：將mean(X)返回矩陣X各列元素之和的平均值賦予行向量Y。若X為向量，則Y為單變量。Y=mean(X,DIM)：按數(shù)組X的第DIM維的方向的元素求其平均值賦予向量Y。若DIM=1，為按列操作；若DIM=2，為按行操作。若X為二維數(shù)組，Y為一個(gè)向量；若X為一維數(shù)組，則Y為單變量。36例18x=[456;148];y=sum(x,1)y=5914y=sum(x,2)y=1513y1=mean(x,1)y1=2.504.507.00y2=mean(x,2)y2=5.00004.3333376、求積Y=prod(X)：將prod(X)返回矩陣X各列元素之積賦予行向量Y。若X為向量，則Y為單變量。Y=prod(X,DIM)：按數(shù)組X的第DIM維的方向的元素求其積賦予向量Y。若DIM=1，為按列操作；若DIM=2，為按行操作。若X為二維數(shù)組，Y為一個(gè)向量；若X為一維數(shù)組，則Y為單變量。38例19x=[456;148];y1=prod(x,1)y1=42048y2=prod(x,2)

y2=12032397、累加和、累乘積cumsum(X)：求向量或矩陣X的累加和。cumprod(X)：求向量或矩陣X的累加積。cumsum(A,dim)：當(dāng)dim為1時(shí)，該函數(shù)等同于cumsum(A)；當(dāng)dim為2時(shí)，其第i行是A的第i行的累加和向量。cumprod(A,dim)：同上。408、標(biāo)準(zhǔn)方差std(X)：求向量或矩陣X的標(biāo)準(zhǔn)方差。Y=std(A,flag,dim)：當(dāng)dim=1時(shí)，求各列元素的標(biāo)準(zhǔn)方差；當(dāng)dim=2時(shí)，求各行元素的標(biāo)準(zhǔn)方差。flag取0或1，當(dāng)flag=0時(shí)，按σ=1計(jì)算；當(dāng)flag=1時(shí)，按σ=2計(jì)算。缺省flag=0，dim=1。419、相關(guān)系數(shù)

corrcoef(X)：返回從矩陣X形成的一個(gè)相關(guān)系數(shù)矩陣。此相關(guān)系數(shù)矩陣的大小與矩陣X一樣。它把矩陣X的每列作為一個(gè)變量，然后求它們的相關(guān)系數(shù)。

corrcoef(X,Y)：在這里，X,Y是向量，它們與corrcoef([X,Y])的作用一樣。

42例20

生成滿足正態(tài)分布的10000×5隨機(jī)矩陣，然后求各列元素的均值和標(biāo)準(zhǔn)方差，再求這5列隨機(jī)數(shù)據(jù)的相關(guān)系數(shù)矩陣。命令如下：X=randn(10000,5);M=mean(X)D=std(X)R=corrcoef(X)4310、升序排序sort(X)：返回對(duì)向量X中的元素按升序排列的新向量。[Y,I]=sort(A,dim)：對(duì)矩陣A的各列或各行重新排序，其中dim指明對(duì)A的列還是行進(jìn)行排序。若dim=1，則按列排；若dim=2，則按行排。Y是排序后的矩陣，而I記錄Y中的元素在A中位置。441、一維數(shù)據(jù)插值Y1=interp1(X,Y,X1,‘method’)：根據(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)inear(線性插值)、nearest(最近鄰點(diǎn))、cubic(三次函數(shù))、spline(三次樣條)。6.4數(shù)據(jù)插值和曲線擬合(cftool)45例21

某觀測(cè)站測(cè)得某日6:00時(shí)至18:00時(shí)之間每隔2小時(shí)的室內(nèi)外溫度(℃)，用3次樣條插值分別求得該日室內(nèi)外6:30至17:30時(shí)之間每隔2小時(shí)各點(diǎn)的近似溫度(℃)。h=6:2:18;%時(shí)間變量t=[18,20,22,25,30,28,24;15,19,24,28,34,32,30]’;%已知數(shù)據(jù)XI=6.5:2:17.5%欲插值點(diǎn)YI=interp1(h,t,XI,‘spline’)%3次樣條插值46Z1=interp2(X,Y,Z,X1,Y1,‘method’)：向量X,Y分別是描述兩個(gè)參數(shù)的采樣點(diǎn)，Z是與采樣點(diǎn)對(duì)應(yīng)的函數(shù)值，X1,Y1是兩個(gè)向量或標(biāo)量，描述欲插值的點(diǎn)。Z1是根據(jù)相應(yīng)的插值方法得到的插值結(jié)果。method的取值與一維插值函數(shù)相同。X,Y,Z也可以是矩陣形式。同樣，X1,Y1的取值范圍不能超出X,Y的給定范圍，否則，會(huì)給出“NaN”錯(cuò)誤。2、二維數(shù)據(jù)插值47例22

某實(shí)驗(yàn)對(duì)一根長(zhǎng)10米的鋼軌進(jìn)行熱源的溫度傳播測(cè)試。x表示測(cè)量點(diǎn)0:2.5:10(米)，h表示測(cè)量時(shí)間0:30:60(秒)，T表示測(cè)試所得各點(diǎn)的溫度(℃)。試用線性插值求出在一分鐘內(nèi)每隔20秒、鋼軌每隔1米處的溫度TI。x=0:2.5:10;h=[0:30:60]';T=[95,14,0,0,0;88,48,32,12,6;67,64,54,48,41];xi=[0:10];hi=[0:20:60]';TI=interp2(x,h,T,xi,hi)483、曲線擬合[P,S]=polyfit(X,Y,m)：根據(jù)采樣點(diǎn)X和采樣點(diǎn)函數(shù)值Y，產(chǎn)生一個(gè)m次多項(xiàng)式P及其在采樣點(diǎn)的誤差向量S。其中X,Y是等長(zhǎng)的向量，P是一個(gè)長(zhǎng)