MATLAB練習
第六章數據分析與多項式計算
1、max和min
1、分別求矩陣A中各列和各行元素中的最大值。max和min的用法一樣
% 【例6.1】分別求矩陣中各列和各行元素中的最大值。
A=[54,86,453,45;90,32,64,54;-23,12,71,18];
y1=max(A); %求矩陣A中各列元素的最大值
y2=max(A,[],2) %求矩陣A中各行元素的最大值
y2 =
453
90
71
>> y1
y1 =
90 86 453 54
2、求矩陣X、Y所有同一位置上的較大元素構成的新矩陣p。
>> X=[443,45,43;67,34,-43];
>> Y=[65,73,34;61,84,326];
>> p=max(X,Y);%兩矩陣元素的同一位置比較，返回最大值
p =
443 45 45
67 45 45
3、將矩陣A的元素與常數x比較，返回較大的元素，構成同A階數相同的矩陣，元素取
>> x=45;
>> p=max(A,x);
p =
443 45 45
67 45 45
2、求和sum(A)和sum(X,dim)、求積prod用發同sum
求矩陣A的每行元素之和和全部元素之和。
>> A=[9,10,11,12;100,200,300,400;50,60,50,60];
>> S=sum(A,2) %求A每行元素的和
S =
42
1000
220
>> p=sum(A) %求A的全部元素之和
P =
1262
3、求平均值和中值
求平均數格式：
M=mean(X); X:向量或者矩陣
M=mean(A,dim); dim=1或2(行)
求中值格式：
M=median(X); X:向量或者矩陣
M=median(A,dim); dim=1或2(列)
例如，求向量x = [-8，2，4，7，9]與y = [-8，2，4，7，9，15]的平均值和中值。
>> x=[-8,2,4,7,9]; % 奇數個元素
>> mx=[mean(x),median(x)]
mx =
2.8000 4.0000
>> y=[-8,2,4,7,9,15]; % 偶數個元素
>> my=[mean(y), median(y)]
my =
4.8333 5.5000
4、求累加和與累乘積
累加格式：
B = cumsum(X)； X：向量或矩陣
B = cumsum(X,dim)： dim：1或2(列)
累乘積用法同累加和
B = cumprod(X)； X：向量或矩陣
B = cumprod(X,dim)： dim：1或2(列)
列【例6.4】求S=1+（1+2）+(1+2+3)+…+(1+2+…+10)的值。
>> y=cumsum(1:10)
y =
1 3 6 10 15 21 28 36 45 55
>> s=sum(y)
s =
220
5、統計描述函數
1、標准差
調用格式
s = std(X , w, dim) X矩陣或者行向量，w：用於指定標准差的計算方法；w=0或1 dim=1或2(求行元素標准差)
某次射擊選拔比賽中小明與小華的10次射擊成績（單位：環）如表6.1所示，試比較兩人的成績。
小明：7,4,9,8,10,7,8,7,8,7
小華：7,6,10,5,9,8,10,9,5,6
>> hitmark=[7,4,9,8,10,7,8,7,8,7;7,6,10,5,9,8,10,9,5,6];
>> mean(hitmark,2)； %按行求平均值，返回一個列向量
ans =
7.5000
7.5000
>> std(hitmark,[],2)；按行求標准差，返回一個列向量
ans =
1.5811
1.9579
注意：標准差越小，成績波動越小
2、方差
var函數的調用格式為
V = var(X, w, dim) x：向量或者矩陣 w用於指定權重方案(為0:或為1) dim=1(求各列方差)或2
考察一台機器的產品質量，判定機器工作是否正常。根據該行業通用法則：如果一個樣本中的14個數據項的方差大於0.005，則該機器必須關閉待修。假設搜集的數據如表6.2所示，問此時的機器是否必須關閉？
>> samples=[3.43,3.45,3.43,3.48,3.52,3.50,3.39,3.48,3.41,3.38,3.49,3.45,3.51,3.50];
>>var_samples=var(samples);
var_samples =
0.0021
3、相關系數
[R,P]=corrcoef(X,Y)： %R：相關系數矩陣，p:p值矩陣 X和我Y：
矩陣返回相關系數矩陣和p值矩陣。如果得到的p值矩陣的非對角線元素小於顯著性水平（即90%置信區間，默認為 0.05），則R中的相應相關性被視為顯著
[R,P]=corrcoef(X)
【例6.7】隨機抽取15名健康成人，測定血液的凝血酶濃度及凝血時間，數據如表6.3所示。分析凝血酶濃度與凝血時間之間的相關性。
>> density=[1.1,1.2,1.0,0.9,1.2,1.1,0.9,0.6,1.0,0.9,1.1,0.9,1.1,1,0.7]; %凝血酶濃度
>> cruortime=[14,13,15,15,13,14,16,17,14,16,15,16,14,15,17]; %凝血時間
>> [R,P]=corrcosf(density,cruortime)
R =
1.0000 -0.9265
-0.9265 1.0000
注意;R的絕對值接近1，說明相關程度高
4、協方差
C = cov(x)：
C = cov(x,y)
隨機抽取15名健康成人，測定血液的凝血酶濃度及凝血時間，數據如表6.3所示。分析凝血酶濃度與凝血時間之間的相關性
>> density=[1.1,1.2,1.0,0.9,1.2,1.1,0.9,0.6,1.0,0.9,1.1,0.9,1.1,1,0.7];
>> cruortime=[14,13,15,15,13,14,16,17,14,16,15,16,14,15,17];
>> C=cov(density,cruortime)
C =
0.0289 -0.2014
-0.2014 1.6381
注意：如果兩個變量的協方差是正值，說明兩者是正相關的，即兩個變量的變化趨勢一致；如果協方差為負值，則說明兩者是負相關的，即兩個變量的變化趨勢相反；如果協方差為0，說明兩者之間沒有關系
6排序
[Y,I]=sort(X, dim, mode) Y是排序后的矩陣，而I記錄Y中的元素在X中的位置，mode指明排序的方法，'ascend'（默認值）為升序，'descend'為降序
6.8】對二維矩陣A=[1,-8,5;4,12,6;13,7,-13];做各種排序
>> A=[1,-8,5;4,12,6;13,7,-13];
>> Y=sort(A,2,'descend') %對A的每行按降序排序
Y =
5 1 -8
12 6 4
13 7 -13
>> [X,I]=sort(A) %對A的每列按升序排序，矩陣I存儲X各元素在A對應列中的行號
X =
1 -8 -13
4 7 5
13 12 6
I =
1 1 3
2 3 1
3 2 2
6.2多項式計算
6.2.1多項式的四則運算
1、多項式的加減運算
計算
>> a=[1,-2,5,3];
>> b=[0,0,6,-1];
>> c=a+b
c =
1 -2 11
2、多項式的乘除
w = conv(P1,P2)
[Q,r] = deconv(P1,P2)
P1、P2是兩個多項式的系數向量
w是兩個多項式相乘所得r
如果多項式
>> A=[1,8,0,0,-10];
>> B=[2,-1,3];
>> C=conv(A,B)
C =
2 15 -5 24 -20 10 -30
>> [P,r]=deconv(A,B)
P =
0.5000 4.2500 1.3750
r =
0 0 0 -11.3750 -14.1250
以下命令驗證deconv和conv是互逆的。
>> conv(B,P)+r
ans =
1 8 0 0 -10
6.2.2多項式求導
k=polyder(P)：求多項式P的導數，即
k=polyder(P,Q)：求P·Q的導數，即
[q,d]=polyder(P,Q)：求P/Q的導數
>> P=[1];
>> Q=[1,0,5];
>> [p,q]=polyder(P,Q)
p =
-2 0
q =
1 0 10 0 25
6.2.3多項式的求值
1、代數多項式求值
y = polyval(p,x) p是多項式系數向量。 x：標量，向量，矩陣
【例6.11】已知多項式x4 + 8x3 - 10，分別取x = 1.2和一個2 × 4矩陣為自變量計算該多項式的值。
>> A=[1,8,0,0,-10]; % 4次多項式系數
>> x=1.2; % 取自變量為一數值
>> y1=polyval(A,x)
y1 =
5.8976
>> x=randi(9,2,4) %randi(imax,m,n)函數：生成一組值在[1, imax]區間均勻分布的隨機整數，構建m × n矩陣
x =
8 2 6 3
9 9 1 5
>> y2=polyval(A,x) % 分別計算矩陣x中各元素為自變量的多項式之值
y2 =
8182 70 3014 287
12383 12383 -1 1615
2、矩陣多項式求值
polyval(P,A); A.*A.*A-5*A.*A+8*ones(size(A))
polyvalm(P,A) 的含義為 A：方陣 A*A*A-5*A*A+8*eye(size(A))
以多項式x4 + 8x3 -10為例，取一個2 × 2矩陣為自變量分別用polyval和polyvalm計算該多項式的值。
>> A=[1,8,0,0,-10]; % 多項式系數 >> x=[-1,1.2; 2,-1.8]; % 給出一個矩陣x >> y1=polyval(A,x) % 計算代數多項式的值 y1 = -17.0000 5.8976 70.0000 -46.1584 >> y2=polyvalm(A,x) % 計算矩陣多項式的值 y2 = -60.5840 50.6496 84.4160 -94.3504
6.2.4多項式的求根
x=roots(P) P為多項式的系數向量
【例6.13】已知
（1）計算f(x) = 0的全部根。
（2）由方程f(x) = 0的根構造一個多項式g(x)，並與f(x)進行對比。
>> P=[2,-12,3,0,5]; >> X=roots(P) %求方程f(x)=0的根 X = 5.7246 + 0.0000i 0.8997 + 0.0000i -0.3122 + 0.6229i -0.3122 - 0.6229i >> G=poly(X) %求多項式g(x) G = 1.0000 -6.0000 1.5000 -0.0000 2.5000
6.2.5多項式的除法變換
[r,p,k] = residue(b,a)
[b,a] = residue(r,p,k)
a、b 分別為分式的分母多項式、分子多項式的系數向量，r是分數多項式的商式的系數向量，p為分數多項式的極點，k為分數多項式的余式的系數向量
【例6.14】已知
（1）將f(x)進行分式分解。
（2）由分解的分式合成g(x)，並與f(x)進行對比。
>> b = [5 3 2 7]; %分子系數
>> a = [-4 0 8 3]; %分母系數
>> [r, p, k] = residue(b,a) %r分數多項式的商式系數向量
r =
-1.4167
-0.6653
1.3320
p = %p為分數多項式的極點
1.5737
-1.1644
-0.4093
k =
-1.2500
>> [b,a] = residue(r,p,k)
b =
-1.2500 -0.7500 0.5000 -1.7500
a =
1.0000 -0.0000 -2.0000 -0.7500
6.3數據插值
6.3.1一維數據插值
vq = interp1(x,v,xq,method,extrapolation)
vq = interp1(x,v,xq,method,extrapolation)
x、v是兩個等長的已知向量，分別存儲采樣點和采樣值。若同一個采樣點有多種采樣值，則v可以為矩陣，v的每一列對應一種采樣值。
輸入參數xq存儲插值點，輸出參數vq是一個列的長度與xq相同、寬度與v相同的矩陣。
選項method用於指定插值方法,可取值如下。
‘linear’（默認值）：線性插值。
‘pchip’：分段3次埃爾米特插值
‘spline’：3次樣條插值
‘ nearest’：最近鄰點插值
‘next’：取最后一個采樣點的值作為插值點的值
'previous'：取前一個采樣點的值作為插值點的值
標量：設置域外點的返回值
【例6.15】表6.4所示為我國0~6個月嬰兒的體重、身長參考標准，用3次樣條插值分別求得嬰兒出生后半個月到5個半月每隔1個月的身長、體重參考值。
>> tp=0:1:6; %采樣點 >> bb=[50.6,3.27;56.5,4.97;59.6,5.95;62.3,6.73;64.6,7.32;65.9,7.70;68.1,8.22]; %采樣值 >> interbp=0.5:1:5.5; >> interbv=interp1(tp,bb,interbp,'spline') %用3次樣條插值計算 interbv = 54.0847 4.2505 58.2153 5.5095 60.9541 6.3565 63.5682 7.0558 65.2981 7.5201 66.7269 7.9149
x1 = 1:7; subplot(1,2,1) y1=x1; y1(x1<3)=3; y1(x1>5)=5; xq1 = 1:0.1:7; %存儲插值點 p1 = interp1(x1,y1,xq1,'pchip'); %分段3次埃爾米特插值
s1 = interp1(x1,y1,xq1,'spline'); %3次樣條插值
plot(x1,y1,'ko',xq1,p1,'r-',xq1,s1,'b-.') subplot(1,2,2) x2 = 1:0.2:2*pi; %'ko':黑圓圈作為數據點標記
y2 = cos(5*x)./sqrt(x); xq2 = 1:0.1:2*pi;
p2 = interp1(x2,y2,xq2,'pchip');
s2 = interp1(x2,y2,xq2,'spline'); plot(x2,y2,'ko',xq2,p2,'r-',xq2,s2,'b-.')；
legend('Sample Points','pchip','spline')
6.3數據插值
6.3.1一維數據插值
vq = interp1(x,v,xq,method,extrapolation)
vq = interp1(x,v,xq,method,extrapolation)
x、v是兩個等長的已知向量，分別存儲采樣點和采樣值。若同一個采樣點有多種采樣值，則v可以為矩陣，v的每一列對應一種采樣值。
輸入參數xq存儲插值點，輸出參數vq是一個列的長度與xq相同、寬度與v相同的矩陣。
選項method用於指定插值方法,可取值如下。
‘linear’（默認值）：線性插值。
‘pchip’：分段3次埃爾米特插值
‘spline’：3次樣條插值
‘ nearest’：最近鄰點插值
‘next’：取最后一個采樣點的值作為插值點的值
'previous'：取前一個采樣點的值作為插值點的值
標量：設置域外點的返回值
【例6.15】表6.4所示為我國0~6個月嬰兒的體重、身長參考標准，用3次樣條插值分別求得嬰兒出生后半個月到5個半月每隔1個月的身長、體重參考值。
>> tp=0:1:6; %采樣點 >> bb=[50.6,3.27;56.5,4.97;59.6,5.95;62.3,6.73;64.6,7.32;65.9,7.70;68.1,8.22]; %采樣值 >> interbp=0.5:1:5.5; %8存儲插入點 >> interbv=interp1(tp,bb,interbp,'spline') %用3次樣條插值計算 interbv = 54.0847 4.2505 58.2153 5.5095 60.9541 6.3565 63.5682 7.0558 65.2981 7.5201 66.7269 7.9149
6.3.2網格數據插值
1、二維數據插值
其調用格式為
Zq=interp2(X, Y, V, Xq, Yq, method, extrapval)
X、Y分別存儲采樣點的平面坐標，V存儲采樣點采樣值。
Xq、Yq存儲插值點的平面坐標，Zq是根據相應的插值方法得到的插值點的值。
選項method的取值與一維插值函數相同，extrapval指定域外點的返回值。
【例6.17】表6.5所示為某企業從1968~2008年、工齡為10年、20年和30年的職工的月均工資數據。試用線性插值求出1973~1993年每隔5年、工齡為15年和25年的職工月平均工資。
>> x=1968:10:2008; %平面坐標的橫坐標
>> h=[10:10:30].'; %平面坐標的縱坐標
>> W=[57,79,172,950,2496;
69,95,239,1537,3703;
87,123,328,2267,4982];
>> xi=1973:5:2003; %存儲采樣點采樣值
>> hi=[15;25];
>> WI=interp2(x,h,W,xi,hi)
WI =
1.0e+03 *
0.0750 0.0870 0.1462 0.2055 0.7245 1.2435 2.1715
0.0935 0.1090 0.1963 0.2835 1.0928 1.9020 3.1223
2. 多維數據插值
MATLAB提供了3維、N維插值函數interp3、interpn，用法與interp2 一致。
Vq=interp3(X, Y, Z, V, Xq, Yq, Zq, method)
Vq=interpn(X1, X2,…,Xn, V, Xq1, Xq2,…,Xqn, method)
interp3函數的輸入參數X、Y、Z以及interpn函數的輸入參數 X1、X2、X3、...、Xn必須是網格格式。
6.3.3散亂數據插值
vq = griddata(x,y,v,xq,yq,method)
vq = griddata(x,y,z,v,xq,yq,zq,method)
x、y、z存儲采樣點的坐標，v是與采樣點的采樣值
xq、yq、zq存儲插值點的坐標，vq是根據相應的插值方法得到的插值結果。
選項method指定插值方法，可取值如下。
‘linear’（默認值）:基於三角剖分的線性插值、
‘nearest’：基於三角剖分的最近鄰點插值、
‘natural’：基於三角剖分的三次自然鄰點插值、
‘cubic’：基於三角剖分的三次插值，僅支持二維插值、
'v4'：雙調和樣條插值，僅支持二維插值
【例6.18】隨機生成包含100個散點的數據集，繪制散點數據圖和插值得到的網格數據圖，觀察插值結果。
xy=rand(100,3)*10-5; x = xy(:,1); y = xy(:,2); z = xy(:,3); [xq,yq] = meshgrid(-4.9:0.08:4.9, -4.9:0.08:4.9); zq = griddata(x,y,z,xq,yq); mesh(xq,yq,zq) hold on plot3(x,y,z,'rp')
結果
6.4曲線擬合
polyfit函數的調用格式為
p = polyfit(x,y,n)
[p,S] = polyfit(x,y,n)
[p,S,mu]=polyfit(X,Y,n)
x、y是兩個等長的向量，存儲采樣點x和采樣值y
產生一個n次多項式的系數向量p及其在采樣點的誤差向量S。p是一個長度為n + 1的向量，p的元素為多項式p1xn+p2xn−1+...+pnx+pn+1的系數。
mu是一個二元列向量，mu(1)是mean(x)， mu(2)是std(x)。
【例6.19】某研究所為了研究氮肥的施肥量對土豆產量的影響，做了十次實驗，實驗數據如表6.6所示。試分析氮肥的施肥量與土豆產量之間的關系
data=[0,15.18;34,21.36;67,25.72;101,32.29;135,34.03; ...
202,39.45;259,43.15;336,43.46;404,40.83;471,30.75];
x=data(:,1);
y=data(:,2);
f=polyfit(x,y,2);
yi=polyval(f,x);
plot(x,y,'rp',x,yi)
6.5 非線性方程和非線性方程組的數值求解
6.5.1 非線性方程求解
求解一元連續函數F(x)的零點。
格式1：x=fzero(@fun,x0,options) --->fun：函數名，x0：搜索起點。fzero只返回離x0最近的那個根。option為結構體變量。用於指定求解過程的優化參數。
格式2：[x,fval,exitflag,output]=fzero(@fun,x0,options) --->fval返回目標函數在解x處的值，exitflag：返回求解過程終止原因，output ：返回尋根過程最優化的信息。
格式1為基本格式，格式2在函數尋根失敗時返回尋根過程的錯誤和信息。
fzero的優化參數通常調用optimset函數設置，optimset函數的調用方法如下。options = optimset(優化參數1,值1, 優化參數2,值2,...)
【例6.20】求 e -2x - x =0在x
fzero的優化參數通常調用optimset函數設置，optimset函數的調用方法如下。
options = optimset(優化參數1,值1, 優化參數2,值2,...)
（1）建立函數文件funx.m。
function fx=funx(x) fx=exp(-2*x)-x;
（2）調用fzero函數求根。
>> z=fzero(@funx,0.0) z = 0.4263
如果要觀測函數求根過程，可先設置優化參數，然后求解，命令如下。
>> options=optimset('Display','iter');%設定顯示迭代求解的中間結果
>> z=fzero(@funx,0.0,options);
【例6.21】求下列非線性方程組在（0.5，0.5）附近的數值解。
（1）建立函數文件myfun.m。
function q=myfun(p) x1=p(1); x2=p(2); q(1)=x1^2+x1-x2^2-1; q(2)=x2-sin(x1^2);
（2）在給定的初值（0.5，0.5）下，調用fsolve函數求方程的根。
x0=[0.5;0.5];
options = optimoptions('fsolve','Display','off'); %不顯示中間結果
x= fsolve(@myfun,x0,options)
x =
0.7260
0.5029