matlab基礎(chǔ)及應(yīng)用(教程全書)part2

第4章程序設(shè)計(jì)

在前面我們已經(jīng)看到，MATLAB不但可以在命令窗直接輸入命令并運(yùn)行，而且還可以生

成自己的程序文件，這就是我們通常說(shuō)的一類以M為后綴的M文件，本章我們就來(lái)研究這類

文件的形成方法。

M文件可分分為兩大類，一是命令式M文件（也稱為腳本文件，script）,二是函數(shù)式M

文件（function）。兩類文件的區(qū)別在于：

（1）命令式文件可以直接運(yùn)行，函數(shù)式文件不能直接運(yùn)行，只能調(diào)用。

（2）命令式文件運(yùn)行時(shí)沒(méi)有輸入輸出參量，函數(shù)式文件在調(diào)用時(shí)需要進(jìn)行輸入輸出

參量設(shè)置。

（3）命令式文件運(yùn)行中可以調(diào)用工作空間的數(shù)據(jù)，運(yùn)行中產(chǎn)生的所有變量為全局變

量。

（4）函數(shù)式文件不能調(diào)用工作空間的數(shù)據(jù)，運(yùn)行中產(chǎn)生的所有變量為局部變量。命

令式文件運(yùn)行中產(chǎn)生的所有變量為全局變量，可以調(diào)用和存儲(chǔ)到工作空間的數(shù)據(jù)。

4.1MATLAB的程序文件-M文件

4.1.1腳本文件（Scripts）

當(dāng)我們需要在命令窗進(jìn)行大量的命令集合運(yùn)行時(shí),直接從命令窗口輸入比較麻煩，這時(shí)

就可以將這些命令集合存放在一個(gè)腳本文件（Scripts）中，運(yùn)行時(shí)只需要輸入其文件名就

可以自動(dòng)執(zhí)行這些命令集合。需要注意的是，腳本文件運(yùn)行所產(chǎn)生的變量都駐留在MATLAB

的工作空間中，同時(shí)腳本文件也可以調(diào)用工作空間中的數(shù)據(jù)。因此，腳本文件所涉及的變量

是全局變量。前兒章所涉及到的M文件都是這類腳本文件。

編輯一個(gè)腳本文件可以直接在命令窗口的左上角打開(kāi)編輯窗進(jìn)行編輯。

4.1.2函數(shù)文件（function）

函數(shù)式文件（function）的構(gòu)成

（1）函數(shù)定義行：

Function［輸出參量］=gauss（輸入?yún)⒘浚?/p>

（2）：

完成函數(shù)的功能。

（3）函數(shù)說(shuō)明。

（4）函數(shù)行注。

從上面構(gòu)成的情況看，函數(shù)式文件實(shí)際上是完成輸入?yún)⒘颗c輸出參量的轉(zhuǎn)換，這樣的轉(zhuǎn)換是

由函數(shù)文件名為gauss的文件來(lái)完成的。函數(shù)體的功能必須說(shuō)明清楚輸入?yún)⒘颗c輸出參量的

關(guān)系。函數(shù)說(shuō)明是用來(lái)解釋該函數(shù)的功能的，函數(shù)行注是對(duì)程序行進(jìn)行說(shuō)明的。上面（1）

和（2）是必須的。

【例47】分析下面函數(shù)文件。

%一個(gè)數(shù)列，任意項(xiàng)等于前兩項(xiàng)之和，輸入項(xiàng)數(shù)可以給出這個(gè)數(shù)列

function[a]=sul(n)

ifn==1

a=1;

elseifn==2

a=2;

end

"1)=1;

b(2)=2;

fori=3:n

b(i)=b(i-2)+b(i-1);

end

a=b;

end

該函數(shù)文件的文件名為sul.m,在第一行給出了該函數(shù)的功能，即輸入項(xiàng)數(shù)就可以自動(dòng)給出

一個(gè)滿足條件的數(shù)列。定義行規(guī)定了輸入?yún)?shù)是該數(shù)列的項(xiàng)數(shù)，輸出參數(shù)是該數(shù)列。從第三

行起，是該函數(shù)的主體，主要說(shuō)明了輸入?yún)?shù)與輸出參數(shù)的關(guān)系。當(dāng)我們?cè)诿畲盎蚰_本文

件中調(diào)用該函數(shù)就會(huì)有結(jié)果，如

?p=sul(5)

P=

12358

【例4-2】分析下面函數(shù)文件。

%一個(gè)數(shù)列的通項(xiàng)公式為a(n+1尸a(n)+M2,給定任意項(xiàng)的值，求這個(gè)數(shù)列的后10項(xiàng)，并畫

離散函數(shù)圖

function[b]=shulie(n,zhi)%zhi為初值

b(1)=zhi

form=1:n-1

b(1+m)=b(m)+mA2

end

forj=n:n+9

b(1+j戶b(j)+『2

end

x=n:(n+9)

stem(x,b(n:n+9))

end

如果在命令窗輸入

shulie(5,10)

結(jié)果為

b=

10

x=

567891011121314

ans=

Columns1through5

1011152440

Columns6through10

65101150214295

Columns11through15

3955166608291025

以及

4.2程序的流程控制

4.2.1循環(huán)控制

1.硬循環(huán)語(yǔ)句(for-end)

所謂硬循環(huán)是指無(wú)條件的循環(huán)，其結(jié)構(gòu)為

for(循環(huán)變量)

end

【例4-3］用硬循環(huán)語(yǔ)句形成一個(gè)5階對(duì)角線元素為2的矩陣。

clear

n=5;%每個(gè)for和if都必須對(duì)應(yīng)一個(gè)end

fori=l:n%i是循環(huán)變量

forj=l:n

ifi==j

a(i,j)=2;

else

a(i,j)-0;

end

end

end

結(jié)果為

a=

20000

02000

00200

00020

00002

根據(jù)結(jié)果，如果要兩條主對(duì)角線元素都為2,如何調(diào)此程序?

注意，硬循環(huán)語(yǔ)句也可用break來(lái)跳出循環(huán)。

4.2.2條件控制

1.條件分支語(yǔ)句(if-else-end)如果-那么-否則

結(jié)構(gòu)

If條件1

語(yǔ)句1

elseif條件2

語(yǔ)句2

else

語(yǔ)句

end

注意這一結(jié)構(gòu)的條件優(yōu)先問(wèn)題，當(dāng)有多個(gè)條件時(shí)，如果條件1不滿足，再判斷elseif

后面的條件2…，如果所有的條件都不滿足，在則執(zhí)行else后面的語(yǔ)句。

【例4-5】分析下面條件語(yǔ)句。

x=input('x=?')

if((x+3)<2)

y=x*2;

elseif(x<2)

y二x'2;

else(x<l)

y=sqrt(x);

end

y

運(yùn)行結(jié)果

x=?l.5

x=

1.5000

y=

2.2500

2.條件循環(huán)語(yǔ)句(while-end)

當(dāng)條件滿足時(shí)執(zhí)行下一語(yǔ)句，否則就到end返回while。其結(jié)構(gòu)為

while(表達(dá)式)

OOO

end

【例4-4】求階乘大于或等于99-99的最小整數(shù)。

clear

n=l;

whileprod(l:n)<99^99;%prod()向量元素的乘積

n=n+l；

end

n

結(jié)果為

n=

120

4.2.3開(kāi)關(guān)與試探控制

分支開(kāi)關(guān)語(yǔ)句(switch-case-otherwise-end)

結(jié)構(gòu)switch(開(kāi)關(guān)量)，相當(dāng)于滿足什么條件做什么事。

【例4-5］根據(jù)開(kāi)關(guān)量的情況來(lái)決定計(jì)算機(jī)計(jì)算機(jī)執(zhí)行那個(gè)語(yǔ)句。

a=input('a=?f)

switcha

case1

disp('Itisraning1)

case0

disp('Itdonotkonw,)

case-1

disp(,Itisnotraning，)

otherwise

disp('Itisraning?，)

end

4.3程序設(shè)計(jì)

4.3.1程序的調(diào)試

1.程序的錯(cuò)誤一般分為語(yǔ)法和邏輯兩種錯(cuò)型。語(yǔ)法錯(cuò)誤時(shí)程序無(wú)法運(yùn)行，并在命令窗給予

提示。而邏輯錯(cuò)誤只能根據(jù)計(jì)算的具體情況來(lái)判斷。

2.根據(jù)出錯(cuò)信息調(diào)試

根據(jù)命令窗的提示，注意一般情況不加號(hào)調(diào)試。

3.設(shè)置斷點(diǎn)來(lái)判斷

dbstop格式：dbstopinM文件名at行號(hào)(也可用鼠標(biāo))

4.利用keyboard命令在文件中來(lái)判斷

當(dāng)出現(xiàn)k?時(shí)returrio

5.變量的鼠標(biāo)觀測(cè)法

當(dāng)程序代碼運(yùn)行后，每個(gè)變量的數(shù)據(jù)都可以用鼠標(biāo)放在相應(yīng)的位置處去觀察該變量的數(shù)據(jù)情

況，以判斷程序運(yùn)行的具體執(zhí)行情況。

6.代碼運(yùn)行的計(jì)時(shí)方法

(1)整段程序代碼的計(jì)時(shí)

tic…toe表示計(jì)算tic與toe之間的時(shí)間

【例4-6】計(jì)算機(jī)程序運(yùn)行時(shí)間計(jì)算。

tic

a=rand(300);

inv(a);%計(jì)算逆矩陣

toe

%注：程序中不需要顯示結(jié)果的就不顯示，這可以節(jié)約機(jī)時(shí)。

結(jié)果為

?

elapsedtime=

0.3900

(2)也可以用來(lái)計(jì)算tl,t2之間的時(shí)間差來(lái)完成上述功能。

【例4-7】計(jì)算機(jī)程序運(yùn)行時(shí)間計(jì)算。

t0=clock

a=rand(300);

inv(a);%計(jì)算逆矩陣

elapsedtime=etime(clock,tO)

%注：程序中不需要顯示結(jié)果的就不顯示，這可以節(jié)約機(jī)時(shí)。

結(jié)果仍為

?

elapsedtime=

0.3900

(3)也可以用eputime變量來(lái)完成上述功能。

【例4-8】計(jì)算機(jī)程序運(yùn)行時(shí)間計(jì)算。

t0=eputime

a=rand(300);

inv(a);%計(jì)算逆矩陣

e1apsed_time=cputime-t0

%注：程序中不需要顯示結(jié)果的就不顯示，這可以節(jié)約機(jī)時(shí)。

結(jié)果為

to=

9.5787e+003

elapsed_time=

0.2710

4.3.2M文件的性能優(yōu)化

1.程序代碼的向量化

【例4-8］比較下面兩個(gè)程序段，觀察其運(yùn)行的時(shí)間。

程序1

t0二eputime

n=100000;

total=0;

fori=l:n

total=total+l/i

end

total

tl=cputime-tO

運(yùn)行結(jié)果為

total=

12.0901

tl=

3.1950

程序2

tic

n=100000;

a=l:n;

total二sum(1./a)

toe

運(yùn)行結(jié)果為

total=

12.0901

elapsedtime=

0.0200

發(fā)現(xiàn)，程序1和程序2運(yùn)行結(jié)果完全相同，但運(yùn)行時(shí)間程序1是程序2的160倍！

2.用矩陣結(jié)構(gòu)進(jìn)行運(yùn)算

一般情況下，完全采用矩陣運(yùn)行的方式，MATLAB的程序與C語(yǔ)言的運(yùn)行基本相同。這

必須對(duì)矩陣非常熟練，例如

x=[l23;121]

a=[456]

如果希望將a中的每一個(gè)元素乘以x的每一列，用diag(x)。

3.陣的預(yù)先配置

在一些大的程序運(yùn)行時(shí)，為了變量不至于讓計(jì)算機(jī)不斷向內(nèi)存要空間位置而耽誤時(shí)間，

運(yùn)行前需要給變量留出?定的空間。

【例4-9]比較下面兩個(gè)程序段，觀察其運(yùn)行的時(shí)間。

程序1

tic

a=[l23;456;789];

fori=l:100

y(i)=det(a"i)

end

toe

結(jié)果為

elapsed_time=

0.1100

程序2

tic

a=[l23;456;789];

y=zeros(l,100);%分配內(nèi)存

fori=l:100

y(i)=det(a*i);

end

toe

結(jié)果為

elapsedtime=

0.0100

練習(xí)4

i.已知函數(shù)

y=/一ax+1

請(qǐng)研究參數(shù)a在11:0.1:1］內(nèi)函數(shù)的根(實(shí)部與虛部)隨a的變化關(guān)系。(用圖形表示)

2.編制?個(gè)函數(shù)，函數(shù)的輸入?yún)?shù)為一個(gè)任意矩陣或向量，輸出參數(shù)為該矩陣中不相同的

元素個(gè)數(shù)。

3.編制?個(gè)函數(shù)文件，函數(shù)的輸入?yún)?shù)為一個(gè)任意矢性場(chǎng)

(p=f(x,y)ex+g(x,y)ey

和這個(gè)矢性場(chǎng)自變量的取值范圍，要求這個(gè)函數(shù)能自動(dòng)畫出矢性場(chǎng)的平面場(chǎng)分布。

4.設(shè)計(jì)一個(gè)函數(shù)文件，要求輸入一個(gè)任意方陣，該函數(shù)能統(tǒng)計(jì)出方陣中大于零、等于零和

小于零的元素個(gè)數(shù)。

5.設(shè)計(jì)一個(gè)能實(shí)現(xiàn)下面功能的函數(shù)文件：輸入兩個(gè)參數(shù)：(1)任意2Xn矩陣(n>8,且第一

行為自變量，第二行為對(duì)應(yīng)的函數(shù)值)；(2)介于任意兩自變量之間的值。輸出為該自變量

對(duì)應(yīng)的函數(shù)值。(用9次多項(xiàng)式擬和)

6.一個(gè)MATLAB的遞歸函數(shù)fibo.m來(lái)計(jì)算fibo數(shù)列，定義如下：

fibo(n+2)=(fibo(n+1)+fibo(n))/2

此數(shù)列的初條件為

fibo(1)=0,fibo(2)=1

n的最大數(shù)為100,要求：

(1)保存你的fibo.m文件，當(dāng)在命令窗調(diào)用fibo函數(shù)時(shí)，不論輸入任何整數(shù)有正確的

輸出。

(2)做出fibo的二維離散函數(shù)圖，n取1到10,圖的函數(shù)值處用小圓圈并涂為黑色，請(qǐng)

保存你的圖形。

(3)用三次樣條插值的方法對(duì)(2)中的10個(gè)點(diǎn)進(jìn)行插值，自變量的分辨率為0.01,

請(qǐng)保存你的圖形。

(4)將完成(3)工作的插值函數(shù)保存為fib.m文件，，當(dāng)在命令窗調(diào)用fib函數(shù)時(shí)，不論輸

入任何具有兩位小數(shù)且小于10大于0的數(shù)(如5.45)時(shí)有正確的輸出。

7.設(shè)電子粒子束流從恒定磁場(chǎng)中某點(diǎn)以相同速率發(fā)射，發(fā)射的方向與磁場(chǎng)方向的夾角很小,

觀察不同方向入射的粒子束流的運(yùn)動(dòng)軌道。(設(shè)磁場(chǎng)沿Z方向)

數(shù)學(xué)模型：

粒子流的速度初值為

匕=%

匕=%6sin(a)

vv=%6cos(a)

標(biāo)準(zhǔn)化方程

歸一化方程

MATLAB的標(biāo)準(zhǔn)化方程

運(yùn)行調(diào)試。此題考慮磁場(chǎng)沿Z相變化的情況：8（z）=80exp（z-0）并研究。變化的影

響。

8.函數(shù)

y=/-a%+1

研究參數(shù)a在[T:0.1:1]內(nèi)函數(shù)的根（實(shí)部與虛部）隨a的變化關(guān)系。（用圖形表示（exno54））

9.典型的蟲(chóng)口混沌問(wèn)題

%=一怎一嫡

X”“和XM分別為次年和當(dāng)年的蟲(chóng)口數(shù)。a=[2.03.23.53.8]xl=O.5n=50?研究系數(shù)

a的變化對(duì)x0“變化的影響，用圖形的方法，橫、縱坐標(biāo)分別為n和

10.5個(gè)學(xué)生A、B、C、D、E參加一項(xiàng)比賽。甲、乙兩觀眾猜測(cè)比賽結(jié)果。甲猜的名次順序

為A、B、C、D、E,結(jié)果一個(gè)也不對(duì)，也沒(méi)一對(duì)相鄰名次正確。乙猜的名次順序?yàn)镈、A、E、

C、B,結(jié)果猜對(duì)了兩個(gè)學(xué)生的名次并猜對(duì)了兩對(duì)學(xué)生相鄰名次，問(wèn)比賽結(jié)果？

數(shù)學(xué)模型：A、B、C、1）、E的編號(hào)為1,2,3,4,5,名次的變量為Cl,C2,C3,C4,C5

相鄰問(wèn)題DAAEECCB用41,15,53,32表示。

11.制一個(gè)函數(shù)文件要求：函數(shù)能根據(jù)輸入?yún)?shù)的個(gè)數(shù)來(lái)決定程序的走向：

a）輸入零參數(shù)時(shí)給出相關(guān)信息；

b）輸入1個(gè)參數(shù)時(shí)，以這個(gè)數(shù)為邊長(zhǎng)畫正方形；

c）輸入2個(gè)參數(shù)時(shí)，以這兩個(gè)數(shù)分別為邊長(zhǎng)畫長(zhǎng)方型；

d）輸入3個(gè)參數(shù)時(shí)，以前兩個(gè)數(shù)為坐標(biāo)原點(diǎn)，后一個(gè)數(shù)為半徑畫院。

第5章MATLAB在不同計(jì)算方法中的應(yīng)用

MATLAB擁有最強(qiáng)大的計(jì)算功能，本章主要介紹MATLAB在不同計(jì)算方法中的應(yīng)用。

5.1離散數(shù)據(jù)的擬合

我們需要經(jīng)常關(guān)心離散數(shù)據(jù)的連續(xù)性規(guī)律，包括對(duì)數(shù)據(jù)發(fā)展的趨勢(shì)分析等等，而在

MATLAB中這是比較方便的事。

5.1.1多項(xiàng)式及其運(yùn)算

1.多項(xiàng)式的創(chuàng)建

(1)由1XN的向量

aa

p=[420%2--n}

表示

勺=gx"+qx"1+出+x"一…+

多項(xiàng)式，如用poly2sym()可以查看這個(gè)多項(xiàng)式。

(2)由函數(shù)poly(A)定義

如果A為二維矩陣，poly(A)給出A的特征多項(xiàng)式。如果A為一維矩陣，poly(A)表

示由A的元素為多項(xiàng)式的根所確定的多項(xiàng)式。

【例5-1】產(chǎn)生多項(xiàng)式的方法。

clear

%方法一(山多項(xiàng)式的系數(shù)確定的多項(xiàng)式)

p=[l-23]%直接給出多項(xiàng)式p

poly2sym(p)%給出p多項(xiàng)式的表達(dá)式

%方法二(由矩陣所確定的多項(xiàng)式)

a=[l2;-24]

ps=poly(a)%給出a的特征多項(xiàng)式

poly2sym(ps)%給出ps多項(xiàng)式的表達(dá)式

%方法三(由多項(xiàng)式的根確定的多項(xiàng)式)

x=[-l2]

px=poly(x)%以x的元素為多項(xiàng)式的根確定的多項(xiàng)式。

poly2sym(px)%給出ps多項(xiàng)式的表達(dá)式

運(yùn)行結(jié)果為

P

1-23

ans=

x八2-2*x+3

a=

12

-24

ps=

1-58

ans=

xA2-5*x+8

x=

-12

px=

1-1-2

ans=

xA2-x-2

2.多項(xiàng)式函數(shù)的引用

我們可以很方便地引用多項(xiàng)式函數(shù)(即求多項(xiàng)式的函數(shù)值)

引用格式

Y=polyval(px,x)

這里，引用函數(shù)為polyval。括號(hào)中，px為多項(xiàng)式的名，x為多項(xiàng)式自變量取值，Y為對(duì)應(yīng)

的函數(shù)值。

【例5-2】多項(xiàng)式函數(shù)的引用

clear

d=[-l2]

px=poly(d)

y=polyval(px,4)%求多項(xiàng)式px在自變量等于4時(shí)的函數(shù)值

x=-4:0.5:8

yx=polyval(px,x)%求多項(xiàng)式px在自變量等于x序列時(shí)的函數(shù)值序列

plot(x,yx)%作出兩個(gè)變量的函數(shù)圖

a=roots(px)%求多項(xiàng)式px的根

運(yùn)行結(jié)果為

d=

-12

px=

1-1-2

y=

io

x=

Columns1through6

-4.0000-3.5000-3.0000-2.5000-2.0000-1.5000

Columns7through12

-1.0000-0.500000.50001.00001.5000

Columns13through18

2.00002.50003.00003.50004.00004.5000

Columns19through24

5.00005.50006.00006.50007.00007.5000

Column25

8.0000

yx=

Columns1through6

18.000013.750010.00006.75004.00001.7500

Columns7through12

0-1.2500-2.0000-2.2500-2.0000-1.2500

Columns13through18

01.75004.00006.750010.000013.7500

Columns19through24

18.000022.750028.000033.750040.000046.7500

Column25

54.0000

a=

2

-1

3.分式多項(xiàng)式的展開(kāi)

(1)傳遞函數(shù):本質(zhì)是將時(shí)域上的微分或積分方程進(jìn)行Laplace變換，結(jié)果是將時(shí)域問(wèn)

題變?yōu)轭l域問(wèn)題求解，數(shù)學(xué)變換的關(guān)鍵是

d

dt

以及

【力1

式中

S=j3

于是，傳遞函數(shù)一般是S的多項(xiàng)式。

【例5-3]求一個(gè)RC低通濾波器的幅頻與相頻特性圖和轉(zhuǎn)折頻率。

%低通r=100千歐c=l微法

x=0:100;

y=l./(j*0.l*x+l);

A=abs(y);

P=angle(y);

g=abs(A-0.707)；

[a,b]=min(g)

x0=x(b)

P0=P(b)

subplot(221)

plot(x,A)

subplot(222)

plot(x,P)

%轉(zhuǎn)折頻率為1/RC

運(yùn)行結(jié)果為

a=

1.0678e-004

b=

11

xO=

10

P0=

-0.7854

幅頻

(2)分子、分母多項(xiàng)式的單項(xiàng)展開(kāi)

留數(shù)定理：設(shè)函數(shù)在D域內(nèi)除有限個(gè)奇點(diǎn)外解析，在閉域D+C上除這些點(diǎn)外連續(xù)，則有

,/(z)dz=2%/6/?(%)

k=\

分子、分母多項(xiàng)式的單項(xiàng)展開(kāi)

在控制系統(tǒng)的分析中經(jīng)常需要將由分母、分子多項(xiàng)式構(gòu)成的傳遞函數(shù)進(jìn)行部分展開(kāi)，

如

A(s)a\alan.

—=-------+--------+......+---------+k(zs)x

B(s)s-bls-b2s-bn

這時(shí)可以用

[a,b,k]=residue(AN,BN)

來(lái)進(jìn)行分解。這里，A和B為多項(xiàng)式，a和b是展開(kāi)式的多項(xiàng)式，分別稱為留數(shù)和殘數(shù)。AN

和BN是A和B的系數(shù)。K為直行向量。這對(duì)分析函數(shù)奇點(diǎn)非常有用。

【例5-4]請(qǐng)將

(S+1)(5+2)

s(s+3)(s+4)

進(jìn)行部分分式展開(kāi)。

?AN=[132];

?BN=[17120];

?[r,p,k]=residue(AN,BN)

1.5000

-0.6667

0.1667

P=

-4

-3

0

k=

D

相當(dāng)于原式為

1.5-0.66670.1667

-------+-------------+-----------

5+4s+3

4.多項(xiàng)式的乘除與微分運(yùn)算

乘:conv（卷積）除：deconv（解卷）polyder（微分）

【例5-5】計(jì)算x（2x+3）（x+18）

clear

al=[l0];

a2=[23];

a3=[l18];

pl二conv(al,a2)

p2=conv(pl,a3)

[p3,r]=deconv(p2,a3)

conv(p3,a3)+r

運(yùn)行結(jié)果為

pl=

230

p2=

239540

p3=

230

r=

0000

ans=

239540

>>poly2sym(ans)

ans=

2*x-3+39*x-2+54*x

5,多項(xiàng)式的求根

n次多項(xiàng)式有n個(gè)根，它們可以是實(shí)數(shù)、虛數(shù)或共扼復(fù)數(shù)。MATLAB中roots用來(lái)求全

部根。如

?A=[61031]

A=

61031

?roots(A)

ans

0.4414+0.6980i

0.4414-0.6980i

-0.7006

-0.3488

5.1.2多項(xiàng)式的曲線擬合

1.用多項(xiàng)式函數(shù)去模擬?個(gè)離散數(shù)據(jù)的方法，稱為多項(xiàng)式的曲線擬合。

2.方法：

1)找出函數(shù)上的已知點(diǎn)系列。

2)由已知點(diǎn)系列確定多項(xiàng)式，即

P=polyfit(x,>1,//)

式中，P為模擬的多項(xiàng)式，polyfit為調(diào)用函數(shù)，x和y是已知點(diǎn)系列，n是多項(xiàng)式的階次。

(一般n越大越精確)

【例5-6】用多項(xiàng)式去模擬一個(gè)正弦函數(shù)

clear

x=0:0.1:6;

y=sin(x);

xx=0:6;

yy=sin(xx);

a1=polyfit(xx,yy,3)；

yl=polyval(al,x);

a2=polyfit(xx,yy,4)；

y2=polyval(a2,x);

a3=polyfit(xx,yy,5);

y3=polyval(a3,x);

subplot(231)

plot(x,y；；x,yl,?)

subplot(232)

plot(x,y；-',x,y2,7)

subplot(233)

plot(x,y；-',x,y3,7)

結(jié)果為

11

可見(jiàn)，與模擬情況與多項(xiàng)式的階次有關(guān)。

【例5-7】函數(shù)的趨勢(shì)分析

clear

loadcensus.mat

x=cdate

y=pop

plot(cdate,pop/o*)

a1=polyfit(cdate,pop,2)

a2=polyfit(cdate,pop,3)

a3=polyfil(cdate,pop,4)

a4=polyfit(cdate,pop,5)

aa1=polyval(a1,2000)

aa2=polyval(a2,2000)

aa3=polyval(a3,2000)

aa4=polyval(a4,2000)

運(yùn)行結(jié)果為

aal=

274.6221

aa2=

276.9632

aa3=

279.0252

aa4=

281.0268

200

150

100

50

圖5-7

很明顯，在2000年時(shí)得到的數(shù)據(jù)是不一樣的，應(yīng)怎樣來(lái)確定這個(gè)數(shù)，我們將在后面的

章節(jié)進(jìn)行闡述。

5.2插值

在離散數(shù)據(jù)之間根據(jù)數(shù)據(jù)的規(guī)律選擇適當(dāng)?shù)拿芏戎萌霐?shù)據(jù)的方法稱為插值。與擬和的方

法不同，插值后形成的曲線一定過(guò)已知數(shù)據(jù)點(diǎn)。

5.2.1—維內(nèi)插

一維內(nèi)插的基本語(yǔ)句

Yi=interpl(x,y,xi,)方法參數(shù)')

(Dx是原始數(shù)據(jù)自變量單增向量，y是原始數(shù)據(jù)函數(shù)向量或矩陣，若為矩陣，插值按列進(jìn)

行.

(2)xi是需要插值的自變量矩陣，yi是由插值方法計(jì)算出的函數(shù)矩陣。

(3)插值方法參數(shù)有四個(gè)。Nearset(近點(diǎn)插值)，1inear(線性)Spline(三次樣條)cublic

(內(nèi)插)

【例5-8】分析下面程序，學(xué)習(xí)插值方法和插值方法參數(shù)的影響。

x=0:2*pi

y=sin(x)

xi=0:0.1:8

yil=interpl(x,y,xi,linear')

yi2=interpl(x,y,xi,'nearest*)

yi3=interp1(x,y,xi,'spline')

yi4=interpl(x,y,xi/cubic')

p=polyfit(x,y,5)

yy=polyval(p,xi)

subplot(3,2,l)

plot(x,y,'o')

subplot(3,2,2)

plot(x,y,'o',xi,yy)

subplot(3,2,3)

plot(x,y,'o',xi,yil)

subplot(3,2,4)

plot(x,y,'o',xi,yi2)

subplot(3,2,5)

plot(x,y,'o',xi,yi3)

subplot(3,2,6)

plot(x,y；o',xi,yi4)

運(yùn)行結(jié)果為

例5-8圖

值得注意的是，擬和形成了?條連續(xù)的多項(xiàng)式曲線，而插值后得到更密集的離散數(shù)據(jù)，

另外不同插值的方法可以得到不同的離散數(shù)據(jù)群。

5.2.2二維內(nèi)插

二維插值的基本格式為

zi=interp2(x,y,z,xi,yi,'方法參數(shù)')

(l)x和y是原始數(shù)據(jù)自變量單增向量生成的二維矩陣，z為原始數(shù)據(jù)矩陣。

(2)xi和yi是需要插值的自變量向量生成的二維矩陣,zi是由插值方法計(jì)算出的函數(shù)矩

陣。

(3)插值方法參數(shù)有四個(gè)。(同前)

【例5-9】學(xué)習(xí)二維內(nèi)插數(shù)據(jù)的方法。

clear

[x,y]=meshgrid(-3:3);%建立原史數(shù)據(jù)數(shù)據(jù)點(diǎn)

z=peaks(x,y)；

[xi,yi]=meshgrid(-3:0.1:3);%給出自變量的插值點(diǎn)

zi=interp2(x,y,z,xi,yi「spline')%根據(jù)原史數(shù)據(jù)數(shù)據(jù)點(diǎn)和插值方法計(jì)算插值函數(shù)點(diǎn)

surf(x,y,z)%原史數(shù)據(jù)數(shù)據(jù)描點(diǎn)

holdon

surf(xi,yi,zi+20)%插值點(diǎn)數(shù)據(jù)數(shù)據(jù)描點(diǎn)

holdoff

axistight

運(yùn)行結(jié)果為

例5-9圖

5.2.3二維散布點(diǎn)內(nèi)插

如果自變量不是整齊分布在珊格內(nèi)，而是散布點(diǎn)分布，這時(shí)內(nèi)插函數(shù)用

zi二griddata(x,y,z,xi,yi)

其中，x,y,z是散布點(diǎn)坐標(biāo)。

【例5-10】散布點(diǎn)坐標(biāo)內(nèi)插

clear

x=6*rand(100,1)?3;%[-3,3]之間的個(gè)均勻分布的隨機(jī)數(shù)

y=5*rand(100,1)-3;

z=peaks(x,y)；

[xi,yi]=meshgrid(-3:0.2:3);

zi=griddata(x,y,z,xiJyi)；

closeall

mesh(xi,yi,zi)

holdon;plot3(x,y,z,'o');holdoff

axistight;hiddenoff;

運(yùn)行結(jié)果為

例5-10圖

另外，還有三維珊格點(diǎn)內(nèi)插

vi=interp3(x,y,z,v,xi,yi,zi,'method?)

x,y,z,v是三維矩陣，x,y,z,是坐標(biāo)，y一般由ndgrid所產(chǎn)生。v是數(shù)據(jù)(比如粒子密

度等),xi,yi,zi是內(nèi)插點(diǎn)。方法有nearst,linear,spline,cubic.

高維內(nèi)插為

vi=interpn(xl,x2,x3…，v,yl,y2,y3…'methodJ)

這些函數(shù)在讀者需要用的時(shí)候可以查相關(guān)函數(shù)說(shuō)明。

5.2.4三角內(nèi)插與計(jì)算幾何

在計(jì)算兒何(ComputationalGeometry)方面，ATLAB也提供了一?系列命令，可解決三角

化(Triangulation)、鄰近點(diǎn)(NearestNeighbom)等方面的問(wèn)題。何謂三角化(Triangulation),就

是給定--組x-y平面上的點(diǎn)，可用Delaunay命令可返回■組由這些點(diǎn)所形成的三角形，而

且任一點(diǎn)均不會(huì)位于任一個(gè)三角形的外接圓之內(nèi)。例如，先載如一組平面數(shù)據(jù)

?loadseamount;

?closeall

?plot(x,y,\,)

得到

-47.951?????????

??

??

?48?-

???4??

-48.05■.**?-

-481

-48.15

-48.2

■48.25

-48.3???-

-48.35■?-

-48.4-?-

???

-48.45---------1---------1----------1---------1----------1---------1---------1----------1---------1---------

210.8210.9211211.1211.2211.3211.4211.5211.6211.7211.8

對(duì)這些數(shù)據(jù)實(shí)現(xiàn)三角化，并將圖疊在原數(shù)據(jù)點(diǎn)上

tri=delaunay(x,y);

holdon;trimesh(tri,x,y,z);holdoff

得至U

函數(shù)功能

Delaunay執(zhí)行三角化

利用三角內(nèi)插產(chǎn)生曲面

?trisurf(tri,x,y,z)

關(guān)于計(jì)算幾何的些函數(shù)如下

Delaunayn執(zhí)行N-DDelaunay三角化

desearch搜尋三角化的最近點(diǎn)

Tsearch搜尋給定點(diǎn)的三角形

Convhull計(jì)算最小凸多邊形

Voronoi計(jì)算Voronoi圖形

Inpolygon決定…點(diǎn)是否位于多邊形的內(nèi)部

Recyint兩個(gè)(或多個(gè))長(zhǎng)方形的面積

polyarea多方形的面積

5.3回歸分析

在曲線擬合問(wèn)題中，所建立的數(shù)學(xué)模型(即多項(xiàng)式的系數(shù))是山MATLAB自動(dòng)確定的，

現(xiàn)在我們需要研究這些系數(shù)的不同取值所形成的多項(xiàng)式對(duì)模擬數(shù)值造成的誤差，這就是回歸

分析。如果輸出與這些系數(shù)成線性關(guān)系，則稱為線性回歸；如果輸出與這些系數(shù)成非線性關(guān)

系，則稱為非線性回歸。

5.3.1線性回歸

我們?nèi)砸悦绹?guó)自1790至1990年(以10年為一個(gè)位位)的總?cè)丝陔S時(shí)間變化的離散數(shù)

據(jù)為例來(lái)分析。加載文件census,mat

loadcensus.mat

x=cdate

y=pop

plot(x,y,'o')

能得到圖5-7o

由該圖直觀觀察得知，通過(guò)這些點(diǎn)的曲線可能是一條二次的拋物線，所以我們可以假設(shè)

所需的數(shù)學(xué)模型為：

2

y=f(x,a0,at,a2)=a0+atx+a2x

其中y(人口總數(shù))為此模型的輸出，x(年度)為此模型的輸入，ao、ai及a2則為此

模型的參數(shù)(Parameters)。由于這些參數(shù)相對(duì)于輸出y呈線性關(guān)系，所以此模型稱為“具有

線性參數(shù)”(Linear-in-the-paramelers)的模型，現(xiàn)在我們需要找出最佳的參數(shù)值，使得模

型輸出與實(shí)際數(shù)據(jù)越接近越好。假設(shè)我們的觀察數(shù)據(jù)可寫成(xi,yi),1=1-21?當(dāng)輸入為

xi時(shí)，實(shí)際輸出為yi,但模型的預(yù)測(cè)值為/(七,即，4，電)=4)+。1七+42后，因此平方誤

差為[yi—f(xi)]2,而總平方誤差為：

21

/■=|

顯然，E是各系數(shù)的函數(shù)。求出E對(duì)各系數(shù)的導(dǎo)數(shù)，令其等于零就可以求出各系數(shù)的值。另

外，21個(gè)點(diǎn)帶入二次方程中可以得到組線性方程組，這時(shí)，三個(gè)系數(shù)就是這個(gè)方程組對(duì)

應(yīng)的矩陣方程的未知數(shù)。只要找到一組系數(shù)使方程兩邊的差值最小，問(wèn)題就算解決。

MATLAB利用左除的方法來(lái)獲得，代碼如下：

?A=[ones(size(cdate)),cdate,cdate.A2];

?y=pop;

?theta=A\y

theta=

1.0e+004*

2.1130

-0.0024n

0.0000

這就是一組最佳的系數(shù)值。找出組系數(shù)后我們就可以確定一個(gè)二次式?？梢宰C明它與原

始數(shù)據(jù)最接近。實(shí)際上，前面講到的polyfit就是利用這個(gè)原理來(lái)編制的函數(shù)。當(dāng)然，線性

回歸的效果，與所選取的數(shù)學(xué)模型有很大的關(guān)系。模型所含的參數(shù)越多，平方誤差就越小，

如果參數(shù)個(gè)數(shù)等于數(shù)據(jù)點(diǎn)個(gè)數(shù)型，則在一般情況下，平方誤差會(huì)等于零。但這并不表示表示

預(yù)測(cè)會(huì)最準(zhǔn)，因?yàn)閿?shù)據(jù)點(diǎn)含有噪聲，完全吻合數(shù)據(jù)的模型也代表此模型受噪聲的影響最大，

預(yù)測(cè)的準(zhǔn)確度也會(huì)很差，因此，“模型復(fù)雜度”(即可變參數(shù)的個(gè)數(shù))和“預(yù)測(cè)準(zhǔn)確度”是相

互抗衡的兩個(gè)因素，兩者無(wú)法兼得，因此我們必須根據(jù)對(duì)問(wèn)題本身的了解來(lái)找到一個(gè)平衡點(diǎn)。

有些問(wèn)題，可以用“多輸入、單輸出”的線性回歸模型來(lái)表達(dá)

y=/(x)=aafQ(x)+aJi(x)+…+a?fn(x)

由此，我們可以將一個(gè)函數(shù)按照上面的方法組合成多輸入(即左、力等等)的線性。下面我

們舉例來(lái)說(shuō)明，如果關(guān)系，只要通過(guò)線性回歸分析找出上述的最佳系數(shù)就可以按相關(guān)理論來(lái)

求解。這里，我們將攜用稱為本底函數(shù)。下面我們舉例來(lái)說(shuō)明，如果本底函數(shù)含有噪聲，可

以通過(guò)線性回歸分析來(lái)重現(xiàn)本底函數(shù)。

假定本底函數(shù)是peaks函數(shù)

z=3*(1-x).A2,*exp(-(x.A2)-(y+1).八2)-10*(x/5-x.八3-y.A5).*exp(-x.A2-y.A2)-

l/3*exp(-(x+l),A2-y.A2)

上述函數(shù)可以寫為

z=3fl(x)-10f2(x)-l/3f3(x)

假如加上噪聲n,多輸入的系數(shù)未知，則

z=al*fl(x)+a2*f2(x)+a3*f3(x)+n

下面我們先看加上正態(tài)分布的噪聲n的peaks函數(shù)圖，輸入

?p=10;

?[xx,yy,zz]=peaks(p);

?zz=zz+randn(size(zz));%加上正態(tài)分布的噪聲

?surf(xx,yy,zz);axistight

可以得到下圖

顯然，由于噪聲的影響這與原函數(shù)差別很大。F面我們利用已知的數(shù)據(jù)點(diǎn)，在知道基底函數(shù)

的情況下，通過(guò)回歸分析來(lái)找出最佳的al、a2和a3去擬合曲面。

x=xx(:);

y=yy(：);

z=zz(:);

A=[(l-x).A2.*exp(-(x.A2)-(y+l).A2),(x/5-x.A3-y.A5).*exp(-x.A2-y.A2),exp(-(x+1).A2-y.A2)];

?tq=A\z

tq=

3.0794

-9.1175

0.4095

這與原函數(shù)的系數(shù)很接近！此時(shí)我們可以輸入較密集的數(shù)據(jù)點(diǎn)來(lái)得到回歸曲面。代碼如下

pp=31;

?[xx,yy]=meshgrid(linspace(-3,3,pp))Jinspace(-3,3,pp))；

?lxx,yy]=meshgrid(linspace(-3,3,pp),linspace(-3,3,pp))；

?x=xx(:);

?y=yy(：);

?A=[(l-x).A2.*exp(-(x.A2)-(y+l).A2),(x/5-x.A3-y.A5).*exp(-x.A2-y.A2),exp(-(x+l).A2-y.A2)];

?zz=reshape(A*tq,pp,pp);

?surf(xx,yy,zz);

?surf(xx,yy,zz);axistight

得到如下曲面

從以上例子可以看見(jiàn)，我們實(shí)際上是在含有噪聲的信息中提取的多輸入的系數(shù)以及本底函

數(shù)，模擬出了函數(shù)的曲面。顯然，本底函數(shù)的正確與否或者與真實(shí)情況接近的程度以及噪聲

的性質(zhì)決定了模擬曲面的正確性。本例我們是借用了正態(tài)分布的噪聲和原函數(shù)的本底函數(shù)，

使模擬得到了成功。

5.3.2非線性回歸

由于多方面的原因，非線性回歸的分析要困難的多。我們先來(lái)認(rèn)識(shí)一下。例如，數(shù)學(xué)

模型為

y=a}e+a2e~

這里，%、的是線性的，4、但友是非線性的，此數(shù)學(xué)模型是非線性的。平方誤差為

￡1(%,%,44)=+。2如')2

1-1

可以使用以下函數(shù)來(lái)處理，使E最小

FunctionE=fun_w(theza,data)

x=data(:,l);

y=data(:,2);

model_y=theta(1)*exp(theta(3)*x)+theta(2)*exp(theta(4)*x)；

E=suml((y-model_y).A2);

這里，theta是參數(shù)向量，包含了a”a?：、4和％,data則是觀察到的數(shù)據(jù)點(diǎn)，返回

的值E則是總平方誤差，要求E的最小值，我們可使用fminsearch命令，例如我們?nèi)杂们?/p>

面的人口數(shù)據(jù)

data=

1.0e+003*

1.79000.0039

1.80000.0053

1.81000.0072

1.82000.0096

1.83000.0129

1.84000.0171

1.85000.0231

1.86000.0314

1.87000.0386

1.88000.0502

1.89000.0629

1.90000.0760

1.91000.0920

1.92000.1057

1.93000.1228

1.94000.1317

1.95000.1507

1.96000.1790

1.97000.2050

1.98000.2265

1.99000.2487

結(jié)合下面的代碼

loadcensus.mat

data=[cdatepop]

s=[0000];

theta=fminsearch(,fun_w\s,[],data);

x=data(:,l);

y=data(:,2);

esti=theta(1)*exp(theta(3)*x)+theta(2)*exp(theta(4)*x);

pk>t(x,y,To',x,esti,'b-');

legendCsmpledata1,'regression',turve')

就能得到回歸曲線。當(dāng)然，這是一個(gè)一般方法，深入的研究讀者可以查相關(guān)文獻(xiàn)。

5.4常微分方程

常微分方程(OrdinaryDifferentialEquation(ODE))一般情況下是沒(méi)有解析解的,

通常只能是數(shù)值求解，現(xiàn)在我們來(lái)研究這一問(wèn)題。在MATLAB中，非剛性常微分方程(組)

通常用ode23和0de45這兩個(gè)函數(shù)去求解，它們分別分別采用了二階三級(jí)和四階五級(jí)RKF

方法，即為變步長(zhǎng)積分。

5.4.1ODE函數(shù)的一般用法

(1)建立標(biāo)準(zhǔn)矩陣微分方程(組)

1=F(y,t)

dt

目的：將所有不同變量歸為同一變量矩陣。高階微分方程歸為一階微分方程(組)。例如,

二階微分方程

y”+2y-y=O

令

M=>

y2=/

原式變?yōu)橐浑A微分方程(組)

>2=M-2y2

(2)建立ODE相應(yīng)的函數(shù)M文件(寫微分方程)，格式如

functionff=fun(t,y)

ff=[F(y,t)]

(3)ODE的調(diào)用格式

可在命令窗進(jìn)行也可在M文件中進(jìn)行。

［t,y］=ode45('fun',［自變量的范圍］，［初值列陣］)

【例5-11】求解

—=—2y+2x2+lxy(0)=1

dx

微分方程的函數(shù)文件

functionyy=fun1(x,y)

yy=[-2*y+2*x.A2+2*x]

end

調(diào)用該函數(shù)求解微分方程

clear

[x,y]=ode23('funl',[0,0.51,1)

plot(x,y)

結(jié)果為

095

09■\-

085■\-

08■\-

075■\■

07?-

065■-

06)-------1--------1--------1----------1--------1----------1------------1------1------------1------

00050101502025030350404505

X

例5Tl圖

【例5-12]求解

:=x-0.01xy

dy=-y+0.02盯

dt

x(0)=30

y(0)=20

這是一個(gè)微分方程組，微分方程組的函數(shù)文件

functionf=fun2(t,y)

f=[y(l)-0.01*y(l)*y(2);-y(2)+0.02*y⑴*y(2)]

end

這里，y(l)和y(2)分別表示x和y,調(diào)用該函數(shù)求解微分方程的代碼

clear

ts=0:0.1:20

y0=[30,20];

[t,y]=ode45(，fun2*,ts,yO);

subplot⑵2,1)

plot(t,y(:,1),'r',t,y(:,2),'b'),gtext('y(l)'),gtext('y(2)')

subplot(2,2,2)

plot(y(:,1),y(:,2)),gtext(*y⑴'),gtext(*y⑵')

結(jié)果為

400400