第5讲 MATLAB数值计算

1、第五讲 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)完全相同

2、。例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

3、，其用法和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函数的用法和ma

4、x完全相同。例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

5、为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)

6、：当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时，该函数等同于cums

7、um(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时，求各列元素的标准

8、方差；当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 生成满足正态

9、分布的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)使其通过或

10、近似通过有限序列的资料点(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=lins

11、pace(0,2*pi,50);Y=sin(X);plot(X,Y)P=polyfit(X,Y,3) %得到3次多项式的系数和误差%观察拟合的好坏k=polyval(P,X)hold onplot(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的逆函数，即有P

12、1=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函数用来求代数多项式的值，其调用格式为

13、：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.*

14、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函数建

15、立起该多项式，其调用格式为：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=di

16、ff(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的13阶差分。命令如下：X=linspace(0,2*pi,10);Y=sin(X);DY=diff(Y); %计算Y的一阶差分D2Y=diff(Y,2); %计算Y的二阶差分，也可用 命令diff(DY)计算D3Y=diff(

17、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,to

18、l,trace) quadl(filename,a,b,tol,trace)例5.20 用两种不同的方法求定积分。先建立一个函数文件ex.m：function ex=ex(x)ex=exp(-x.2);然后在MATLAB命令窗口，输入命令：format longI=quad('ex',0,1) %注意函数名应加字符引号I = 0.74682418072642I=quadl('ex',0,1)I = 0.74682413398845也可不建立关于被积函数的函数文件，而使用语句函数(内联函数)求解，命令如下：g=inline('exp(-x.2)')

19、; %定义一个语句函数g(x)=exp(-x2)I=quadl(g,0,1) %注意函数名不加'号I = 0.74682413398845format short(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

20、用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：function f=fxy(x,y)global ki;ki=ki+1; %ki用

21、于统计被积函数的调用次数f=exp(-x.2/2).*sin(x.2+y);(2) 调用dblquad函数求解。global ki;ki=0;I=dblquad('fxy',-2,2,-1,1)kiI = 1.57449318974494ki = 1038 四、线性方程组求解（一）直接解法1利用左除运算符的直接解法对于线性方程组Ax=b，可以利用左除运算符“”求解： x=Ab例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=Ab（二）用矩阵求逆方法求解线性方程组在

22、线性方程组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最近的那个根。t

23、ol控制结果的相对精度，缺省时取tol=eps，trace指定迭代信息是否在运算中显示，为1时显示，为0时不显示，缺省时取trace=0。例5.33 求f(x)在x0=-5和x0=1作为迭代初值时的零点。先建立函数文件fz.m：function f=fz(x)f=x-1/x+5;然后调用fzero函数求根。：fzero('fz',-5) %以-5作为迭代初值ans = -5.1926fzero('fz',1) %以1作为迭代初值ans = 0.1926 （二）无约束最优化问题求解在实际应用中，许多科学研究和工程计算问题都可以归结为一个最小化问题，如能量最小、时间

24、最短等。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)的最大值，所

25、以fminbnd(-f,x1,x2)返回函数f(x)在区间(x1,x2)上的最大值。例5.36 求函数在区间(-10，1)和(1，10)上的最小值点。首先建立函数文件f.m：function f=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：function f=fxyz(u)x=u(1);y=u(2);z=u(3);f=x+y.2./x/4+z.2./y+2./z;在MALAB命令窗口，输入