第5讲 MATLAB数值计算Word文件下载.docx

第5讲 MATLAB数值计算Word文件下载.docx

《第5讲 MATLAB数值计算Word文件下载.docx》由会员分享，可在线阅读，更多相关《第5讲 MATLAB数值计算Word文件下载.docx（13页珍藏版）》请在冰豆网上搜索。

U=max（A,B）：

A,B是两个同型的向量或矩阵，结果U是与A,B同型的向量或矩阵，U的每个元素等于A,B对应元素的较大者。

U=max（A,n）：

n是一个标量，结果U是与A同型的向量或矩阵，U的每个元素等于A对应元素和n中的较大者。

min函数的用法和max完全相同。

例3：

求两个2×

3矩阵x,y所有同一位置上的较大元素构成的新矩阵p。

2.求矩阵的平均值和中值

求数据序列平均值的函数是mean，求数据序列中值的函数是median。

两个函数的调用格式为：

mean（X）：

返回向量X的算术平均值。

median（X）：

返回向量X的中值。

mean（A）：

返回一个行向量，其第i个元素是A的第i列的算术平均值。

median（A）：

返回一个行向量，其第i个元素是A的第i列的中值。

mean（A,dim）：

当dim为1时，该函数等同于mean（A）；

当dim为2时，返回一个列向量，其第i个元素是A的第i行的算术平均值。

median（A,dim）：

当dim为1时，该函数等同于median（A）；

当dim为2时，返回一个列向量，其第i个元素是A的第i行的中值。

3.矩阵元素求和与求积

数据序列求和与求积的函数是sum和prod，其使用方法类似。

设X是一个向量，A是一个矩阵，函数的调用格式为：

sum（X）：

返回向量X各元素的和。

prod（X）：

返回向量X各元素的乘积。

sum（A）：

返回一个行向量，其第i个元素是A的第i列的元素和。

prod（A）：

返回一个行向量，其第i个元素是A的第i列的元素乘积。

sum（A,dim）：

当dim为1时，该函数等同于sum（A）；

当dim为2时，返回一个列向量，其第i个元素是A的第i行的各元素之和。

prod（A,dim）：

当dim为1时，该函数等同于prod（A）；

当dim为2时，返回一个列向量，其第i个元素是A的第i行的各元素乘积。

例4求矩阵A的每行元素的乘积和全部元素的乘积。

4.矩阵元素累加和与累乘积

在MATLAB中，使用cumsum和cumprod函数能方便地求得向量和矩阵元素的累加和与累乘积向量，函数的调用格式为：

cumsum（X）：

返回向量X累加和向量。

cumprod（X）：

返回向量X累乘积向量。

cumsum（A）：

返回一个矩阵，其第i列是A的第i列的累加和向量。

cumprod（A）：

返回一个矩阵，其第i列是A的第i列的累乘积向量。

cumsum（A,dim）：

当dim为1时，该函数等同于cumsum（A）；

当dim为2时，返回一个矩阵，其第i行是A的第i行的累加和向量。

cumprod（A,dim）：

当dim为1时，该函数等同于cumprod（A）；

当dim为2时，返回一个向量，其第i行是A的第i行的累乘积向量。

5．求标准方差

（方差=每个元素与平均值的差的平方的平均值，标准差为方差的平方根）。

在MATLAB中，提供了计算数据序列的标准方差的函数std。

对于向量X，std（X）返回一个标准方差。

对于矩阵A，std（A）返回一个行向量，它的各个元素便是矩阵A各列或各行的标准方差。

std函数的一般调用格式为：

Y=std（A,flag,dim）

其中dim取1或2。

当dim=1时，求各列元素的标准方差；

当dim=2时，则求各行元素的标准方差。

flag取0或1，当flag=0时，按σ1所列公式计算标准方差，当flag=1时，按σ2所列公式计算标准方差。

缺省flag=0，dim=1。

例5对二维矩阵A，从不同维方向求出其标准方差。

6．相关系数

MATLAB提供了corrcoef函数，可以求出数据的相关系数矩阵。

corrcoef函数的调用格式为：

corrcoef（X）：

返回从矩阵X形成的一个相关系数矩阵。

此相关系数矩阵的大小与矩阵X一样。

它把矩阵X的每列作为一个变量，然后求它们的相关系数。

corrcoef（X,Y）：

在这里，X,Y是向量，它们与corrcoef（[X,Y]）的作用一样。

例6生成满足正态分布的10000×

5随机矩阵，然后求各列元素的均值和标准方差，再求这5列随机数据的相关系数矩阵。

X=randn（10000,5）;

M=mean（X）

D=std（X）

R=corrcoef（X）

7.排序

MATLAB中对向量X是排序函数是sort（X），函数返回一个对X中的元素按升序排列的新向量。

sort函数也可以对矩阵A的各列或各行重新排序，其调用格式为：

[Y,I]=sort（A,dim）

其中dim指明对A的列还是行进行排序。

若dim=1，则按列排；

若dim=2，则按行排。

Y是排序后的矩阵，而I记录Y中的元素在A中位置。

（三）曲线拟合

（定义：

推求一个解析函数y=f（x）使其通过或近似通过有限序列的资料点（xi，yi），通常用多项式函数通过最小二乘法求得此拟合函数。

）

在MATLAB中，用polyfit函数来求得最小二乘拟合多项式的系数，再用polyval函数按所得的多项式计算所给出的点上的函数近似值。

polyfit函数的调用格式为：

[P,S]=polyfit（X,Y,m）

函数根据采样点X和采样点函数值Y，产生一个m次多项式P及其在采样点的误差向量S。

其中X,Y是两个等长的向量，P是一个长度为m+1的向量，P的元素为多项式系数。

polyval函数的功能是按多项式的系数计算x点多项式的值。

例5.11用一个3次多项式在区间[0，2π]内逼近函数y=sin（X）。

X=linspace（0,2*pi,50）;

Y=sin（X）;

plot（X,Y）

P=polyfit（X,Y,3）

%得到3次多项式的系数和误差

%观察拟合的好坏

k=polyval（P,X）

holdon

plot（X,k）

（四）多项式计算

1.多项式的四则运算

（1）多项式的加减运算

（2）多项式乘法运算

函数conv（P1,P2）用于求多项式P1和P2的乘积。

这里，P1、P2是两个多项式系数向量。

（3）多项式除法

函数[Q,r]=deconv（P1,P2）用于对多项式P1和P2作除法运算。

其中Q返回多项式P1除以P2的商式，r返回P1除以P2的余式。

这里，Q和r仍是多项式系数向量。

deconv是conv的逆函数，即有P1=conv（P2,Q）+r。

2.多项式的导函数

对多项式求导数的函数是：

p=polyder（P）：

求多项式P的导函数

p=polyder（P,Q）：

求P·

Q的导函数

[p,q]=polyder（P,Q）：

求P/Q的导函数，导函数的分子存入p，分母存入q。

上述函数中，参数P,Q是多项式的向量表示，结果p,q也是多项式的向量表示。

3.多项式求值

MATLAB提供了两种求多项式值的函数：

polyval与polyvalm，它们的输入参数均为多项式系数向量P和自变量x。

两者的区别在于前者是代数多项式求值，而后者是矩阵多项式求值。

（1）代数多项式求值

polyval函数用来求代数多项式的值，其调用格式为：

Y=polyval（P,x）

若x为一数值，则求多项式在该点的值；

若x为向量或矩阵，则对向量或矩阵中的每个元素求其多项式的值。

例5.14已知多项式x4+8x3-10，分别取x=1.2和一个2×

3矩阵为自变量计算该多项式的值。

（2）矩阵多项式求值

polyvalm函数用来求矩阵多项式的值，其调用格式与polyval相同，但含义不同。

polyvalm函数要求x为方阵，它以方阵为自变量求多项式的值。

设A为方阵，P代表多项式x3-5x2+8，那么polyvalm（P,A）的含义是：

A*A*A-5*A*A+8*eye（size（A））

而polyval（P,A）的含义是：

A.*A.*A-5*A.*A+8*ones（size（A））

例5.15仍以多项式x4+8x3-10为例，取一个2×

2矩阵为自变量分别用polyval和polyvalm计算该多项式的值。

4.多项式求根

预备知识：

n次多项式具有n个根，当然这些根可能是实根，也可能含有若干对共轭复根。

MATLAB提供的roots函数用于求多项式的全部根，其调用格式为：

x=roots（P），其中P为多项式的系数向量，求得的根赋给向量x，即x

（1）,x

（2）,…,x（n）分别代表多项式的n个根。

例5.16求多项式x4+8x3-10的根。

A=[1,8,0,0,-10];

x=roots（A）

若已知多项式的全部根，则可以用poly函数建立起该多项式，其调用格式为：

P=poly（x）

若x为具有n个元素的向量，则poly（x）建立以x为其根的多项式，且将该多项式的系数赋给向量P。

例5.17已知f（x）=3x5+4x3-5x2-7.2x+5

（1）计算f（x）=0的全部根。

（2）由方程f（x）=0的根构造一个多项式g（x），并与f（x）进行对比。

P=[3,0,4,-5,-7.2,5];

X=roots（P）%求方程f（x）=0的根

G=poly（X）%求多项式g（x）

二、数值微积分

（一）数值微分

2.数值微分的实现

在MATLAB中，没有直接提供求数值导数的函数，只有计算向前差分的函数diff，其调用格式为：

DX=diff（X）：

计算向量X的向前差分，DX（i）=X（i+1）-X（i），i=1,2,…,n-1。

DX=diff（X,n）：

计算X的n阶向前差分。

例如，diff（X,2）=diff（diff（X））。

DX=diff（A,n,dim）：

计算矩阵A的n阶差分，dim=1时（缺省状态），按列计算差分；

dim=2，按行计算差分。

例5.18设x由[0,2π]间均匀分布的10个点组成，求sinx的1~3阶差分。

X=linspace（0,2*pi,10）;

DY=diff（Y）;

%计算Y的一阶差分

D2Y=diff（Y,2）;

%计算Y的二阶差分，也可用命令diff（DY）计算

D3Y=diff（Y,3）;

%计算Y的三阶差分，也可用diff（D2Y）或diff（DY,2）

（二）数值积分

1.数值积分基本原理

求解定积分的数值方法多种多样，如简单的梯形法、辛普生（Simpson）法、牛顿－柯特斯（Newton-Cotes）法等都是经常采用的方法。

它们的基本思想都是将整个积分区间[a,b]分成n个子区间[xi,xi+1]，i=1,2,…,n，其中x1=a，xn+1=b。

这样求定积分问题就分解为求和问题。

积分计算一个函数与自变量轴之间的面积。

2.数值积分的实现

（1）被积函数是一个解析式

MATLAB提供了quad函数和quadl函数来求定积分。

它们的调用格式为：

quad（filename,a,b,tol,trace）

quadl（filename,a,b,tol,trace）

例5.20用两种不同的方法求定积分。

先建立一个函数文件ex.m：

functionex=ex（x）

ex=exp（-x.^2）;

然后在MATLAB命令窗口，输入命令：

formatlong

I=quad（'

ex'

0,1）%注意函数名应加字符引号

I=

0.74682418072642

I=quadl（'

0,1）

0.74682413398845

也可不建立关于被积函数的函数文件，而使用语句函数（内联函数）求解，命令如下：

g=inline（'

exp（-x.^2）'

）;

%定义一个语句函数g（x）=exp（-x^2）

I=quadl（g,0,1）%注意函数名不加'

号

formatshort

（2）被积函数由一个表格定义

在科学实验和工程应用中，函数关系往往是不知道的，只有实验测定的一组样本点和样本值，这时，就无法使用quad函数计算其定积分。

在MATLAB中，对由表格形式定义的函数关系的求定积分问题用trapz（X,Y）函数。

其中向量X、Y定义函数关系Y=f（X）。

X、Y是两个等长的向量：

X=（x1,x2,…,xn），Y=（y1,y2,…,yn），并且x1<

x2<

…<

xn，积分区间是[x1，xn]。

例5.21用trapz函数计算定积分。

在MATLAB命令窗口，输入命令：

X=0:

0.01:

1;

Y=exp（-X.^2）;

trapz（X,Y）

ans=

0.7468

（3）二重积分数值求解

使用MATLAB提供的dblquad函数就可以直接求出上述二重定积分的数值解。

该函数的调用格式为：

I=dblquad（f,a,b,c,d,tol,trace）

该函数求f（x,y）在[a,b]×

[c,d]区域上的二重定积分。

参数tol，trace的用法与函数quad完全相同。

例5.22计算二重定积分。

（1）建立一个函数文件fxy.m：

functionf=fxy（x,y）

globalki;

ki=ki+1;

%ki用于统计被积函数的调用次数

f=exp（-x.^2/2）.*sin（x.^2+y）;

（2）调用dblquad函数求解。

ki=0;

I=dblquad（'

fxy'

-2,2,-1,1）

ki

1.57449318974494

ki=

1038

四、线性方程组求解

（一）直接解法

1．利用左除运算符的直接解法

对于线性方程组Ax=b，可以利用左除运算符“\”求解：

x=A\b

例5.24用直接解法求解下列线性方程组。

A=[2,1,-5,1;

1,-5,0,7;

0,2,1,-1;

1,6,-1,-4];

b=[13,-9,6,0]'

;

x=A\b

（二）用矩阵求逆方法求解线性方程组

在线性方程组Ax=b两边各左乘A-1，有

A-1Ax=A-1b

由于A-1A=I，故得

x=A-1b

例用求逆矩阵的方法解线性方程组。

A=[1,2,3;

1,4,9;

1,8,27];

b=[5,-2,6]'

x=inv（A）*b　

五、非线性方程与最优化问题求解

（一）非线性方程数值求解

1.单变量非线性方程求解

在MATLAB中提供了一个fzero函数，可以用来求单变量非线性方程的根。

z=fzero（'

fname'

x0,tol,trace）

其中fname是待求根的函数文件名，x0为搜索的起点。

一个函数可能有多个根，但fzero函数只给出离x0最近的那个根。

tol控制结果的相对精度，缺省时取tol=eps，trace指定迭代信息是否在运算中显示，为1时显示，为0时不显示，缺省时取trace=0。

例5.33求f（x）在x0=-5和x0=1作为迭代初值时的零点。

先建立函数文件fz.m：

functionf=fz（x）

f=x-1/x+5;

然后调用fzero函数求根。

：

fzero（'

fz'

-5）%以-5作为迭代初值

-5.1926

1）%以1作为迭代初值

0.1926

（二）无约束最优化问题求解

在实际应用中，许多科学研究和工程计算问题都可以归结为一个最小化问题，如能量最小、时间最短等。

MATLAB提供了3个求最小值的函数，它们的调用格式为：

（1）[x,fval]=fminbnd（filename,x1,x2,option）：

求一元函数在（xl,x2）区间中的极小值点x和最小值fval。

（2）[x,fval]=fminsearch（filename,x0,option）：

基于单纯形算法求多元函数的极小值点x和最小值fval。

（3）[x,fval]=fminunc（filename,x0,option）：

基于拟牛顿法求多元函数的极小值点x和最小值fval。

MATLAB没有专门提供求函数最大值的函数，但只要注意到-f（x）在区间（a,b）上的最小值就是f（x）在（a,b）的最大值，所以fminbnd（-f,x1,x2）返回函数f（x）在区间（x1,x2）上的最大值。

例5.36求函数在区间（-10，1）和（1，10）上的最小值点。

首先建立函数文件f.m：

functionf=f（x）

上述函数文件也可用一个语句函数代替：

f=inline（'

x-1/x+5'

再在MATLAB命令窗口，输入命令：

fminbnd（'

f'

-10,-1）

%求函数在（-10，-1）内的最小值点和最小值

1,10）

%求函数在（1，10）内的最小值点。

例5.37求函数f在（0.5,0.5,0.5）附近的最小值。

建立函数文件fxyz.m：

functionf=fxyz（u）

x=u

（1）;

y=u

（2）;

z=u（3）;

f=x+y.^2./x/4+z.^2./y+2./z;

在MALAB命令窗口，输入命令：

[U,fmin]=fminsearch（'

fxyz'

[0.5,0.5,0.5]）%求函数的最小值点和最小值

六、常微分方程的数值求解

（一）龙格－库塔法简介

1.龙格－库塔法的实现

基于龙格－库塔法，MATLAB提供了求常微分方程数值解的函数，一般调用格式为：

[t,y]=ode23（'

tspan,y0）

[t,y]=ode45（'

其中fname是定义f（t,y）的函数文件名，该函数文件必须返回一个列向量。

tspan形式为[t0,tf],表示求解区间。

y0是初始状态列向量。

t和y分别给出时间向量和相应的状态向量。

例5.39设有初值问题，试求其数值解，并与精确解相比较（精确解为y（t）=）。

（1）建立函数文件funt.m。

functionyp=funt（t,y）

yp=（y^2-t-2）/4/（t+1）;

（2）求解微分方程。

t0=0;

tf=10;

y0=2;

[t,y]=ode23（'

funt'

[t0,tf],y0）;

%求数值解

y1=sqrt（t+1）+1;

%求精确解

t'

y'

y1'

y为数值解，y1为精确值，显然两者近似。

例5.40已知一个二阶线性系统的微分方程，绘制系统的时间响应曲线和相平面图。

函数ode23和ode45是对一阶常微分方程组设计的，因此对高阶常微分方程，需先将它转化为一阶常微分方程组，即状态方程。

建立一个函数文件sys.m：

functionxdot=sys（t,x）

xdot=[-2*x

（2）;

x

（1）];

取t0=0，tf=20，求微分方程的解：

tf=20;

[t,x]=ode45（'

sys'

[t0,tf],[1,0]）;

[t,x]

subplot（1,2,1）;

plot（t,x（:

2））;

%解的曲线，即t-x

subplot（1,2,2）;

plot（x（:

2）,x（:

1））%相平面曲线，即x-x’

axisequal