matlab教程

发布时间：2023-06-13 作者：admin 来源：文学

matlab教程

2023年3月18日发(作者：南京中医药大学网络教学平台)

MatLab简介

MATLAB是什么？

典型的使用包括：

数学和计算

算术发展模型，

模拟，和原型

数据分析，开发，和可视化

科学和工程图学

应用发展包括图形用户界面设计

MATLAB表示矩阵实验室。

MATLAB系统

MATLAB系统由5主要的部分构成：

语言。这是高阶的矩阵/数组语言,带控制流动陈述，函数，数据结构，输入/输出，而且面向

目标的编程特点。

Ops操作符和特殊字符。

Lang程序设计语言作。

strfun字符串。

iofun输入/输出。

timefun时期和标有日期。

datatypes数据类型和结构。

工作环境。这是你作为MATLAB用户或程序编制员的一套工具和设施。

3.制图这是MATLAB制图系统。它为2维上，而且三维的数据可视化，图象处理，动画片制作和表示图形

包括高阶的指令在内。它也为包括低阶的指令在内,允许你建造完整的图形用户界面（GUIs），MATLAB应

用。制图法功能在MATLAB工具箱中被组织成5文件夹：

graph2d2-的维数上的图表。

graph3d三维的图表。

specgraph专业化图表。

graphics制图法。

uitools图形用户界面工具。

的数学的函数库。数学和分析的功能在MATLAB工具箱中被组织成8文件夹。

elmat初步矩阵，和矩阵操作。

elfun初步的数学函数。

specfun专门的数学函数。

matfun矩阵函数－用数字表示的线性的代数。

datafun数据分析和傅立叶变换。

polyfun插入物，并且多项式。

funfun功能函数。

sparfun稀少矩阵。

应用程序接口（API）。这是允许你写C、Fortran语言与MATLAB交互。

关于Simulink

Simulink？MATLAB为做非线性的动态的系统的模拟实验的交互式的系统。它是允许你通过把方框图拉到

屏幕，灵活地窜改它制作系统的模型的用图表示的鼠标驱动的程序。实时工作室？允许你产生来自你的图

表块的C代码，使之能用于各种实时系统。

关于工具箱

工具箱是为了解答特别种类的问题扩展MATLAB环境的MATLAB函数的综合的（M-文件）收集

MatLab工作环境

命令窗口

若输入

A=[123;456;7810]

按下回车键后显示如下

123

456

7810

清除命令窗口

clc

这并不清除工作间，只是清除了显示，仍可按上箭头看到以前发出的命令

数据格式命令

x=[4/31.2345e–6]

formatshort

1.33330.0000

formatshorte

1.3333e+0001.2345e–006

formatshortg

1.33331.2345e–006

formatlong

1.333333333333330.00

formatlonge

1.333333333333333e+0001.2345e–006

formatlongg

1.333333333333331.2345e–006

formatbank

1.330.00

format+

formatrat

4/31/810045

formathex

3ff55555555555553eb4b6231abfd271

若最大的元素大于1000或小于0.001,则显示short或long格式时时会加上一个比例

还有两个格式:

formatcompact

formatloose

禁止结果的显示

在命令后加上分号,则屏幕上不会立即显示出结果,这在运算大的数据量时十分有用,如下命令产生100*100

的幻方:

A=magic(100);

长命令行

如想另起一行输入命令,在末尾加上"..."即可,如:

s=1–1/2+1/3–1/4+1/5–1/6+1/7...

–1/8+1/9–1/10+1/11–1/12;

MatLab工作间

你可用who或whos来察看当前工作间中有哪些变量,如:

whos

NameSizeBytesClass

A4x4128doublearray

D3x5120doublearray

M10x140cellarray

S1x3628structarray

h1x1122chararray

n1x18doublearray

s1x510chararray

v1x1428chararray

Grandtotalis93elementsusing984bytes

若要从工作间中删除所有的变量,用

clear

保存、重载工作间

你可以将工作间保存为一个二进制的M文件，以后还可以恢复回来：

savejune10

也可只保存工作间中的部分变量值

savejune10xyz

重载时只需输入

loadjune10

文件名保存在字符串中

这样可以像调用函数一样调用工作间

save(’myfile’,’VAR1’,’VAR2’)

A=’myfile’;

load(A)

与下面的命令相同

savemyfileVAR1VAR2

loadmyfile

下面的命令把1至10的平方值分别存放在data1至data10中：

file=’data’;

fori=1:10j=i.^2;

save([fileint2str(i)],’j’);

end

查找路径

当你输入“yourpig"时发生了什么呢?

1:察看是否是变量;

2:察看是否是内建函数;

3:察看当前目录下是否有文件:yourpig.m;

4:察看查找目录下是否有文件:yourpig.m;

对于查找路径中的文件,what显示当前目录下的文件,加上路径后可显示输入的路径下所有的MatLab文件.

如:

whatmatlab/elfun

以下二命令分别显示、编辑m文件

typerank

editrank

图像窗口

下面的命令产生一个与命令窗口隔离的图形窗口,

figure

plot函数则会在新的窗口中绘制图形,如

t=0:pi/100:2*pi;

y=sin(t);

plot(t,y)

则有如下图形:

寻求帮助

下面的函数在寻求帮助时十分有用:

help列出你所寻求帮助的函数的功能描述；

lookfor列出所有函数的功能描述中含有你所输入的内容的函数的简介

如:

helpinverse

显示

und.

但如输入

lookforinverse

则显示

INVHILBInverseHilbertmatrix

ACOSHInversehyperboliccosine

ERFINVInverseoftheerrorfunction

INVMatrixinverse

PINVPseudoinverse

IFFTInversediscreteFouriertransform

IFFT2Two–dimensionalinversediscreteFouriertransform

ICCEPSInversecomplexcepstrum

IDCTInversediscretecosinetransform

数据分析和统计

面向列的数据集

这年头似乎十分风行”面向”这个词，这儿故也套用，其英文为"Column-OrientedDataSets",可理

解为MatLab按列的存储方式来分析数据,下面是一个例子:

TimeLocation1Location2Location3

01h0011119

02h0071311

03h00141720

04h0011139

05h00435169

06h00384676

07h

08h

09h003888115

10h00283655

11h00121214

12h00182730

13h00181929

14h00171518

15h00193648

16h00324710

17h00426592

18h005766151

19h00445590

20h

21h00355868

22h00111215

23h0013915

24h001097

以上数据被保存在一个称为的文件中．

11119

71311

141720

11139

435169

384676

61132186

75135180

3888115

283655

121214

182730

181929

171518

193648

324710

426592

5766151

445590

114145257

355868

111215

13915

1097

下面，我们调入此文件，并看看文件的一些参数

[n,p]=size(count)

创建一个时间轴后，我们可以把图画出来：

t=1:n;

set(0,'defaultaxeslinestyleorder’,’-|--|-.’)

set(0,'defaultaxescolororder’,[000])

plot(t,count),legend('Location1','Location2','Location3',0)

xlabel('Time'),ylabel('VehicleCount'),gridon

足以证明，以上是对３个对象的２４次观测．

基本数据分析函数

（一定注意是面向列的）

继续用上面的数据，其每列最大值．均值．及偏差分别为：

mx=max(count)

mu=mean(count)

sigma=std(count)

mx=

114145257

mu=

32.000046.541765.5833

sigma=

25.370341.405768.0281

重载函数，还可以定位出最大．最小值的位置

[mx,indx]=min(count)

mx=

797

indx=

22324

试试看，你能看懂下面的命令是干什么的吗？

[n,p]=size(count)

e=ones(n,1)

x=count–e*mu

点这看看答案！

下面这句命令则找出了整个矩阵的最小值:

min(count(:))

ans=

协方差及相关系数

下面,我们来看看第一列的方差:

cov(count(:,1))

ans=

643.6522

cov()函数作用于矩阵,则会计算其协方差矩阵.

corrcoef()用于计算相关系数,如:

corrcoef(count)

ans=

1.00000.93310.9599

0.93311.00000.9553

0.95990.95531.0000

数据的预处理

未知数据

NaN(NotaNumber--不是一个数)被定义为未经定义的算式的结果,如0/0.在处理数据中,NaN常用来表示

未知数据或未能获得的数据.所有与NaN有关的运算其结果都是NaN.

a=magic(3);

a(2,2)=NaN

816

3NaN7

492

sum(a)

ans=

15NaN15

在做统计时，常需要将NaN转化为可计算的数字或去掉,以下是几种方法:

注:判断一个值是否为NaN,只能用isnan(),而不可用x==NaN;

i=find(~isnan(x));

x=x(i)

先找出值不是NaN的项的下标,将这些元素保留

x=x(find(~isnan(x)))同上,去掉NaN

x=x(~isnan(x));更快的做法

x(isnan(x))=[];消掉NaN

X(any(isnan(X)’),:)=[];把含有NaN的行都去掉

用此法可以从数据中去掉不相关的数据,看看下面的命令是干什么用的:

mu=mean(count);

sigma=std(count);

[n,p]=size(count)

outliers=abs(count—mu(ones(n,1),:))>3*sigma(ones(n,1),:);

nout=sum(outliers)

nout=

100

count(any(outliers'),:)=[];

点这看看答案

回归与曲线拟合

我们经常需要把观测到的数据表达为函数，假如有如下的对时间的观测：

t=[0.3.81.11.62.3]’;

y=[0.50.821.141.251.351.40]’;

plot(t,y,’o’),

gridon

多项式回归

由图可以看出应该可以用多项式来表达：y=a0+a1*t+a2*t^2

系数a0,a1,a2可以由最小平方拟合来确定，这一步可由反除号＂＼＂来完成

解下面的三元方程组可得：

X=[ones(size(t))tt.^2]

1.000000

1.00000.30000.0900

1.00000.80000.6400

1.00001.10001.2100

1.00001.60002.5600

1.00002.30005.2900

a=Xy

0.53180.9191–0.2387

a即为待求的系数，画图比较可得

T=(0:0.1:2.5)’;

Y=[ones(size(T))TT.^2]*a;

plot(T,Y,'–',t,y,'o',),gridon

结果令人失望，但我们可以增加阶数来提高精确度，但更明智的选择是用别的方法．

线性参数回归

形如：y=a0+a1*exp(-t)+a2*t*exp(-t)

计算方法同上：

X=[ones(size(t))exp(–t)t.*exp(–t)];

a=Xy

1.3974–0.89880.4097

T=(0:0.1:2.5)';

Y=[ones(size(T))exp(–T)(–T)]*a;

plot(T,Y,'–',t,y,'o'),gridon

看起来是不是好多了！

例子研究：曲线拟合

下面我们以美国人口普查的数据来研究一下有关曲线拟合的问题（MatLab是别人的,教学文档是别人

的,例子也是别人的，我只是一个翻译而已．．．）

loadcensus

这样我们得到了两个变量，cdate是1790至1990年的时间列向量(10年一次),pop是相应人口数列向量.

上一小节所讲的多项式拟合可以用函数polyfit()来完成,数字指明了阶数

p=polyfit(cdate,pop,4)

Warning:Matrixisclosetosingularorbadlyscaled.

=5.429790e–20

1.0e+05*

0.0000–0.00000.0000–0.01266.0020

产生警告的原因是计算中的cdata值太大，在计算中的Vandermonde行列式使变换产生了问题，解决的方

法之一是使数据标准化．

预处理：标准化数据

数据的标准化是对数据进行缩放，以使以后的计算能更加精确，一种方法是使之成为0均值：

sdate=(cdate–mean(cdate))./std(cdate)

现在再进行曲线拟合就没事了！

p=polyfit(sdate,pop,4)

0.70470.921023.470673.859862.2285

pop4=polyval(p,sdate);

plot(cdate,pop4,'–',cdate,pop,'+'),gridon

在上面的数据标准化中，也可以有别的算法，如令１７９０年的人口数为０．

余量分析

p1=polyfit(sdate,pop,1);

pop1=polyval(p1,sdate);

plot(cdate,pop1,'–',cdate,pop,'+')

res1=pop–pop1;

figure,plot(cdate,res1,'+')

p=polyfit(sdate,pop,2);

pop2=polyval(p,sdate);

plot(cdate,pop2,'–',cdate,pop,'+')

res2=pop–pop2;

figure,plot(cdate,res2,’+’)

p=polyfit(sdate,pop,4);

pop4=polyval(p,sdate);

plot(cdate,pop4,'–',cdate,pop,'+')

res4=pop–pop4;

figure,plot(cdate,res4,'+')

可以看出，多项式拟合即使提高了阶次也无法达到令人满意的结果

指数拟合

从人口增长图可以发现人数的增长基本是呈指数增加的，因此我们可以用年份的对数来进行拟合，这

儿，年数是标准化后的！

logp1=polyfit(sdate,log10(pop),1);

logpred1=10.^polyval(logp1,sdate);

semilogy(cdate,logpred1,'–',cdate,pop,'+');

gridon

logres2=log10(pop)–polyval(logp2,sdate);

plot(cdate,logres2,'+')

上面的图不令人满意，下面，我们用二阶的对数分析：

logp2=polyfit(sdate,log10(pop),2);

logpred2=10.^polyval(logp2,sdate);

semilogy(cdate,logpred2,'–',cdate,pop,'+');

gridon

r=pop–10.^(polyval(logp2,sdate));

plot(cdate,r,'+')

这种余量分析比多项式拟合的余量分析图案要随机的多（没有很强的规律性），可以预见，随着人数的增

加，余粮所反映的不确定度也在增加，但总的说来，这种拟合方式要强好多！

误差边界

误差边界常用来反映你所用的拟合方式是否适用于数据，为得到误差边界，只需在polyfit()中传递第二

个参数，并将其送入polyval().

下面是一个二阶多项式拟合模型,年份已被标准化,下面的代码用了2σ,对应于95%的可置信度:

[p2,S2]=polyfit(sdate,pop,2);

[pop2,del2]=polyval(p2,sdate,S2);

plot(cdate,pop,'+',cdate,pop2,'g–',cdate,pop2+2*del2,'r:',...

cdate,pop2–2*del2,'r:'),

gridon

差分方程和滤波

MatLab中的差分和滤波基本都是对向量而言的，向量则是存储取样信号或序列的．

函数y=filter(b,a,x)将用a,b描述的滤波器处理向量x,然后将其存储在向量y中,

filter()函数可看为是一差分方程a1y(n)=b1*x(1)+b2*x(2)+...-a2*y(2)-...

如有以下差分方程:y(n)=1/4*x(n)+1/4*x(n-1)+1/4*x(n-2)+1/4*x(n-3),则

a=1;

b=[1/41/41/41/4];

我们载入数据,取其第一列,并计算有:

x=count(:,1);

y=filter(b,a,x);

t=1:length(x);

plot(t,x,'–.',t,y,'–'),

gridon

legend('OriginalData','SmoothedData',2)

实现所表示的就是滤波后的数据,它代表了4小时的平均车流量

MatLab的信号处理工具箱中提供了很多用来滤波的函数,可用来处理实际问题!

快速傅立叶变换(FFT)

傅立叶变换能把信号按正弦展开成不同的频率值,对于取样信号,用的是离散傅立叶变换.

FFT是离散傅立叶变换的一种高速算法,在信号和图像处理中有极大的用处!

fft离散傅立叶变换

fft2二维离散傅立叶变换

fftnn维离散傅立叶变换

ifft离散傅立叶反变换

ifft2二维离散傅立叶反变换

ifftnn维离散傅立叶反变换

abs幅度

angle相角

unwrap相位按弧度展开,大于π的变换为2π的补角

fftshift把零队列移至功率谱中央

cplxpair把数据排成复数对

nextpow2下两个更高的功率

向量x的FFT可以这样求:

x=[437–91000]’

y=fft(x)

6.0000

11.4853–2.7574i

–2.0000–12.0000i

–5.4853+11.2426i

18.0000

–5.4853–11.2426i

–2.0000+12.0000i

11.4853+2.7574i

x虽然是实数,但y是复数,其中,第一个是因为它是常数相加的结果,第五个则对应于奈奎斯特频率,后三个

数是由于负频率的影响,它们是前面三个数的共轭值!

下面,让我们来验证一下太阳黑子活动周期是11年!Wolfer数记录了300年太阳黑子的数量及大小:

year=sunspot(:,1);

wolfer=sunspot(:,2);

plot(year,wolfer)

title(’SunspotData’)

现在来看看其FFT:

Y=fft(wolfer);

Y的幅度是功率谱,画出功率谱和频率的对应关系就得出了周期图,去掉第一点,因为他只是所有数据的和,

画图有:

N=length(Y);

Y(1)=[];

power=abs(Y(1:N/2)).^2;

nyquist=1/2;

freq=(1:N/2)/(N/2)*nyquist;

plot(freq,power),

gridon

xlabel(’cycles/year’)

title(’Periodogram’)

上面的图看起来不大方便,下面我们画出频谱-周期图

period=1./freq;

plot(period,power),

axis([04002e7]),

gridon

ylabel(’Power’)

xlabel(’Period(Years/Cycle)’)

为了得出精确一点的解,如下:

[mpindex]=max(power);

period(index)

ans=

11.0769

变换后的幅度和相位

abs()和angle()是用来计算幅度和相位的

先创建一信号,再进行分析,unwarp()把相位大于π的变换为2π的补角:

t=0:1/99:1;

x=sin(2*pi*15*t)+sin(2*pi*40*t);

y=fft(x);

m=abs(y);

p=unwrap(angle(y));

f=(0:length(y)–1)'*99/length(y);

subplot(2,1,1),plot(f,m),

ylabel('ude'),gridon

subplot(2,1,2),plot(f,p*180/pi)

ylabel('Phase[Degrees]'),gridon

xlabel('Frequency[Hertz]')

可以发现幅度曲线关于奈奎斯特频率对称,只有0-50Hz的信息是有用的!

FFT的长度与速度

可以为FFT加上第二个参数,告诉MatLab这是n点FFT.如y=fft(x,n),若x长度大于n,软件自动补0,

否则截取x.

若:

1.n为2的幂,软件将执行基2快速傅立叶算法,这时的运算速度是最快的

2.n为合数,软件将n分解为素数来算,计算量与n的值有关.n为1013将比1000点的速度慢的多!

3.n为素数,软件执行DFT的公式,此时最慢

矩阵和线性代数

MatLab中的矩阵

MatLab中有好多函数可以产生不同的矩阵,下面就让我们产生两个3*3的矩阵,这一章中,我们的学习就靠

她们了!!!

A=pascal(3)

111

123

136

B=magic(3)

816

357

492

还有一个3*2的随机矩阵:

C=fix(10*rand(3,2))

看看列矩阵，行矩阵，以及常数的表达：

u=[3;1;4]

v=[20—1]

s=7

产生的矩阵是：

20—1

加减法

X=A+B

927

4710

5128

Y=X–A

816

357

492

若二矩阵维数不统一，则会出错！

X=A+C

Errorusing==>+

Matrixdimensionsmustagree.

向量的乘积与转置

x=v*u

X=u*v

60—3

20—1

80—4

X=B'

834

159

672

如ｘ与ｙ均是列向量，则x*y无解,但下二表达式却可以:

x'*y

y'*x

称内积或点积.

下面的语句产生单位矩阵

eye(m,n)

若用eye(n)则产生n*n的方阵

解线性方程

情况一:

x=Au

—12

又如：

X=AB

19–3—1

—17413

60—6

情况二；y是不同时刻t时的观测值:

t=[0.3.81.11.62.3]';

y=[.82.72.63.60.55.50]';

若函数形式是:y(t)=c1+c2*exp(t);

构造矩阵:

E=[ones(size(t))exp(–t)]

1.00001.0000

1.00000.7408

1.00000.4493

1.00000.3329

1.00000.2019

1.00000.1003

则可求得系数c1及c2

c=Ey

0.47600.3413

表明：y(t)=0.4760+0.3413*exp(t)

画图如下:

T=(0:0.1:2.5)';

Y=[ones(size(T))exp(–T)]*c;

plot(T,Y,'–',t,y,'o')

转置与行列式

若A是方阵,且是非奇异的,则:

d=det(A)

X=inv(A)

3—31

—35—2

1—21

若c不是方阵,则用pinv:

X=pinv(C)

0.1159—0.07290.0171

—0.05340.11520.0418

那么我们可以发现，下面３个命令具有同样的功效（Ａ是ｍ*ｎ的矩阵，ｍ＞ｎ）：

x=Ab

x=pinv(A)*b

x=inv(A’*A)*A’*b

.及Cholesky分解

MatLab求解线性方程建立在以下三个分解之上:

Cholesky分解

Guass(高斯)分解

正交分解

Cholesky分解

A=p*p'

让我们临时把A变一变:

A=pascal(6)

111111

123456

136101521

1410203556

A是二项式系数,每一项是其左方与上方系数之和,求其Cholesky分解系数有:

R=chol(A)

111111

012345

0013610

0001410

000015

000001

Ｒ认识二项式系数．

这样对于线性方程便可化简：

A*x=b

R'*R*x=b

x=R(R'b)

复杂度由O(n^3)变为O(n^2);

LU分解

A=LU

其中,L时下三角阵,U是上三角阵,如:

[L,U]=lu(B)

1.000000

0.37500.54411.0000

0.50001.00000

8.00001.00006.0000

08.5000—1.0000

005.2941

同样：

A*x=b可以解为

x=U(Lb)

QR分解

正交阵有如下性质:

Q'Q=I

正交阵的好处在于,她保持了原阵的长度,角度,并且在计算的过程中不会扩大误差.

RQ分解如下:

A=QR或AP=QR

其中,Q是正交阵,R是上三角阵.

矩阵的幂与指数

若A是方阵,p是正数,则

X=A^2

3610

61425

102546

若Ａ是方阵，且是非奇异的，则X=A^(-P)将inv(A)P次方,如:

Y=B^(–3)

0.0053—0.00680.0018

—0.00340.00010.0036

—0.00160.0070—0.0051

分数词幂将由Ａ的特征值决定．

若是对矩阵的每个元素进行幂，用.^,如

X=A.^2

111

149

1936

sqrtm(A)计算A^(1/2)，但要更精确，而

sqrt(A)则计算A.^(1/2),是一个元素一个元素的算．

dx/dt=Ax,可以表示为x(t)=exp(tA)*x(0);

下面来看看如何计算:--expm(A)

0—6—1

62—16

—520—10

x0=

计算如下：

X=[];

fort=0:.01:1

X=[Xexpm(t*A)*x0];

end

作图有：

plot3(X(1,:),X(2,:),X(3,:),'–o')

特征值

Av=λv

若L是特阵值矩阵,则特征向量是V:

AV=VL;

如下:

0—6—1

62—16

—520—10

lambda=eig(A)

lambda=

—3.0710

—2.4645+17.6008i

—2.4645-17.6008i

由exp(λt)可以看出exp(At)(见上小节)

若用二参数调用函数eig(),则返回特征向量及特征值矩阵:

[V,D]=eig(A)

—0.8326—0.1203+0.2123i—0.1203–0.2123i

—0.35530.4691+0.4901i0.4691–0.4901i

—0.42480.6249–0.2997i0.6249+0.2997i

—3.071000

0—2.4645+17.6008i0

00—2.4645—17.6008i

对于下面的矩阵：

61219

—9—20—33

4915

0.47410.4082—0.4082

—0.8127—0.81650.8165

0.33860.4082—0.4082

—1.000000

01.00000

001.0000

可以看出，有二特征值是一样的，其特征向量仅差一个符号，在SymbolicMathToolbox中提供了Jordan

标准型的函数,如下:

[X,J]=jordan(A)

—1.75001.50002.7500

3.0000—3.0000—3.0000

—1.25001.50001.2500

—100

011

001

常微分方程

常微分方程(OdinaryDifferentialEquations---ODE)

类别函数描述

解常微分方程

ode45

ode23

ode113

ode15s

ode23s

ode23t

ode23tb

常微分方程选项

odeset

常微分方程输出选

项

odeplot

odephas2

odephas3

odeprint

如何表述问题

多项式与插值

多项式

多项式的表达

MatLab中用按降幂排列的多项式系数组成的行向量表示多项式,如:

p(x)=x^3-2x-5被表示为:

p=[10–2–5];

多项式的根

r=roots(p)

2.0946

–1.0473+1.1359i

–1.0473–1.1359i

根被储存为列向量．

若要由方程的根构造多项式，则

p2=poly(r)

p2=

18.8818e-16–2–5

多项式估计

可以用多项式估计出多项式在某一点的值：

polyval(p,5)

ans=

110

同样也可以估计矩阵多项式的值p(X)=X^3–2X–5I,

X=[245;–103;715];

Y=polyvalm(p,X)

377179439

11181136

490253639

卷积

多项式相乘是一个卷积的过程,conv()

a=[123];

b=[456];

c=conv(a,b)

413282718

多项式相除是其逆过程，用deconv():

[q,r]=deconv(c,a)

456

00000

多项式曲线逼近

polyfit(x,y,n)能用多项式逼近由x,y向量提供的数据,n是其阶数,如:

x=[12345];

y=[5.543.1128290.7498.4];

p=polyfit(x,y,3)

–0.191731.5821–60.326235.3400

将图画出

x2=1:.1:5;

y2=polyval(p,x2);

plot(x,y,’o’,x2,y2)

gridon

分式多项式分解

residue()可将分式多项式分解如下:

对于下式

分解为:

b=[–48];

a=[168];

[r,p,k]=residue(b,a)

–128

–4–2

[]

重载此函数可以完成分式多项式相加：

[b2,a2]=residue(r,p,k)

b2=

–48

a2=

168

插值

插值是在已知的数据列中，估计别点的函数值．

一维插值

一维插值在MatLab中有两种方法:

@多项式插值

@建立在FFT上的插值

多项式插值

yi=interp1(x,y,xi,method)

x是坐标向量,y是数据向量,xi是待估计点向量,method是插值方法,

method有四种:

t寻找最近数据点,由其得出函数值;

线性插值(该函数的默认方法);

样条插值,数据点处光滑--左导等于右导;

三次插值

以上四种方法得出的数据值一个比一个精确,而所需内存及计算时间也一个比一个要大要长.

建立在FFT上的插值

这种方法利用了快速傅立叶变换

y=interpft(x,n),其中,x含有周期性的函数值.

二维插值

ZI=interp2(X,Y,Z,XI,YI,method)

method有三种:

t寻找最近数据点,由其得出函数值;

二维线性插值

二维三次插值

下面来看看二维插值的例子:

先创造数据点:

[x,y]=meshgrid(–3:1:3);

z=peaks(x,y);

surf(x,y,z)

再比较一下不同的插值

[xi,yi]=meshgrid(–3:0.25:3);

zi1=interp2(x,y,z,xi,yi,'nearest');

zi2=interp2(x,y,z,xi,yi,'bilinear');

zi3=interp2(x,y,z,xi,yi,'bicubic');

数值分析

本节介绍的是关于函数的函数(functionfunctions),这些函数是用来处理函数而非数值的;

类别函数描述

绘图

优化

求解

fplot画出函数

fminbnd由一有范围限制的变量找出函数的最小值

fminsearch由几个变量找出函数的最小值

fzero找出函数的解(零值)

数值积分

quad低阶数值估计积分

quad8高阶数值估计积分

dblquad二重积分

数值微分见下一章

MatLab中的函数表达

MatLab中用M文件来表示函数,设有如下函数:

他别表示为一称为hump.m的文件中:

functiony=humps(x)

y=1./((x–0.3).^2+0.01)+1./((x–0.9).^2+0.04)–6;

这个函数文件可用于数值分析的函数中.

第二种方法就是创造一个行内对象(inline()),方法如下:

f=inline(‘1./((x–0.3).^2+0.01)+1./((x–0.9).^2+0.04)–6’);

用了上面的方法创造了函数文件,我们就可以找出函数在2的值:

f(2.0)

ans=

–4.8552

用创造行内对象的方法还可以创造多参数的函数,如下:

f=inline('y*sin(x)+x*cos(y)','x','y')

f(pi,2*pi)

ans=

3.1416

把函数画出来

fplot()可画出在给定范围内的函数值,如下

fplot('humps',[–55])

gridon

可通过限制y轴来放大图形

fplot('humps',[–55–1025])

gridon

你也可直接在fplot()中传递表达式,如:

fplot('2*sin(x+3)',[–11])

更可在一附图中画多个函数,如下

fplot('[2*sin(x+3),humps(x)]',[–11])

式中,[2*sin(x+3),humps(x)]组成了一个矩阵,每一列都是对应于x的函数

函数的最小值与解

找出一变量的函数的极值

x=fminbnd(’humps’,0.3,1)

0.6370

你可通过向fminbnd()函数传递一个函数optimset()作为参数来把此过程显示为列表形式:

x=fminbnd(’humps’,0.3,1,optimset(’Display’,’iter’))

Func-countxf(x)Procedure

10.56737612.9098initial

20.73262413.7746golden

30.46524825.1714golden

40.64441611.2693parabolic

50.641311.2583parabolic

60.63761811.2529parabolic

70.63698511.2528parabolic

80.63701911.2528parabolic

90.63705211.2528parabolic

0.6370

多变量函数极值

先创造一个ｍ文件，three_var.m:

functionb=three_var(v)

x=v(1);

y=v(2);

z=v(3);

b=x.^2+2.5*sin(y)–z^2*x^2*y^2;

现在,以x=–0.6,y=–1.2,z=0.135为起始点找出函数的极值:

v=[–0.6–1.20.135];

a=fminsearch('three_var',v)

0.0000–1.57080.1803

设置寻找极值的参数

x=fminbnd(fun,x1,x2,options)或

x=fminsearch(fun,x0,options)

其中，options是优化工具箱中(OptimizationToolbox)中的函数所用的一个结构,可如下设置

options=optimset('Display','iter');

y用来设置是否显示中间过程,如为:"iter"则显示,为"off"则不显示,为"final"则只显示

最后结果;

1X设置结果的误差范围,默认值是:1.e–4.

Eval设置函数运行次数的上限,默认fminbnd()是500次,fminsearch()是200*length(x0)

次

找出函数的解(零点值)

fzero()找出函数的零点值,你可以给出一个起始点,函数会从点开始搜索直到找到一个异号的值,最终给出

解;

如果你知道函数会于哪两点异号,你可以给出一个两点的向量,表明起始值和起始搜索步长

a=fzero('humps',–0.2)

–0.1316

验证一下，此函数值的确很接近０，

humps(a)

ans=

8.8818e–16

再看看下面的命令，看看你是否能看懂！

humps(1)

ans=

humps(–1)

ans=

–5.1378

options=optimset('Display','iter');

a=fzero('humps',[–11],options)

Func-countxf(x)Procedure

1–1–5.13779initial

1116initial

2–0.513876–4.02235interpolation

30.24306271.6382bisection

4–0.473635–3.83767interpolation

5–0.1152870.414441bisection

6–0.150214–0.423446interpolation

7–0.132562–0.0226907interpolation

8–0.131666–0.0011492interpolation

9–0.1316181.88371e–07interpolation

10–0.131618–2.7935e–11?interpolation

11–0.1316188.88178e–16interpolation

12–0.131618–9.76996e–15interpolation

–0.1316

最优化问题常常需要多次的运算才能汇集到一点，因此，起始点的选择显得至关重要，它不仅能提高

效率，更可使我们找到的值不是局部的极值，而是全局的最值；

下面是可能遇见的问题及解决的方法

得出的结果不是全局最值用不同的起始值或不同的起始步长去寻找

对应与x,无法计算出f改变你的函数使之含有罚函数,给f更大的数值空间

计算似乎进入了无限循环或

返回的值不是最小值

函数返回的是Inf,NaN,或复数值时,可在M文件中加上一些判断避免情况的

发生,如:isreal(),isfinite()等函数

数值积分

某区域内函数所围的区间可由数值积分来确定，MatLab中的一维函数积分用quad(),quad8().如:

q=quad(’humps’,0,1)

29.8583

函数可有第4个参数,指明误差范围,还可有第5个参数,用来画出图形.

例子:计算曲线的长度

如有一曲线如下确定:

x(t)=sin(2t);y(t)=cos(t);z(t)=t;

此函数可如下画出:

t=0:0.1:3*pi;

plot3(sin(2*t),cos(t),t)

hcurce()用来计算曲线长度:

functionf=hcurve(t)

f=sqrt(4*cos(2*t).^2+sin(t).^2+1);

下面就是计算结果:

len=quad(’hcurve’,0,3*pi)

len=

1.7222e+01

二重积分

result=dblquad(’integrnd’,xmin,xmax,ymin,ymax);

默认状态下,dblquad是用的quad()来计算,为更精确,可

result=dblquad(’integrnd’,xmin,xmax,ymin,ymax,[],’quad8’);

三维及多维插值

列出函数,其余从略

VI=interp3(X,Y,Z,V,XI,YI,ZI,method)

VI=interpn(X1,X2,X3...,V,Y1,Y2,Y3,...,method)

此文出于:/非本人所写!

👁️ 阅读量：0

🔖 本文标签：

🔥 最新发布文章