matlab基础及应用(教程全书)part2

第4章程序设计

在前面我们已经看到，MATLAB不但可以在命令窗直接输入命令并运行，而且还可以生

成自己的程序文件，这就是我们通常说的一类以M为后缀的M文件，本章我们就来研究这类

文件的形成方法。

M文件可分分为两大类，一是命令式M文件（也称为脚本文件，script）,二是函数式M

文件（function）。两类文件的区别在于：

（1）命令式文件可以直接运行，函数式文件不能直接运行，只能调用。

（2）命令式文件运行时没有输入输出参量，函数式文件在调用时需要进行输入输出

参量设置。

（3）命令式文件运行中可以调用工作空间的数据，运行中产生的所有变量为全局变

量。

（4）函数式文件不能调用工作空间的数据，运行中产生的所有变量为局部变量。命

令式文件运行中产生的所有变量为全局变量，可以调用和存储到工作空间的数据。

4.1MATLAB的程序文件-M文件

4.1.1脚本文件（Scripts）

当我们需要在命令窗进行大量的命令集合运行时,直接从命令窗口输入比较麻烦，这时

就可以将这些命令集合存放在一个脚本文件（Scripts）中，运行时只需要输入其文件名就

可以自动执行这些命令集合。需要注意的是，脚本文件运行所产生的变量都驻留在MATLAB

的工作空间中，同时脚本文件也可以调用工作空间中的数据。因此，脚本文件所涉及的变量

是全局变量。前儿章所涉及到的M文件都是这类脚本文件。

编辑一个脚本文件可以直接在命令窗口的左上角打开编辑窗进行编辑。

4.1.2函数文件（function）

函数式文件（function）的构成

（1）函数定义行：

Function［输出参量］=gauss（输入参量）

（2）：

完成函数的功能。

（3）函数说明。

（4）函数行注。

从上面构成的情况看，函数式文件实际上是完成输入参量与输出参量的转换，这样的转换是

由函数文件名为gauss的文件来完成的。函数体的功能必须说明清楚输入参量与输出参量的

关系。函数说明是用来解释该函数的功能的，函数行注是对程序行进行说明的。上面（1）

和（2）是必须的。

【例47】分析下面函数文件。

%一个数列，任意项等于前两项之和，输入项数可以给出这个数列

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

该函数文件的文件名为sul.m,在第一行给出了该函数的功能，即输入项数就可以自动给出

一个满足条件的数列。定义行规定了输入参数是该数列的项数，输出参数是该数列。从第三

行起，是该函数的主体，主要说明了输入参数与输出参数的关系。当我们在命令窗或脚本文

件中调用该函数就会有结果，如

»p=sul(5)

P=

12358

【例4-2】分析下面函数文件。

%一个数列的通项公式为a(n+1尸a(n)+M2,给定任意项的值，求这个数列的后10项，并画

离散函数图

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)

结果为

b=

10

x=

567891011121314

ans=

Columns1through5

1011152440

Columns6through10

65101150214295

Columns11through15

3955166608291025

以及

4.2程序的流程控制

4.2.1循环控制

1.硬循环语句(for-end)

所谓硬循环是指无条件的循环，其结构为

for(循环变量)

end

【例4-3］用硬循环语句形成一个5阶对角线元素为2的矩阵。

clear

n=5;%每个for和if都必须对应一个end

fori=l:n%i是循环变量

forj=l:n

ifi==j

a(i,j)=2;

else

a(i,j)-0;

end

end

end

结果为

a=

20000

02000

00200

00020

00002

根据结果，如果要两条主对角线元素都为2,如何调此程序?

注意，硬循环语句也可用break来跳出循环。

4.2.2条件控制

1.条件分支语句(if-else-end)如果-那么-否则

结构

If条件1

语句1

elseif条件2

语句2

else

语句

end

注意这一结构的条件优先问题，当有多个条件时，如果条件1不满足，再判断elseif

后面的条件2…，如果所有的条件都不满足，在则执行else后面的语句。

【例4-5】分析下面条件语句。

x=input('x=?')

if((x+3)<2)

y=x*2;

elseif(x<2)

y二x'2;

else(x<l)

y=sqrt(x);

end

y

运行结果

x=?l.5

x=

1.5000

y=

2.2500

2.条件循环语句(while-end)

当条件满足时执行下一语句，否则就到end返回while。其结构为

while(表达式)

OOO

end

【例4-4】求阶乘大于或等于99-99的最小整数。

clear

n=l;

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

n=n+l；

end

n

结果为

n=

120

4.2.3开关与试探控制

分支开关语句(switch-case-otherwise-end)

结构switch(开关量)，相当于满足什么条件做什么事。

【例4-5］根据开关量的情况来决定计算机计算机执行那个语句。

a=input('a=?f)

switcha

case1

disp('Itisraning1)

case0

disp('Itdonotkonw,)

case-1

disp(,Itisnotraning，)

otherwise

disp('Itisraning?，)

end

4.3程序设计

4.3.1程序的调试

1.程序的错误一般分为语法和逻辑两种错型。语法错误时程序无法运行，并在命令窗给予

提示。而逻辑错误只能根据计算的具体情况来判断。

2.根据出错信息调试

根据命令窗的提示，注意一般情况不加号调试。

3.设置断点来判断

dbstop格式：dbstopinM文件名at行号(也可用鼠标)

4.利用keyboard命令在文件中来判断

当出现k»时returrio

5.变量的鼠标观测法

当程序代码运行后，每个变量的数据都可以用鼠标放在相应的位置处去观察该变量的数据情

况，以判断程序运行的具体执行情况。

6.代码运行的计时方法

(1)整段程序代码的计时

tic…toe表示计算tic与toe之间的时间

【例4-6】计算机程序运行时间计算。

tic

a=rand(300);

inv(a);%计算逆矩阵

toe

%注：程序中不需要显示结果的就不显示，这可以节约机时。

结果为

»

elapsedtime=

0.3900

(2)也可以用来计算tl,t2之间的时间差来完成上述功能。

【例4-7】计算机程序运行时间计算。

t0=clock

a=rand(300);

inv(a);%计算逆矩阵

elapsedtime=etime(clock,tO)

%注：程序中不需要显示结果的就不显示，这可以节约机时。

结果仍为

»

elapsedtime=

0.3900

(3)也可以用eputime变量来完成上述功能。

【例4-8】计算机程序运行时间计算。

t0=eputime

a=rand(300);

inv(a);%计算逆矩阵

e1apsed_time=cputime-t0

%注：程序中不需要显示结果的就不显示，这可以节约机时。

结果为

to=

9.5787e+003

elapsed_time=

0.2710

4.3.2M文件的性能优化

1.程序代码的向量化

【例4-8］比较下面两个程序段，观察其运行的时间。

程序1

t0二eputime

n=100000;

total=0;

fori=l:n

total=total+l/i

end

total

tl=cputime-tO

运行结果为

total=

12.0901

tl=

3.1950

程序2

tic

n=100000;

a=l:n;

total二sum(1./a)

toe

运行结果为

total=

12.0901

elapsedtime=

0.0200

发现，程序1和程序2运行结果完全相同，但运行时间程序1是程序2的160倍！

2.用矩阵结构进行运算

一般情况下，完全采用矩阵运行的方式，MATLAB的程序与C语言的运行基本相同。这

必须对矩阵非常熟练，例如

x=[l23;121]

a=[456]

如果希望将a中的每一个元素乘以x的每一列，用diag(x)。

3.阵的预先配置

在一些大的程序运行时，为了变量不至于让计算机不断向内存要空间位置而耽误时间，

运行前需要给变量留出•定的空间。

【例4-9]比较下面两个程序段，观察其运行的时间。

程序1

tic

a=[l23;456;789];

fori=l:100

y(i)=det(a"i)

end

toe

结果为

elapsed_time=

0.1100

程序2

tic

a=[l23;456;789];

y=zeros(l,100);%分配内存

fori=l:100

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

end

toe

结果为

elapsedtime=

0.0100

练习4

i.已知函数

y=/一ax+1

请研究参数a在11:0.1:1］内函数的根(实部与虚部)随a的变化关系。(用图形表示)

2.编制•个函数，函数的输入参数为一个任意矩阵或向量，输出参数为该矩阵中不相同的

元素个数。

3.编制•个函数文件，函数的输入参数为一个任意矢性场

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

和这个矢性场自变量的取值范围，要求这个函数能自动画出矢性场的平面场分布。

4.设计一个函数文件，要求输入一个任意方阵，该函数能统计出方阵中大于零、等于零和

小于零的元素个数。

5.设计一个能实现下面功能的函数文件：输入两个参数：(1)任意2Xn矩阵(n>8,且第一

行为自变量，第二行为对应的函数值)；(2)介于任意两自变量之间的值。输出为该自变量

对应的函数值。(用9次多项式拟和)

6.一个MATLAB的递归函数fibo.m来计算fibo数列，定义如下：

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

此数列的初条件为

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

n的最大数为100,要求：

(1)保存你的fibo.m文件，当在命令窗调用fibo函数时，不论输入任何整数有正确的

输出。

(2)做出fibo的二维离散函数图，n取1到10,图的函数值处用小圆圈并涂为黑色，请

保存你的图形。

(3)用三次样条插值的方法对(2)中的10个点进行插值，自变量的分辨率为0.01,

请保存你的图形。

(4)将完成(3)工作的插值函数保存为fib.m文件，，当在命令窗调用fib函数时，不论输

入任何具有两位小数且小于10大于0的数(如5.45)时有正确的输出。

7.设电子粒子束流从恒定磁场中某点以相同速率发射，发射的方向与磁场方向的夹角很小,

观察不同方向入射的粒子束流的运动轨道。(设磁场沿Z方向)

数学模型：

粒子流的速度初值为

匕=%

匕=%6sin(a)

vv=%6cos(a)

标准化方程

归一化方程

MATLAB的标准化方程

运行调试。此题考虑磁场沿Z相变化的情况：8（z）=80exp（z-0）并研究。变化的影

响。

8.函数

y=/-a%+1

研究参数a在[T:0.1:1]内函数的根（实部与虚部）随a的变化关系。（用图形表示（exno54））

9.典型的虫口混沌问题

%=一怎一嫡

X”“和XM分别为次年和当年的虫口数。a=[2.03.23.53.8]xl=O.5n=50„研究系数

a的变化对x0“变化的影响，用图形的方法，横、纵坐标分别为n和

10.5个学生A、B、C、D、E参加一项比赛。甲、乙两观众猜测比赛结果。甲猜的名次顺序

为A、B、C、D、E,结果一个也不对，也没一对相邻名次正确。乙猜的名次顺序为D、A、E、

C、B,结果猜对了两个学生的名次并猜对了两对学生相邻名次，问比赛结果？

数学模型：A、B、C、1）、E的编号为1,2,3,4,5,名次的变量为Cl,C2,C3,C4,C5

相邻问题DAAEECCB用41,15,53,32表示。

11.制一个函数文件要求：函数能根据输入参数的个数来决定程序的走向：

a）输入零参数时给出相关信息；

b）输入1个参数时，以这个数为边长画正方形；

c）输入2个参数时，以这两个数分别为边长画长方型；

d）输入3个参数时，以前两个数为坐标原点，后一个数为半径画院。

第5章MATLAB在不同计算方法中的应用

MATLAB拥有最强大的计算功能，本章主要介绍MATLAB在不同计算方法中的应用。

5.1离散数据的拟合

我们需要经常关心离散数据的连续性规律，包括对数据发展的趋势分析等等，而在

MATLAB中这是比较方便的事。

5.1.1多项式及其运算

1.多项式的创建

(1)由1XN的向量

aa

p=[420%2--n}

表示

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

多项式，如用poly2sym()可以查看这个多项式。

(2)由函数poly(A)定义

如果A为二维矩阵，poly(A)给出A的特征多项式。如果A为一维矩阵，poly(A)表

示由A的元素为多项式的根所确定的多项式。

【例5-1】产生多项式的方法。

clear

%方法一(山多项式的系数确定的多项式)

p=[l-23]%直接给出多项式p

poly2sym(p)%给出p多项式的表达式

%方法二(由矩阵所确定的多项式)

a=[l2;-24]

ps=poly(a)%给出a的特征多项式

poly2sym(ps)%给出ps多项式的表达式

%方法三(由多项式的根确定的多项式)

x=[-l2]

px=poly(x)%以x的元素为多项式的根确定的多项式。

poly2sym(px)%给出ps多项式的表达式

运行结果为

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.多项式函数的引用

我们可以很方便地引用多项式函数(即求多项式的函数值)

引用格式

Y=polyval(px,x)

这里，引用函数为polyval。括号中，px为多项式的名，x为多项式自变量取值，Y为对应

的函数值。

【例5-2】多项式函数的引用

clear

d=[-l2]

px=poly(d)

y=polyval(px,4)%求多项式px在自变量等于4时的函数值

x=-4:0.5:8

yx=polyval(px,x)%求多项式px在自变量等于x序列时的函数值序列

plot(x,yx)%作出两个变量的函数图

a=roots(px)%求多项式px的根

运行结果为

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.分式多项式的展开

(1)传递函数:本质是将时域上的微分或积分方程进行Laplace变换，结果是将时域问

题变为频域问题求解，数学变换的关键是

d

dt

以及

【力1

式中

S=j3

于是，传递函数一般是S的多项式。

【例5-3]求一个RC低通滤波器的幅频与相频特性图和转折频率。

%低通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)

%转折频率为1/RC

运行结果为

a=

1.0678e-004

b=

11

xO=

10

P0=

-0.7854

幅频

(2)分子、分母多项式的单项展开

留数定理：设函数在D域内除有限个奇点外解析，在闭域D+C上除这些点外连续，则有

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

k=\

分子、分母多项式的单项展开

在控制系统的分析中经常需要将由分母、分子多项式构成的传递函数进行部分展开，

如

A(s)a\alan.

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

B(s)s-bls-b2s-bn

这时可以用

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

来进行分解。这里，A和B为多项式，a和b是展开式的多项式，分别称为留数和残数。AN

和BN是A和B的系数。K为直行向量。这对分析函数奇点非常有用。

【例5-4]请将

(S+1)(5+2)

s(s+3)(s+4)

进行部分分式展开。

»AN=[132];

»BN=[17120];

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

1.5000

-0.6667

0.1667

P=

-4

-3

0

k=

D

相当于原式为

1.5-0.66670.1667

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

5+4s+3

4.多项式的乘除与微分运算

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

【例5-5】计算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

运行结果为

pl=

230

p2=

239540

p3=

230

r=

0000

ans=

239540

>>poly2sym(ans)

ans=

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

5,多项式的求根

n次多项式有n个根，它们可以是实数、虚数或共扼复数。MATLAB中roots用来求全

部根。如

»A=[61031]

A=

61031

»roots(A)

ans

0.4414+0.6980i

0.4414-0.6980i

-0.7006

-0.3488

5.1.2多项式的曲线拟合

1.用多项式函数去模拟•个离散数据的方法，称为多项式的曲线拟合。

2.方法：

1)找出函数上的已知点系列。

2)由已知点系列确定多项式，即

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

式中，P为模拟的多项式，polyfit为调用函数，x和y是已知点系列，n是多项式的阶次。

(一般n越大越精确)

【例5-6】用多项式去模拟一个正弦函数

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)

结果为

11

可见，与模拟情况与多项式的阶次有关。

【例5-7】函数的趋势分析

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)

运行结果为

aal=

274.6221

aa2=

276.9632

aa3=

279.0252

aa4=

281.0268

200

150

100

50

图5-7

很明显，在2000年时得到的数据是不一样的，应怎样来确定这个数，我们将在后面的

章节进行阐述。

5.2插值

在离散数据之间根据数据的规律选择适当的密度置入数据的方法称为插值。与拟和的方

法不同，插值后形成的曲线一定过已知数据点。

5.2.1—维内插

一维内插的基本语句

Yi=interpl(x,y,xi,)方法参数')

(Dx是原始数据自变量单增向量，y是原始数据函数向量或矩阵，若为矩阵，插值按列进

行.

(2)xi是需要插值的自变量矩阵，yi是由插值方法计算出的函数矩阵。

(3)插值方法参数有四个。Nearset(近点插值)，1inear(线性)Spline(三次样条)cublic

(内插)

【例5-8】分析下面程序，学习插值方法和插值方法参数的影响。

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)

运行结果为

例5-8图

值得注意的是，拟和形成了•条连续的多项式曲线，而插值后得到更密集的离散数据，

另外不同插值的方法可以得到不同的离散数据群。

5.2.2二维内插

二维插值的基本格式为

zi=interp2(x,y,z,xi,yi,'方法参数')

(l)x和y是原始数据自变量单增向量生成的二维矩阵，z为原始数据矩阵。

(2)xi和yi是需要插值的自变量向量生成的二维矩阵,zi是由插值方法计算出的函数矩

阵。

(3)插值方法参数有四个。(同前)

【例5-9】学习二维内插数据的方法。

clear

[x,y]=meshgrid(-3:3);%建立原史数据数据点

z=peaks(x,y)；

[xi,yi]=meshgrid(-3:0.1:3);%给出自变量的插值点

zi=interp2(x,y,z,xi,yi「spline')%根据原史数据数据点和插值方法计算插值函数点

surf(x,y,z)%原史数据数据描点

holdon

surf(xi,yi,zi+20)%插值点数据数据描点

holdoff

axistight

运行结果为

例5-9图

5.2.3二维散布点内插

如果自变量不是整齐分布在珊格内，而是散布点分布，这时内插函数用

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

其中，x,y,z是散布点坐标。

【例5-10】散布点坐标内插

clear

x=6*rand(100,1)・3;%[-3,3]之间的个均匀分布的随机数

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;

运行结果为

例5-10图

另外，还有三维珊格点内插

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

x,y,z,v是三维矩阵，x,y,z,是坐标，y一般由ndgrid所产生。v是数据(比如粒子密

度等),xi,yi,zi是内插点。方法有nearst,linear,spline,cubic.

高维内插为

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

这些函数在读者需要用的时候可以查相关函数说明。

5.2.4三角内插与计算几何

在计算儿何(ComputationalGeometry)方面，ATLAB也提供了一•系列命令，可解决三角

化(Triangulation)、邻近点(NearestNeighbom)等方面的问题。何谓三角化(Triangulation),就

是给定--组x-y平面上的点，可用Delaunay命令可返回■组由这些点所形成的三角形，而

且任一点均不会位于任一个三角形的外接圆之内。例如，先载如一组平面数据

»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

对这些数据实现三角化，并将图叠在原数据点上

tri=delaunay(x,y);

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

得至U

函数功能

Delaunay执行三角化

利用三角内插产生曲面

»trisurf(tri,x,y,z)

关于计算几何的些函数如下

Delaunayn执行N-DDelaunay三角化

desearch搜寻三角化的最近点

Tsearch搜寻给定点的三角形

Convhull计算最小凸多边形

Voronoi计算Voronoi图形

Inpolygon决定…点是否位于多边形的内部

Recyint两个(或多个)长方形的面积

polyarea多方形的面积

5.3回归分析

在曲线拟合问题中，所建立的数学模型(即多项式的系数)是山MATLAB自动确定的，

现在我们需要研究这些系数的不同取值所形成的多项式对模拟数值造成的误差，这就是回归

分析。如果输出与这些系数成线性关系，则称为线性回归；如果输出与这些系数成非线性关

系，则称为非线性回归。

5.3.1线性回归

我们仍以美国自1790至1990年(以10年为一个位位)的总人口随时间变化的离散数

据为例来分析。加载文件census,mat

loadcensus.mat

x=cdate

y=pop

plot(x,y,'o')

能得到图5-7o

由该图直观观察得知，通过这些点的曲线可能是一条二次的抛物线，所以我们可以假设

所需的数学模型为：

2

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

其中y(人口总数)为此模型的输出，x(年度)为此模型的输入，ao、ai及a2则为此

模型的参数(Parameters)。由于这些参数相对于输出y呈线性关系，所以此模型称为“具有

线性参数”(Linear-in-the-paramelers)的模型，现在我们需要找出最佳的参数值，使得模

型输出与实际数据越接近越好。假设我们的观察数据可写成(xi,yi),1=1-21«当输入为

xi时，实际输出为yi,但模型的预测值为/(七,即，4，电)=4)+。1七+42后，因此平方误

差为[yi—f(xi)]2,而总平方误差为：

21

/■=|

显然，E是各系数的函数。求出E对各系数的导数，令其等于零就可以求出各系数的值。另

外，21个点带入二次方程中可以得到组线性方程组，这时，三个系数就是这个方程组对

应的矩阵方程的未知数。只要找到一组系数使方程两边的差值最小，问题就算解决。

MATLAB利用左除的方法来获得，代码如下：

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

»y=pop;

»theta=A\y

theta=

1.0e+004*

2.1130

-0.0024n

0.0000

这就是一组最佳的系数值。找出组系数后我们就可以确定一个二次式。可以证明它与原

始数据最接近。实际上，前面讲到的polyfit就是利用这个原理来编制的函数。当然，线性

回归的效果，与所选取的数学模型有很大的关系。模型所含的参数越多，平方误差就越小，

如果参数个数等于数据点个数型，则在一般情况下，平方误差会等于零。但这并不表示表示

预测会最准，因为数据点含有噪声，完全吻合数据的模型也代表此模型受噪声的影响最大，

预测的准确度也会很差，因此，“模型复杂度”(即可变参数的个数)和“预测准确度”是相

互抗衡的两个因素，两者无法兼得，因此我们必须根据对问题本身的了解来找到一个平衡点。

有些问题，可以用“多输入、单输出”的线性回归模型来表达

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

由此，我们可以将一个函数按照上面的方法组合成多输入(即左、力等等)的线性。下面我

们举例来说明，如果关系，只要通过线性回归分析找出上述的最佳系数就可以按相关理论来

求解。这里，我们将携用称为本底函数。下面我们举例来说明，如果本底函数含有噪声，可

以通过线性回归分析来重现本底函数。

假定本底函数是peaks函数

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)

上述函数可以写为

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

假如加上噪声n,多输入的系数未知，则

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

下面我们先看加上正态分布的噪声n的peaks函数图，输入

»p=10;

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

»zz=zz+randn(size(zz));%加上正态分布的噪声

»surf(xx,yy,zz);axistight

可以得到下图

显然，由于噪声的影响这与原函数差别很大。F面我们利用已知的数据点，在知道基底函数

的情况下，通过回归分析来找出最佳的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

这与原函数的系数很接近！此时我们可以输入较密集的数据点来得到回归曲面。代码如下

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

得到如下曲面

从以上例子可以看见，我们实际上是在含有噪声的信息中提取的多输入的系数以及本底函

数，模拟出了函数的曲面。显然，本底函数的正确与否或者与真实情况接近的程度以及噪声

的性质决定了模拟曲面的正确性。本例我们是借用了正态分布的噪声和原函数的本底函数，

使模拟得到了成功。

5.3.2非线性回归

由于多方面的原因，非线性回归的分析要困难的多。我们先来认识一下。例如，数学

模型为

y=a}e+a2e~

这里，%、的是线性的，4、但友是非线性的，此数学模型是非线性的。平方误差为

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

1-1

可以使用以下函数来处理，使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是参数向量，包含了a”a?：、4和％,data则是观察到的数据点，返回

的值E则是总平方误差，要求E的最小值，我们可使用fminsearch命令，例如我们仍用前

面的人口数据

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

结合下面的代码

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')

就能得到回归曲线。当然，这是一个一般方法，深入的研究读者可以查相关文献。

5.4常微分方程

常微分方程(OrdinaryDifferentialEquation(ODE))一般情况下是没有解析解的,

通常只能是数值求解，现在我们来研究这一问题。在MATLAB中，非刚性常微分方程(组)

通常用ode23和0de45这两个函数去求解，它们分别分别采用了二阶三级和四阶五级RKF

方法，即为变步长积分。

5.4.1ODE函数的一般用法

(1)建立标准矩阵微分方程(组)

1=F(y,t)

dt

目的：将所有不同变量归为同一变量矩阵。高阶微分方程归为一阶微分方程(组)。例如,

二阶微分方程

y”+2y-y=O

令

M=>

y2=/

原式变为一阶微分方程(组)

>2=M-2y2

(2)建立ODE相应的函数M文件(写微分方程)，格式如

functionff=fun(t,y)

ff=[F(y,t)]

(3)ODE的调用格式

可在命令窗进行也可在M文件中进行。

［t,y］=ode45('fun',［自变量的范围］，［初值列阵］)

【例5-11】求解

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

dx

微分方程的函数文件

functionyy=fun1(x,y)

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

end

调用该函数求解微分方程

clear

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

plot(x,y)

结果为

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

这是一个微分方程组，微分方程组的函数文件

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,调用该函数求解微分方程的代码

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⑵')

结果为

400400