《MATLAB實(shí)踐教程》課件第4章

4.1線性代數(shù)運(yùn)算

4.2多項(xiàng)式運(yùn)算

4.3數(shù)據(jù)插值

4.4數(shù)據(jù)擬合

4.5數(shù)據(jù)分析

4.6數(shù)值積分與數(shù)值微分

4.7解常微分方程和常微分方程組

第4單元MATLAB典型應(yīng)用

4.1.1矢量和矩陣的范數(shù)

矢量x的p范數(shù)表征x的大小。矩陣A的p范數(shù)表征矢量變換Ax范數(shù)與矢量x范數(shù)的最大比值。矢量范數(shù)和矩陣范數(shù)分別定義如下：

4.1線性代數(shù)運(yùn)算norm(x,p)函數(shù)給出矢量x的p范數(shù)，norm(A,p)函數(shù)給出矩陣A的p范數(shù)，程序如下：

x=[20-1-57]

A=[9472;8635;9781]

x_fanshu=[norm(x,1)norm(x)norm(x,inf)]

A_fanshu=[norm(A,1)norm(A)norm(A,inf)]

程序輸出如下：

x=

20-1-57

A=

9472

8635

9781

x_fanshu=

15.00008.88827.0000

A_fanshu=

26.000021.329925.00004.1.2解線性方程組

MATLAB解線性方程組，需先將線性方程組變換為矩陣方程，然后用系數(shù)矩陣、右端項(xiàng)向量、矩陣左除運(yùn)算符“\”或右除運(yùn)算符“/”編程求解。

MATLAB解矩陣方程，系數(shù)矩陣非奇異則給出唯一解；方程個數(shù)大于自變量個數(shù)則給出最小二乘解，此時相當(dāng)于最小二乘擬合；系數(shù)矩陣奇異則給出無窮解或不定解。

問題1：試驗(yàn)檢測了某物理量的時間響應(yīng)，結(jié)果見表4-1。試用延遲函數(shù)擬合表4-1所示的試驗(yàn)數(shù)據(jù)。表4-1時間和響應(yīng)的測定結(jié)果解問題1的程序如下：

clc;clearall;

t=[00.30.81.11.62.3]';y=[0.820.720.630.600.550.50]' %試驗(yàn)數(shù)據(jù)

A=[ones(size(t))exp(-t)];

c=A\y

T=[0:0.1:2.5]';A=[ones(size(T))exp(-T)];Y=A*c;%計(jì)算擬合曲線的數(shù)據(jù)

h=plot(t,y,'ok',T,Y,'--r');boxoff;%繪制試驗(yàn)數(shù)據(jù)點(diǎn)和擬合曲線，觀察比較

set(h,'LineWidth',2,'MarkerSize',8);

xlabel('x','FontName','Times','FontSize',16);

ylabel('y','FontName','Times','FontSize',16);程序執(zhí)行結(jié)果如下：

y=

0.8200

0.7200

0.6300

0.6000

0.5500

0.5000

A=

1.00001.0000

1.00000.7408

1.00000.4493

1.00000.3329

1.00000.2019

1.00000.1003

c=

0.4760

0.3413

擬合曲線與實(shí)驗(yàn)數(shù)據(jù)如圖4-1所示。圖4-1擬合曲線與試驗(yàn)數(shù)據(jù)點(diǎn)問題2：解下面的線性方程組。

解：系數(shù)矩陣非奇異，變換如下：

解問題2的程序如下：

clc;clearall;

A=[23-5;403;017]

B=[2;6;9]

x=[A\B(B'/A')']程序執(zhí)行結(jié)果如下：

A=

23-5

403

017

B=

2

6

9

x=

0.73640.7364

1.87271.8727

1.01821.01824.1.3矩陣求逆

矩陣求逆可用于解線性方程組，被求逆的矩陣必須是方陣。需要說明的是，解方程組推薦使用前面介紹的矩陣除法，原因是矩陣除法速度快、精度高。

矩陣非奇異，inv(A)函數(shù)給出A逆陣的唯一解，奇異則給出無窮解或不定解。程序如下：

clc;clearall;

A=[23-5;403;017]

B=[2;6;9]

An=inv(A)%求A的逆陣An

x1=det([BA(:,2:3)])/det(A);%克萊姆法則求x1的解，det(A)求A的行列式值

x2=det([A(:,1)BA(:,3)])/det(A); %克萊姆法則求x2

x3=det([A(:,1:2)B])/det(A);%克萊姆法則求x3

x=[A\BAn*B[x1;x2;x3]]%左除、逆陣乘、克萊姆等三法的方程組解比較

程序執(zhí)行結(jié)果如下：

A=

23-5

403

017

B=

2

6

9

An=

?0.0273 ?0.2364 -0.0818

?0.2545 -0.1273 ?0.2364

-0.0364 ?0.0182 ?0.1091

x=

0.73640.73640.7364

1.87271.87271.8727

1.01821.01821.01824.1.4Doolittle分解(LU分解)

使用高斯消去法將任意矩陣A分解為下三角矩陣L與上三角矩陣U，這兩個矩陣的乘積就是原矩陣，稱作LU分解或Doolittle分解，即任意矩陣的三角分解。

表達(dá)式[L,U]=lu(A)給出矩陣A的LU分解，用A-L*U驗(yàn)證分解結(jié)果的正確性。程序如下：

clc;clearall;

A=[3-2-0.90;-2410;00-10;-0.5-0.50.11]

[L,U]=lu(A)

E=A-L*U%計(jì)算LU分解的誤差

程序執(zhí)行結(jié)果如下：

A=

3.0000 -2.0000 -0.9000 0

-2.0000 4.0000 1.0000 0

0 ?0 -1.0000 0

-0.5000 -0.5000 0.1000 ?1.0000

L=

1.00000?0 0

-0.6667 1.00000 0

00 ?1.0000 0

-0.1667 -0.3125-0.0750 ?1.0000

U=

3.0000-2.0000 -0.9000 0

0?2.6667 0.4000 0

00 -1.0000 0

000?1.0000

E=

1.0e-016*

0000

0000

0000

0-0.5551004.1.5Cholesky分解

使用高斯消去法將對稱正定矩陣X分解為一個上三角矩陣R和它的轉(zhuǎn)置矩陣(下三角矩陣)，R的轉(zhuǎn)置與R的乘積就是原矩陣，稱作Cholesky分解，即對稱正定矩陣的三角分解。

函數(shù)chol(X)給出對稱正定矩陣X的Cholesky分解，用X-R’*R驗(yàn)證分解結(jié)果的正確性。程序如下：

clc;clearall;

X=pascal(5)

R=chol(X)

E=X-R’*R %計(jì)算Cholesky分解的誤差

程序執(zhí)行結(jié)果如下：

X=

11111

12345

1361015

14102035

15153570

R=

11111

01234

00136

00014

00001

E=

00000

00000

00000

00000

000004.1.6QR分解

使用高斯消去法將任意矩陣X分解為一個一般矩陣Q和上三角矩陣R，Q與R的乘積就是原矩陣，稱作QR分解，即任意矩陣的正交三角分解。

表達(dá)式[Q,R]=qr(X)函數(shù)給出任意矩陣X的QR分解，用X-Q*R驗(yàn)證分解結(jié)果的正確性。程序如下：

clc;clearall;

X=[1324;5876;1291110]

[Q,R]=qr(X)

E=X-Q*R %計(jì)算QR分解的誤差程序執(zhí)行結(jié)果如下：

X=

1324

5876

1291110

Q=

-0.0767 -0.4737

-0.8774

-0.3835 -0.7982

0.4645

-0.9204

0.3721 -0.1204

R=

-13.0384 -11.5812 -12.9617 -11.8113

0 -4.4583 -2.4422 -2.9634

0 0

0.1720 -1.9267

E=

1.0e-014*

0.1110 0.17760.53290.4441

0

0.08880.17760.0888

0 0.35530.53290.35534.1.7特征值和特征向量

特征值和特征向量可用于求函數(shù)極值和搜尋主軸的坐標(biāo)系變換。方陣A的特征值矩陣λ和特征向量矩陣X滿足特征方程AX?=?Xλ，其中成立Axj=?λjxj，而λj為第j個特征值，xj為對應(yīng)λj的第j個特征向量且位于X矩陣的第j列。特征方程AX?=?Xλ展開如下：方陣A的廣義特征方程為AX?=?BXλ。

若A為實(shí)對稱陣，則λ為實(shí)數(shù)陣，否則λ為復(fù)數(shù)陣；如果A陣中有較小的元素，在計(jì)算特征值或者特征向量時增加選項(xiàng)“nobalance”以減少計(jì)算誤差。

表達(dá)式[X,lambda]=eig(A)給出矩陣A的特征陣lambda和特征向量陣X，表達(dá)式[X,lambda]=eig(A,b)給出矩陣A的廣義特征陣lambda和廣義特征向量陣X，程序如下：

clc;clearall;

A=[3-2-0.90;-2410;00-10;-0.5-0.50.11]

B=pascal(4)

[X1,lambda1]=eig(A,’nobalance’) %求矩陣A的特征陣和特征向量陣

[X2,lambda2]=eig(A,B) %求矩陣A的廣義特征陣和廣義特征向量陣

E1=A*X1-X1*lambda1 %計(jì)算特征方程誤差

E2=A*X2-B*X2*lambda2 %計(jì)算廣義特征方程誤差程序執(zhí)行結(jié)果如下：

A=

3.0000-2.0000-0.90000

-2.0000

4.0000

1.0000 0

00-1.0000 0

-0.5000-0.5000

0.1000

1.0000

B=

1111

1234

13610

141020

X1=

0.7808-0.4924 -0.0000 -0.1563

-1.0000-0.3845 -0.0000

0.1375

0?0 ?0 -1.0000

0.0240

1.0000 -1.0000

0.0453lambda1=

5.5616 0 0

0

0

1.4384 0

0

0 0

1.0000

0

0 0 0-1.0000

X2=

0.4963

0.0588?1.0000

0.0269

-1.0000

0.4716

0.9615-0.1360

0.7400-1.0000-0.4977

0.4337

-0.1947

0.4018-0.0435-1.0000lambda2=

67.9381 0 0 0

0 -1.9674 0 0

0 0 ?1.0736 0

0 0 0

0.0558

E1=

1.0e-014*

-0.2665 -0.0333 -0.0353-0.0167

0.4441 0.1221 0.0294 ?0.0028

?0 ?0 ?0 0

0.0333 -0.0222 0.0111 0

E2=

1.0e-013*

?0.2309 ?0.0269 0.0022 -0.0073

-0.0977 -0.0189 -0.0333 0.0056

-0.2076 -0.0044 -0.0078 ?0.0072

-0.5235 0.0236 0.0111 -0.00224.1.8奇異值分解

奇異值分解可應(yīng)用于最小二乘問題、廣義逆矩陣計(jì)算和其它許多技術(shù)領(lǐng)域。若A為實(shí)矩陣，則它可分解為正交陣P、奇異值對角陣λ、酉矩陣Q(滿足QTQ?=?I單位陣)轉(zhuǎn)置三因子的乘積，表示為下面的奇異值分解式為：表達(dá)式[P,lambda,Q]=svd(A)給出矩陣A奇異值分解式中的三個矩陣P、λ和Q。程序如下：

clc;clearall;

A=[12;34;56;78]

[P,lambda,Q]=svd(A) %求矩陣A奇異值分解式中的各個矩陣

E=A-P*lambda*Q‘ %計(jì)算奇異值分解式誤差

程序執(zhí)行結(jié)果如下：

A=

12

34

56

78P=

-0.1525 -0.8226 -0.3945-0.3800

-0.3499 -0.4214 0.2428 ?0.8007

-0.5474 -0.0201 0.6979 -0.4614

-0.7448 0.3812 -0.5462 ?0.0407

lambda=

14.2691?0

?00.6268

??0?0

?0?0

Q=

-0.6414?0.7672

-0.7672-0.6414E=

1.0e-014*

-0.0888-0.4885

0-0.1776

0-0.2665

-0.0888-0.5329考察下面的例子，一個最高冪次等于4的4階多項(xiàng)式，可轉(zhuǎn)換為如下的矩陣形式：4.2多?項(xiàng)?式?運(yùn)?算4.2.1多項(xiàng)式表達(dá)

在MATLAB中，多項(xiàng)式用降冪排列的系數(shù)行向量P表示，注意它包含零系數(shù)的項(xiàng)。

用方括弧[]生成多項(xiàng)式。程序如下：

clc;clearall;

P1=[1-12025116] %創(chuàng)建4階多項(xiàng)式

P2=[123] %創(chuàng)建2階多項(xiàng)式

程序執(zhí)行結(jié)果如下：

P1=

1-12025116

P2=

1234.2.2多項(xiàng)式算術(shù)運(yùn)算與微分

roots(P)函數(shù)解出多項(xiàng)式P的根，poly(R)函數(shù)解出以R為根的多項(xiàng)式，polyder(P)函數(shù)給出多項(xiàng)式P的微分，conv(P1,P2)函數(shù)給出多項(xiàng)式P1與P2的乘積結(jié)果，deconv(P1,P2)函數(shù)給出多項(xiàng)式P1除以P2的商和余數(shù)。多項(xiàng)式加減法與矩陣加減相同，但需補(bǔ)0以保證加減的兩矩陣維數(shù)相同。程序如下：

clc;clearall;

P1=[1-12025116]

P2=[123]

R1=roots(P1)

P10=poly(R1)

P3=P1+[zeros(1,2)P2]

P4=P1-[zeros(1,2)P2]

P5=conv(P1,P2)

P6=polyder(P1)

[P12,R12]=deconv(P1,P2)

程序執(zhí)行結(jié)果如下：

P1=

1-12025116

P2=

123

R1=

11.7473

2.7028

-1.2251+1.4672i

-1.2251-1.4672i

P10=

1.0000-12.0000-0.000025.0000116.0000

P3=

1-12127119P4=

1-12-123113

P5=

1-10-21-11166307348

P6=

4-36025

P12=

1-1425

R12=

00017414.2.3多項(xiàng)式值計(jì)算

用polyval(P,x)函數(shù)計(jì)算多項(xiàng)式P在自變量x上的值y=Px。程序如下：

clc;clearall;

P1=[1-12025116];

x=0:0.5:14;

y=polyval(P1,x);

plot(x,y,'--*');boxoff;

xlabel('x'),ylabel('y=Px')

程序執(zhí)行結(jié)果如圖4-2所示。圖4-2多項(xiàng)式值計(jì)算試驗(yàn)測得人工控制變量與其響應(yīng)變量的n組數(shù)據(jù)，若解析描述變量間關(guān)系或?yàn)楣こ虘?yīng)用提供數(shù)據(jù)，常常需要補(bǔ)充沒有試驗(yàn)測定過的變量值，此任務(wù)可通過“數(shù)據(jù)擬合”或“數(shù)據(jù)插值”的方法實(shí)現(xiàn)。

假定試驗(yàn)測定的數(shù)據(jù)是準(zhǔn)確的，在試驗(yàn)數(shù)據(jù)范圍之內(nèi)插入與試驗(yàn)數(shù)據(jù)規(guī)律近似的數(shù)值，稱作數(shù)據(jù)內(nèi)插法；在試驗(yàn)數(shù)據(jù)范圍之外插入與試驗(yàn)數(shù)據(jù)規(guī)律近似的數(shù)值，稱作數(shù)據(jù)外插(外推)法。數(shù)據(jù)的內(nèi)插或外插合稱為數(shù)據(jù)插值(datainterpolation)。一般不推薦采用外插。4.3數(shù)據(jù)插值如圖4-3所示，若試驗(yàn)測定變量(x,y)獲得容量為的一個樣本，則可按某種算法即插值函數(shù)計(jì)算插值點(diǎn)實(shí)現(xiàn)點(diǎn)插值。其中，插值函數(shù)穿過所有試驗(yàn)數(shù)據(jù)點(diǎn)，根據(jù)該函數(shù)類型可將插值方法分成若干類。圖4-3插值原理4.3.1插值方法

表達(dá)式Y(jié)=interp1(x,y,X,method,option)可實(shí)現(xiàn)一元函數(shù)插值。其中，輸入?yún)?shù)x、y分別為自變量和因變量的試驗(yàn)測定值，X為自變量的指定值，Y為因變量依照插值函數(shù)的計(jì)算值，參數(shù)method是插值方法選項(xiàng)字符串，參數(shù)option是外推選項(xiàng)字符串?？蛇x的插值方法有如下幾種：

(1)?method='nearest'：按最鄰近試驗(yàn)點(diǎn)原則在指定點(diǎn)(X,Y)上實(shí)施插值；

(2)?method='linear'：按線性函數(shù)在指定點(diǎn)(X,Y)上實(shí)施插值；

(3)?method='spline'：按分段樣條函數(shù)在指定點(diǎn)(X,Y)上實(shí)施插值；

(4)?method='pchip'：遵照形態(tài)保護(hù)原則按分段三次函數(shù)在指定點(diǎn)(X,Y)上實(shí)施插值；

(5)?method='cubic'：按分段三次函數(shù)在指定點(diǎn)(X,Y)上實(shí)施插值；

(6)?method='v5cubic'：按MATLAB5的三次函數(shù)規(guī)則在指定點(diǎn)(X,Y)上實(shí)施內(nèi)插，若插值點(diǎn)是外插則執(zhí)行樣條插值。

MATLAB通過參數(shù)option控制插值規(guī)則：

(1)?option缺省：只允許內(nèi)插；

(2)?option='extrap'：允許在自變量x試驗(yàn)區(qū)間外的插值點(diǎn)(X,Y)上實(shí)施外插；

(3)?option='pp'：允許在自變量x試驗(yàn)區(qū)間外的插值點(diǎn)(X,Y)上實(shí)施pp形式(分段多項(xiàng)式)外插。4.3.2一元函數(shù)插值

表達(dá)式Y(jié)=interp1(x,y,X,method)根據(jù)試驗(yàn)數(shù)據(jù)點(diǎn)(x,y)和插值方法method在自變量X的指定值處計(jì)算響應(yīng)變量Y的值，即依照插值函數(shù)Y?=?f?(X)插值(X,Y)。下面程序遍歷各種插值方法實(shí)施插值。

clc;closeall;clearall;

x=[00.631.261.892.513.143.774.405.035.656.28] %試驗(yàn)值x

y=[00.201.383.555.334.48-0.89-10.87-22.86-31.84-31.94] %試驗(yàn)值y

x2=0:pi/20:2*pi;y2=x2.^2.*sin(0.85*x2); %原函數(shù)計(jì)算值

X=pi/10:pi/5:2*pi; %指定自變量插值點(diǎn)

Y=interp1(x,y,X,'nearest'); %響應(yīng)變量鄰近點(diǎn)插值

%Y=interp1(x,y,X,'linear'); %響應(yīng)變量線性插值

%Y=interp1(x,y,X,'spline'); %響應(yīng)變量樣條插值%Y=interp1(x,y,X,'pchip');%響應(yīng)變量形態(tài)保護(hù)分段三次函數(shù)插值

%Y=interp1(x,y,X,'cubic');%響應(yīng)變量分段三次函數(shù)插值

%Y=interp1(x,y,X,'v5cubic');%響應(yīng)變量MATLAB5三次函數(shù)插值

xmin=min([min(x),min(X),min(x2)]); %橫軸最小值

xmax=max([max(x),max(X),max(x2)]); %橫軸最大值

ymin=min([min(y),min(Y),min(y2)]); %縱軸最小值

ymax=max([max(y),max(Y),max(y2)]); %縱軸最大值

plot(x,y,'ob',X,Y,'pk',x2,y2,'--r','LineWidth',2,'MarkerSize',10);boxoff;

axis([xminxmaxyminymax]);set(gca,'LineWidth',1,'FontSize',18,'FontName','Times');

xlabel('x','FontSize',18,'FontName','Times');

ylabel('f(x)','FontSize',18,'FontName','Times');

legend('試驗(yàn)數(shù)據(jù)點(diǎn)','插值點(diǎn)','原函數(shù)y=f(x)','FontSize',18,'FontName','Times')

圖4-4展示了試驗(yàn)數(shù)據(jù)點(diǎn)、插值點(diǎn)和原函數(shù)曲線(變量x,y的準(zhǔn)確值)。試比較插值效果。Nearest插值法Linear插值法Spline插值法Pchip插值法Cubic插值法V5cubic插值法圖4-4一元函數(shù)插值的插值效果對照4.3.3二元函數(shù)插值

表達(dá)式Y(jié)=interp2(x1,x2,y,X1,X2,method)根據(jù)試驗(yàn)數(shù)據(jù)點(diǎn)(x1,x2,y)和插值方法method在自變量(X1,X2)指定值處計(jì)算響應(yīng)變量Y的值，即依照插值函數(shù)Y?=?f(X1,X2)插值(X1,X2,Y)。程序如下：

clc;closeall;clearall;

[x1,x2]=meshgrid(-3:0.35:3); %創(chuàng)建自變量的網(wǎng)格點(diǎn)(x1,x2)

y=peaks(x1,x2); %用多峰函數(shù)計(jì)算網(wǎng)格點(diǎn)(x1,x2)上的響應(yīng)值y

mesh(x1,x2,y); %繪制多峰函數(shù)網(wǎng)格圖

[X1,X2]=meshgrid(-3:.125:3); %指定自變量的網(wǎng)格插值點(diǎn)(X1,X2)

Y=interp2(x1,x2,y,X1,X2,'nearest');%在網(wǎng)格插值點(diǎn)上用鄰近法插值Y

%Y=interp2(x1,x2,y,X1,X2,'linear');%在網(wǎng)格插值點(diǎn)上線性插值Y

%Y=interp2(x1,x2,y,X1,X2,'spline');%在網(wǎng)格插值點(diǎn)上用樣條法插值Yholdon,mesh(X1,X2,Y+15);holdoff;%繪制平移15的插值曲面網(wǎng)格圖

axis([-33-33-520]);

xlabel('x_1','FontSize',18,'FontName','Times');

ylabel('x_2','FontSize',18,'FontName','Times');

zlabel('y','FontSize',18,'FontName','Times');

set(gca,'LineWidth',1,'FontSize',18,'FontName','Times');

程序執(zhí)行結(jié)果如圖4-5所示。Nearest插值法Linear插值法Spline插值法Cubic插值法圖4-5二元函數(shù)插值的插值效果對照4.3.4三元函數(shù)插值

表達(dá)式Y(jié)=interp3(x1,x2,x3,y,X1,X2,X3,method)根據(jù)試驗(yàn)數(shù)據(jù)點(diǎn)(x1,x2,x3,y)和插值方法method在自變量(X1,X2,X3)指定值處插值響應(yīng)變量Y，即依照插值函數(shù)Y?=?f(X1,X2,X3)插值(X1,X2,X3,Y)。程序如下：

clc;closeall;clearall;

x1=linspace(0,2*pi,9);x2=x1;x3=x1; %自變量x1,x2,x3的試驗(yàn)值

[x1,x2,x3]=meshgrid(x1,x2,x3); %創(chuàng)建自變量的網(wǎng)格點(diǎn)(x1,x2,x3)

y=sin(x1).*cos(x2).*x3.^2; %自變量網(wǎng)格點(diǎn)(x1,x2,x3)上的響應(yīng)值y

X1=[pi/3pi/6];X2=[pi/3pi/4];X3=[pi/3pi/5]; %自變量插值點(diǎn)(X1,X2,X3)

[X1,X2,X3]=meshgrid(X1,X2,X3);%創(chuàng)建自變量插值網(wǎng)格點(diǎn)(X1,X2,X3)Y1=interp3(x1,x2,x3,y,X1,X2,X3,'nearest')%(X1,X2,X3)上的響應(yīng)變量插值Y1

Y2=interp3(x1,x2,x3,y,X1,X2,X3,'linear')%(X1,X2,X3)上的響應(yīng)變量插值Y2

Y3=interp3(x1,x2,x3,y,X1,X2,X3,'spline')%(X1,X2,X3)上的響應(yīng)變量插值Y3

Y4=interp3(x1,x2,x3,y,X1,X2,X3,'cubic')%(X1,X2,X3)上的響應(yīng)變量插值Y4

程序執(zhí)行結(jié)果如下：

Y1(:,:,1)=

0.30840.3084

0.30840.3084

Y1(:,:,2)=

0.30840.3084

0.30840.3084Y2(:,:,1)=

0.46800.2742

0.70200.4112

Y2(:,:,2)=

0.18720.1097

0.28080.1645

Y3(:,:,1)=

0.47090.2757

0.66920.3918

Y3(:,:,2)=

0.16950.0992

0.24090.1410Y4(:,:,1)=

0.47190.2849

0.66460.4012

Y4(:,:,2)=

0.16990.1026

0.23930.14444.3.5n元函數(shù)插值

n元函數(shù)插值依照二元函數(shù)和三元函數(shù)插值類推。其插值函數(shù)的調(diào)用格式如下：

Y=interpn(x1,x2,x3,...,y,X1,X2,X3,...)

Y=interpn(x1,x2,x3,...,y,X1,X2,X3,...,method)

Y=interpn(x1,x2,x3,...,y,X1,X2,X3,...,method,extrapval)

與二元函數(shù)插值和四元函數(shù)插值一樣，可選的插值方法method有如下4種：

(1)?method='nearest'：按最鄰近試驗(yàn)點(diǎn)原則在指定點(diǎn)(X,Y)上實(shí)施插值；

(2)?method='linear'：按線性函數(shù)在指定點(diǎn)(X,Y)上實(shí)施插值；

(3)?method='spline'：按分段樣條函數(shù)在指定點(diǎn)(X,Y)上實(shí)施插值；

(4)?method='cubic'：按分段三次函數(shù)在指定點(diǎn)(X,Y)上實(shí)施插值。4.3.6一元周期函數(shù)的FFT法插值

語句Y=interpft(y,N)采用FFT算法將一元周期函數(shù)y(長度為M的采樣序列)變換至傅里葉數(shù)域，在該數(shù)域確定N個插值點(diǎn)，然后將N個插值點(diǎn)傅里葉反變換到原數(shù)域?qū)崿F(xiàn)插值Y。

語句Y=interpft(y,N)根據(jù)原周期函數(shù)y執(zhí)行N點(diǎn)FFT法插值Y。程序如下：

clc;closeall;clearall;

K=6;M=15;N=K*M;%指定采樣倍率K、原函數(shù)采樣點(diǎn)數(shù)M和插值點(diǎn)數(shù)N

x=linspace(0,2*pi,M);y=sin(x);%原周期函數(shù)y=f(x)的采樣點(diǎn)(x,y)

R=max(x)-min(x); %原周期函數(shù)自變量x的極差R(x的數(shù)據(jù)范圍)

dx=R/M; %原周期函數(shù)自變量x的采樣間隔dx

dX=dx*M/N; %FFT法插值自變量X的采樣間隔dXX=x;

fori=1:K-1;X=[Xx+i*dX];end; %在x數(shù)據(jù)范圍內(nèi)生成自變量插值點(diǎn)X

X=sort(X);Y=interpft(y,N); %在自變量插值點(diǎn)X上按FFT法插值Y

X=X(1:N-K+1);Y=Y(1:N-K+1); %剔除x數(shù)據(jù)范圍外的插值點(diǎn)

plot(x,y,'ok',X,Y,'*b','MarkerSize',9,'LineWidth',1);boxoff; %繪圖比較

xmin=min([min(x)min(X)]);xmax=max([max(x)max(X)]);

ymin=min([min(y)min(Y)]);ymax=max([max(y)max(Y)]);

axis([xminxmaxyminymax]);

xlabel('x','FontSize',18,'FontName','Times');

ylabel('y','FontSize',18,'FontName','Times');

set(gca,'LineWidth',1,'FontSize',18,'FontName','Times');

legend('Originaldata','Interpolateddata');

程序執(zhí)行結(jié)果如圖4-6所示。圖4-6對照FFT法的插值點(diǎn)“*”與原函數(shù)點(diǎn)“”若響應(yīng)變量y與自變量系之間的關(guān)系可由下面的模型表述：4.4數(shù)據(jù)擬合4.4.1多項(xiàng)式最小二乘擬合

下面以一個節(jié)流元件的流量特性為例討論數(shù)據(jù)擬合問題。節(jié)流元件的壓差流量模型如下：

人工調(diào)節(jié)壓差x的值，測定流量的響應(yīng)值y，獲得表4-2所示的樣本。試用2階多項(xiàng)式做數(shù)據(jù)擬合，即擬合模型為

，它是壓差流量模型的一個近似。表4-2節(jié)流元件流量與壓差的試驗(yàn)測定結(jié)果

polyfit(x,y,p)函數(shù)給出數(shù)據(jù)(x,y)的p階多項(xiàng)式擬合函數(shù)。程序如下：

clc;clearall;

x=0:0.1:1;y=[0.451.983.286.167.087.347.669.569.489.3011.20];

P=polyfit(x,y,2)

x1=0:0.01:1.1;y1=polyval(P,x1);

h=plot(x,y,'*k',x1,y1,'--r');boxoff;

axis([0max(x1)0max(max(y),max(y1))]);

set(gca,'FontSize',16);

set(h,'LineWidth',2,'MarkerSize',10);

xlabel('x');ylabel('y');

legend('試驗(yàn)數(shù)據(jù)點(diǎn)','y=-8.2413x^2+18.1513x+0.4897');

%title('擬合模型y=-9.8147x^2+20.1338x-0.0327');程序執(zhí)行結(jié)果如下：

P=

-8.241318.15130.4897

多項(xiàng)式擬合函數(shù)y?=?-8.2413x2?+?18.1513x?+?0.4897對照試驗(yàn)數(shù)據(jù)點(diǎn)，如圖4-7所示。圖4-7多項(xiàng)式擬合函數(shù)對照試驗(yàn)數(shù)據(jù)點(diǎn)4.4.2一元函數(shù)最小二乘擬合

下面以節(jié)流元件的壓差流量模型為擬合模型，做數(shù)據(jù)的最小二乘擬合。

k=lsqcurvefit(fun,k0,x,y)函數(shù)擬合數(shù)據(jù)(x,y)，并給出模型參數(shù)k的最小二乘解。程序如下：

clc;clearall;

x=0:0.1:1;y=[0.451.983.286.167.087.347.669.569.489.3011.20];

fun=@(k,x)k*sqrt(x);

k0=1;[k,Resnorm,Residual,Exitflag]=lsqcurvefit(fun,k0,x,y)

x1=0:0.01:1.1;y1=k*sqrt(x1);

h=plot(x,y,'*k',x1,y1,'--r');boxoff;

axis([0max(x1)0max(max(y),max(y1))]);

set(gca,'FontSize',16);

set(h,'LineWidth',2,'MarkerSize',10);

xlabel('x');ylabel('y');

legend('試驗(yàn)數(shù)據(jù)點(diǎn)','擬合函數(shù)y=10.4670x^1^/^2')程序執(zhí)行結(jié)果如下：

k=

10.4670

Resnorm=

6.1245

Residual=

-0.45001.33001.4010-0.4270-0.46010.06130.4477-0.8026

-0.11800.6299-0.7330

Exitflag=

1

最小二乘擬合函數(shù)

對照試驗(yàn)數(shù)據(jù)點(diǎn)，如圖4-8所示。圖4-8最小二乘擬合函數(shù)對照試驗(yàn)數(shù)據(jù)點(diǎn)數(shù)據(jù)分析包括簡單統(tǒng)計(jì)和其它數(shù)據(jù)描述工作，執(zhí)行數(shù)據(jù)分析的MATLAB函數(shù)量極多。

表4-3給出了一些常用數(shù)據(jù)分析函數(shù)，使用help命令可查閱詳細(xì)的函數(shù)調(diào)用格式。4.5數(shù)據(jù)分析表4-3MATLAB常用數(shù)據(jù)分析函數(shù)4.5.1數(shù)據(jù)概括

下面通過一個實(shí)例介紹max函數(shù)、min函數(shù)、mean函數(shù)、median函數(shù)和std函數(shù)的用法。程序如下：

clc;closeall;clearall;

x=[48271;35107;16715];y=[1357;2468;11233127];

x_max1=[max(x);max(x,[],1)] %求列向量最大值

x_max2=max(x,[],2) %求行向量最大值

x_max3=max(x,y) %求兩矩陣對應(yīng)元素最大值

x_min1=[min(x);min(x,[],1)] %求列向量最小值

x_min2=min(x,[],2) %求行向量最小值

x_min3=min(x,y) %求兩矩陣對應(yīng)元素最小值

x_mean1=[mean(x);mean(x,1)] %求列向量均值

x_mean2=mean(x,2) 求行向量均值

x_median1=[median(x);median(x,1)] %求列向量中值

x_median2=median(x,2) %求行向量中值

x_std1=[std(x);std(x,1)]

%求列向量標(biāo)準(zhǔn)差，除數(shù)n-1

x_std2=[std(x,0,1);std(x,1,1)]

%求列向量標(biāo)準(zhǔn)差，0除數(shù)n-1，1除數(shù)n

x_std3=[std(x,0,2)std(x,1,2)]

%求行向量標(biāo)準(zhǔn)差，0除數(shù)n-1，1除數(shù)n

x_coef=100*std(x)./mean(x)

%求列向量變異系數(shù)

程序執(zhí)行結(jié)果如下：

x_max1=

481071

481071

x_max2=

71

10

15

x_max3=

48571

35108

11233127x_min1=

1527

1527

x_min2=

2

3

1

x_min3=

1327

2467

16715

x_mean1=

2.66676.33336.333331.0000

2.66676.33336.333331.0000x_mean2=

21.2500

6.2500

7.2500

x_median1=

36715

36715

x_median2=

6.0000

6.0000

6.5000

x_std1=

1.52751.52754.041534.8712

1.24721.24723.299828.4722x_std2=

1.52751.52754.041534.8712

1.24721.24723.299828.4722

x_std3=

33.260328.8043

2.98612.5860

5.79515.0187

x_coef=

57.282224.118863.8124112.48774.5.2協(xié)差陣和相關(guān)陣

協(xié)差陣函數(shù)cov、相關(guān)陣函數(shù)corrcorf和方差函數(shù)var。程序如下：

x=[482716;3510723;1671545]

x_cov=cov(x); %求矩陣x中列向量的協(xié)差陣

x_corr=corrcoef(x) %求矩陣x中列向量的相關(guān)陣

x_var=var(x) %求矩陣x中列向量的方差

程序執(zhí)行結(jié)果如下：

x=

482716

3510723

1671545x_cov=1.0e+003*

0.00230.0012-0.00280.0360-0.0297

0.00120.0023-0.00620.0520-0.0178

?-0.0028-0.00620.0163-0.13600.0442

0.03600.0520-0.13601.2160-0.5160

-0.0297-0.01780.0442-0.51600.3823

x_corr=

1.00000.5000-0.45900.6758-0.9933

0.50001.0000-0.99890.9762-0.5971

?-0.4590-0.99891.0000-0.96500.5589

0.67580.9762-0.96501.0000-0.7568

?-0.9933-0.59710.5589-0.75681.0000

x_var=1.0e+003*

0.00230.00230.01631.21600.38234.5.3矩陣元素的求和與累加和

求和函數(shù)sum與累加和函數(shù)cumsum。程序如下：

x=[482716;3510723;1671545]

sx1=[sum(x);sum(x,1)] %求矩陣x中列向量元素的和

sx2=sum(x,2) %求矩陣x中行向量元素的和

csx1=[cumsum(x)cumsum(x,1)] %求矩陣x中列向量元素的累加和

csx2=cumsum(x,2) %求矩陣x中行向量元素的累加和

程序執(zhí)行結(jié)果如下：

x=

482716

3510723

1671545sx1=

819199374

819199374

sx2=

91

48

74

csx1=

482716482716

713127829713127829

819199374819199374

csx2=

412148591

38182548

171429744.5.4矩陣元素的乘積與累乘積

乘積函數(shù)prod與累乘積函數(shù)cumprod。程序如下：

x=[482716;3510723;1671545]

px1=[prod(x);prod(x,1)] %求矩陣x中列向量元素的乘積

px2=prod(x,2) %求矩陣x中行向量元素的乘積

cpx1=[cumprod(x);cumprod(x,1)] %求矩陣x中列向量元素的累乘積

cpx2=cumprod(x,2) %求矩陣x中行向量元素的累乘積

程序執(zhí)行結(jié)果如下：

x=

482716

3510723

1671545px1=

1224014074556210

1224014074556210

px2=

27264

24150

28350

cpx1=

482716

124020497138

1224014074556210

482716

124020497138

1224014074556210cpx2=

43264454427264

315150105024150

1642630283504.5.5數(shù)據(jù)排序

元素排序函數(shù)sort和向量排序函數(shù)sortrows。程序如下：

x=[482716;3510723;1671545]

y1=sort(x,2) %矩陣x中每行元素按升序方式排序

y2=[sort(x);sort(x,1)] %矩陣x中每列元素按升序方式排序

y3=sort(x,2,'descend') %矩陣x中每行元素按降序方式排序

y4=sortrows(x,1)%矩陣x中每一行向量按第1列元素的升序方式排序

y5=sortrows(x,3)%矩陣x中每一行向量按第3列元素的升序方式排序

程序執(zhí)行結(jié)果如下：

x=

482716

3510723

1671545

y1=

246871

3571023

1671545

y2=

15276

3671523

48107145

15276

3671523

48107145y3=

718642

2310753

4515761

y4=

1671545

3510723

482716

y5=

482716

1671545

35107234.6.1數(shù)值積分

trapz(x,y)函數(shù)計(jì)算梯形面積之和逼近函數(shù)積分。程序如下：

x=0:pi/10:pi;

y=sin(x);

Z1=trapz(x,y) %計(jì)算f(X)dx之和，dx=pi/10

Z2=pi/10*trapz(y) %計(jì)算f(x)×1之和與pi/10的乘積

程序執(zhí)行結(jié)果如下：

Z1=

1.9835

Z2=

1.98354.6數(shù)值積分與數(shù)值微分

quad(fx,a,b,tol)函數(shù)在x區(qū)間［a,b］上對函數(shù)fx采樣并用Simpson遞歸法(ecursiveadaptiveSimpsonquadrature)以允差tol計(jì)算數(shù)值積分逼近函數(shù)積分。程序如下：

fx=@(x)sin(x);

Q=quad(fx,0,pi,1e-2)

程序輸出如下：

Q=

2.0000

trapz(x,y)函數(shù)數(shù)值積分原理，如圖4-9所示。圖4-9trapz(x,y)函數(shù)數(shù)值積分原理

quadgk(fx,a,b,tol)函數(shù)在x區(qū)間［a,b］上對函數(shù)fx采樣并用Gauss-Kronrod遞歸法(adaptiveGauss-Kronrodquadrature)以允差tol計(jì)算數(shù)值積分逼近函數(shù)積分。程序如下：

fx=@(x)sin(x);

Q=quadgk(fx,0,pi,'AbsTol',1e-2)

程序輸出如下：

Q=

2

quadl(fun,a,b,tol)函數(shù)在x區(qū)間［a,b］上對函數(shù)fx采樣并用Lobatto遞歸法(recursiveadaptiveLobattoquadrature)以允差tol計(jì)算數(shù)值積分逼近函數(shù)積分。程序如下：

fx=@(x)sin(x);

Q=quadl(fx,0,pi,1e-2)

程序輸出如下：

Q=

2.0000

quadv(fx,a,b,tol)函數(shù)在x區(qū)間[a,b]上對向量函數(shù)fx(多個函數(shù)組成的向量)采樣，采用矢量化方法(VectorizedAdaptiveQuadrature)以允差tol計(jì)算每個函數(shù)的數(shù)值積分。程序如下：

fx=@(x)sin(x+(0:pi/2:pi));

Qv=quadv(fx,0,pi,1e-2)

程序輸出如下：

Qv=

2.00000.0000-2.0000

quadv函數(shù)對所示三條函數(shù)曲線積分如圖4-10所示。圖4-10quadv函數(shù)對所示三條函數(shù)曲線積分

quad2d(fxy,a,b,c,d,param1,val1,param2,val2,...)函數(shù)在x區(qū)間[a,b]及y區(qū)間[c,d]組成的區(qū)域上對二元函數(shù)f(x,y)采樣，采用矢量化方法(VectorizedAdaptiveQuadrature)計(jì)算二重?cái)?shù)值積分以逼近函數(shù)積分(注意，MATLAB2009a以前的版本無此函數(shù))。程序如下：

fxy=@(x,y)abs(x.^2+y.^2-0.25);

c=@(x)-sqrt(1-x.^2);

d=@(x)sqrt(1-x.^2);

Q2d=quad2d(fun,-1,1,c,d,'AbsTol',1e-5)

程序輸出如下：

Q2d=

0.98174.6.2數(shù)值微分

積分是函數(shù)乘以自變量微分的累積，因此積分結(jié)果對函數(shù)變化不敏感。導(dǎo)數(shù)是函數(shù)微分與自變量微分之比，較小函數(shù)變化有時也會產(chǎn)生巨大的導(dǎo)數(shù)，故微分問題對函數(shù)變化極為敏感。因此，數(shù)值微分是一種可能導(dǎo)致巨大誤差的計(jì)算工作，應(yīng)當(dāng)盡量避免數(shù)值微分，尤其是對試驗(yàn)數(shù)據(jù)直接數(shù)值微分。較合理的做法是，先對試驗(yàn)數(shù)據(jù)進(jìn)行數(shù)據(jù)擬合，再對擬合數(shù)據(jù)執(zhí)行數(shù)值微分。

函數(shù)diff(y)對y作1階差分，diff(y,k)對y作k階差分。程序如下：clc;closeall;clearall;

x=0:pi/10:2*pi;n=length(x);y=exp(-x).*sin(2*x); %函數(shù)y=f(x)和微分點(diǎn)(x,y)

dx=diff(x);dy=diff(y); %計(jì)算x差分和y差分

X=x;P=polyfit(x,y,7);Y=polyval(P,X); %采用7階多項(xiàng)式在X上擬合Y值

dX=diff(X);dY=diff(Y); %計(jì)算X差分和Y差分

續(xù)程序(可省略)：

plot(x,y,'--r',x(1:n-1),dy./dx,'ob',X(1:n-1),dY./dX,'pk');%繪圖比較數(shù)值導(dǎo)數(shù)

ymin=min(dy./dx);ymax=max(dy./dx);

Ymin=min(dY./dX);Ymax=max(dY./dX);

axis([min(x)max(x)min(ymin,Ymin)max(ymax,Ymax)]);boxoff;

set(gca,'FontSize',18,'FontName','Times')

xlabel('x','FontSize',18,'FontName','Times');

ylabel('dy/dx','FontSize',18,'FontName','Times');

legend('被微分函數(shù)y=e^-^x*sin(2x)','直接法數(shù)值導(dǎo)數(shù)dy/dx','擬合法數(shù)值導(dǎo)數(shù)dY/dX');

程序輸出如圖4-11所示。圖4-11直接法與擬合法的數(shù)值微分比較語句[Fx1,Fx2,...]=gradient(F,h1,h2,...)以x1間隔h1、x2間隔h2等計(jì)算函數(shù)F的數(shù)值梯度，語句[Fx1,Fx2,...]=gradient(F)缺省h1=h2=…=1。若函數(shù)F是向量，則為一元函數(shù)，若函數(shù)F是n維矩陣，則為多元函數(shù)，矩陣每一個元素對應(yīng)n個自變量的一組值，維序號與自變量序號一致。程序如下：

clc;closeall;clearall;

v=0:pi/10:2*pi;

[x1,x2]=meshgrid(v);

F=sin(0.85*x1).*cos(0.75*x2);

[Fx1,Fx2]=gradient(F,pi/5,pi/5);

figure,mesh(x1,x2,F);view(-25,45);boxoff;

續(xù)程序(可省略)：

set(gca,'FontSize',18,'FontName','Times');

xlabel('x_1','FontSize',18,'FontName','Times');

ylabel('x_2','FontSize',18,'FontName','Times');

zlabel('F(x_1,x_2)','FontSize',18,'FontName','Times');

figure,contour(v,v,F)

holdon,quiver(v,v,Fx1,Fx2),holdoff

set(gca,'FontSize',18,'FontName','Times');

xlabel('x_1','FontSize',18,'FontName','Times');

ylabel('x_2','FontSize',18,'FontName','Times');

程序輸出如圖4-12所示。圖4-12二元函數(shù)的梯度計(jì)算和結(jié)果展示(a)二元函數(shù)

(b)函數(shù)y的等高線與梯度

MATLAB可解常微分方程(ODE：ordinarydifferentialequations)組，其標(biāo)準(zhǔn)型如下：4.7解常微分方程和常微分方程組或

將常微分方程組編寫成一個函數(shù)，需符合上面所示的標(biāo)準(zhǔn)型輸入輸出，其中：4.7.1解常微分方程

求解下面的一維常微分方程問題。

解：式(4.7.1