第5讲 MATLAB数值计算.docx
《第5讲 MATLAB数值计算.docx》由会员分享，可在线阅读，更多相关《第5讲 MATLAB数值计算.docx（16页珍藏版）》请在冰豆网上搜索。
第5讲MATLAB数值计算
第五讲MATLAB数值计算
数据处理与多项式计算
数值微积分
离散傅立叶变换
线性方程组求解
非线性方程与最优化问题求解
常微分方程的数值求解
一、数据处理与多项式计算
（一）数据统计与分析
1.求矩阵最大元素和最小元素
MATLAB提供的求数据序列的最大值和最小值的函数分别为max和min，两个函数的调用格式和操作过程类似。
（1）求向量的最大值和最小值
y=max（X）：
返回向量X的最大值存入y，如果X中包含复数元素，则按模取最大值。
[y,I]=max（X）：
返回向量X的最大值存入y，最大值的序号存入I，如果X中包含复数元素，则按模取最大值。
求向量X的最小值的函数是min（X），用法和max（X）完全相同。
例1求向量x的最大值。
命令如下：
x=[-43,72,9,16,23,47];
y=max（x）%求向量x中的最大值
[y,l]=max（x）%求向量x中的最大值及其该元素的位置
（2）求矩阵的最大值和最小值
求矩阵A的最大值的函数有3种调用格式，分别是：
max（A）：
返回一个行向量，向量的第i个元素是矩阵A的第i列上的最大值。
[Y,U]=max（A）：
返回行向量Y和U，Y向量记录A的每列的最大值，U向量记录每列最大值的行号。
max（A,[],dim）：
dim取1或2。
dim取1时，该函数和max（A）完全相同；dim取2时，该函数返回一个列向量，其第i个元素是A矩阵的第i行上的最大值。
求最小值的函数是min，其用法和max完全相同。
例2分别求矩阵A中各列和各行元素中的最大值，并求整个矩阵的最大值和最小值。
A=[17,0,1,0,15;23,5,7,14,16;4,0,13,0,22;10,12,19,21,3;...
11,18,25,2,19];
（3）两个向量或矩阵对应元素的比较
函数max和min还能对两个同型的向量或矩阵进行比较，调用格式为：
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）;
Y=sin（X）;
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（'ex',0,1）
I=
0.74682413398845
也可不建立关于被积函数的函数文件，而使用语句函数（内联函数）求解，命令如下：
g=inline（'exp（-x.^2）'）;
%定义一个语句函数g（x）=exp（-x^2）
I=quadl（g,0,1）%注意函数名不加'号
I=
0.74682413398845
formatshort
（2）被积函数由一个表格定义
在科学实验和工程应用中，函数关系往往是不知道的，只有实验测定的一组样本点和样本值，这时，就无法使用quad函数计算其定积分。
在MATLAB中，对由表格形式定义的函数关系的求定积分问题用trapz（X,Y）函数。
其中向量X、Y定义函数关系Y=f（X）。
X、Y是两个等长的向量：
X=（x1,x2,…,xn），Y=（y1,y2,…,yn），并且x1例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函数求解。
globalki;ki=0;
I=dblquad（'fxy',-2,2,-1,1）
ki
I=
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作为迭代初值
ans=
-5.1926
fzero（'fz',1）%以1作为迭代初值
ans=
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=x-1/x+5;
上述函数文件也可用一个语句函数代替：
f=inline（'x-1/x+5'）
再在MATLAB命令窗口，输入命令：
fminbnd（'f',-10,-1）
%求函数在（-10，-1）内的最小值点和最小值
fminbnd（'f',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（'fname',tspan,y0）
[t,y]=ode45（'fname',tspan,y0）
其中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，求微分方程的解：
t0=0;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
升级会员 