第6章MATLAB数值计算61-数据处理与多项式课件

单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,第6章 MATLAB数值计算,6.1 数据处理与多项式计算,6.2 数值微积分,6.3 离散傅立叶变换,6.4 线性方程组求解,6.5 非线性方程与最优化问题求解,6.6 常微分方程的数值求解,6.7 稀疏矩阵,第6章 MATLAB数值计算,1,6.1 数据处理与多项式计算,6.1.1 数据统计与分析,1. 求矩阵最大元素和最小元素,MATLAB提供的求数据序列的最大值和最小值的函数分别为max和min，两个函数的调用格式和操作过程类似。,（1）求向量的最大值和最小值,y=max(X)：返回向量X的最大值存入y，如果X中包含复数元素，则按模取最大值。,6.1 数据处理与多项式计算,2,y,I=max(X)：返回向量X的最大值存入y，最大值的序号存入I，如果X中包含复数元素，则按模取最大值。,求向量X的最小值的函数是min(X)，用法和max(X)完全相同。,例 求向量x的最大值。,命令如下：,x=-43,72,9,16,23,47;,y=max(x) %求向量x中的最大值,y,l=max(x) %求向量x中的最大值及其该元素的位置,y,I=max(X)：返回向量X的最大值存入y，,3,（2）求矩阵的最大值和最小值,求矩阵A的最大值的函数有3种调用格式，分别是：,max(A)：返回一个行向量，向量的第i个元素是矩阵A的第i列上的最大值。,Y,U=max(A)：返回行向量Y和U，Y向量记录A的每列的最大值，U向量记录每列最大值的行号。,（2）求矩阵的最大值和最小值,4,max(A,dim)：dim取1或2。dim取1时，该函数和max(A)完全相同；dim取2时，该函数返回一个列向量，其第i个元素是A矩阵的第i行上的最大值。,求最小值的函数是min，其用法和max完全相同。,例6.1 分别矩阵A中各列和各行元素中的最大值，并求整个矩阵的最大值和最小值。,max(A,dim)：dim取1或2。dim取,5,（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完全相同。,例 求两个23矩阵x, y所有同一位置上的较大元素构成的新矩阵p。,（3）两个向量或矩阵对应元素的比较,6,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行的中值。,2. 求矩阵的平均值和中值,7,3. 矩阵元素求和与求积,数据序列求和与求积的函数是sum和prod，其使用方法类似。设X是一个向量，A是一个矩阵，函数的调用格式为：,sum(X)：返回向量X各元素的和。,prod(X)：返回向量X各元素的乘积。,sum(A)：返回一个行向量，其第i个元素是A的第i列的元素和。,3. 矩阵元素求和与求积,8,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行的各元素乘积。,例6.2 求矩阵A的每行元素的乘积和全部元素的乘积。,prod(A)：返回一个行向量，其第i个元素是A的第i列的元,9,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行的累乘积向量。,4. 矩阵元素累加和与累乘积,10,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。,例6.4 对二维矩阵x，从不同维方向求出其标准方差。,5求标准方差,11,6相关系数,MATLAB提供了corrcoef函数，可以求出数据的相关系数矩阵。corrcoef函数的调用格式为：,corrcoef(X)：返回从矩阵X形成的一个相关系数矩阵。此相关系数矩阵的大小与矩阵X一样。它把矩阵X的每列作为一个变量，然后求它们的相关系数。,corrcoef(X,Y)：在这里，X,Y是向量，它们与corrcoef(X,Y)的作用一样。,6相关系数,12,例6.5 生成满足正态分布的100005随机矩阵，然后求各列元素的均值和标准方差，再求这5列随机数据的相关系数矩阵。,命令如下：,X=randn(10000,5);,M=mean(X),D=std(X),R=corrcoef(X),例6.5 生成满足正态分布的100005随机矩阵，然后求,13,7. 排序,MATLAB中对向量X是排序函数是sort(X)，函数返回一个对X中的元素按升序排列的新向量。,sort函数也可以对矩阵A的各列或各行重新排序，其调用格式为：,Y,I=sort(A,dim),其中dim指明对A的列还是行进行排序。若dim=1，则按列排；若dim=2，则按行排。Y是排序后的矩阵，而I记录Y中的元素在A中位置。,7. 排序,14,6.1.2 数据插值,1. 一维数据插值,在MATLAB中，实现这些插值的函数是interp1，其调用格式为：,Y1=interp1(X,Y,X1,method),函数根据X,Y的值，计算函数在X1处的值。X,Y是两个等长的已知向量，分别描述采样点和样本值，X1是一个向量或标量，描述欲插值的点，Y1是一个与X1等长的插值结果。method是插值方法，允许的取值有linear、nearest、cubic、spline。,6.1.2 数据插值,15,注意：X1的取值范围不能超出X的给定范围，否则，会给出“NaN”错误。,例6.7 给出概率积分的数据表如表6.1所示，用不同的插值方法计算f(0.472)。,例6.8 某检测参数f随时间t的采样结果如表5.1，用数据插值法计算t=2，7，12，17，22，17，32，37，42，47，52，57时的f值。,注意：X1的取值范围不能超出X的给定范围，否则，会给出“Na,16,2. 二维数据插值,在MATLAB中，提供了解决二维插值问题的函数interp2，其调用格式为：,Z1=interp2(X,Y,Z,X1,Y1,method),其中X,Y是两个向量，分别描述两个参数的采样点，Z是与参数采样点对应的函数值，X1,Y1是两个向量或标量，描述欲插值的点。Z1是根据相应的插值方法得到的插值结果。 method的取值与一维插值函数相同。X,Y,Z也可以是矩阵形式。,同样，X1,Y1的取值范围不能超出X,Y的给定范围，否则，会给出“NaN”错误。,2. 二维数据插值,17,例6.9 设z=x,2,+y,2,，对z函数在0,10,2区域内进行插值。,例6.10 某实验对一根长10米的钢轨进行热源的温度传播测试。用x表示测量点0:2.5:10(米)，用h表示测量时间0:30:60(秒)，用T表示测试所得各点的温度()。试用线性插值求出在一分钟内每隔10秒、钢轨每隔0.5米处的温度。,例6.9 设z=x2+y2，对z函数在0,10,2,18,6.1.3 曲线拟合,在MATLAB中，用polyfit函数来求得最小二乘拟合多项式的系数，再用polyval函数按所得的多项式计算所给出的点上的函数近似值。,polyfit函数的调用格式为：,P,S=polyfit(X,Y,m),函数根据采样点X和采样点函数值Y，产生一个m次多项式P及其在采样点的误差向量S。其中X,Y是两个等长的向量，P是一个长度为m+1的向量，P的元素为多项式系数。,polyval函数的功能是按多项式的系数计算x点多项式的值。,6.1.3 曲线拟合,19,例6.11 用一个3次多项式在区间0，2内逼近函数。,命令如下,：,X=linspace(0,2*pi,50);,Y=sin(X);,P=polyfit(X,Y,3) %得到3次多项式的系数和误差,例6.11 用一个3次多项式在区间0，2内逼近函数。,20,6.1.4 多项式计算,1. 多项式的四则运算,（1）多项式的加减运算,（2）多项式乘法运算,函数conv(P1,P2)用于求多项式P1和P2的乘积。这里，P1、P2是两个多项式系数向量。,6.1.4 多项式计算,21,（3）多项式除法,函数Q,r=deconv(P1,P2)用于对多项式P1和P2作除法运算。其中Q返回多项式P1除以P2的商式，r返回P1除以P2的余式。这里，Q和r仍是多项式系数向量。,deconv是conv的逆函数，即有P1=conv(P2,Q)+r。,（3）多项式除法,22,2. 多项式的导函数,对多项式求导数的函数是：,p=polyder(P)：求多项式P的导函数,p=polyder(P,Q)：求PQ的导函数,p,q=polyder(P,Q)：求P/Q的导函数，导函数的分子存入p，分母存入q。,上述函数中，参数P,Q是多项式的向量表示，结果p,q也是多项式的向量表示。,2. 多项式的导函数,23,3. 多项式求值,MATLAB提供了两种求多项式值的函数：polyval与polyvalm，它们的输入参数均为多项式系数向量P和自变量x。两者的区别在于前者是代数多项式求值，而后者是矩阵多项式求值。,3. 多项式求值,24,(1)代数多项式求值,polyval函数用来求代数多项式的值，其调用格式为：,Y=polyval(P,x),若x为一数值，则求多项式在该点的值；若x为向量或矩阵，则对向量或矩阵中的每个元素求其多项式的值。,例6.14 已知多项式x,4,+8x,3,-10，分别取x=1.2和一个23矩阵为自变量计算该多项式的值。,(1)代数多项式求值,25,(2)矩阵多项式求值,polyvalm函数用来求矩阵多项式的值，其调用格式与polyval相同，但含义不同。polyvalm函数要求x为方阵，它以方阵为自变量求多项式的值。设A为方阵，P代表多项式x,3,-5x,2,+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),例6.15 仍以多项式x,4,+8x,3,-10为例，取一个22矩阵为自变量分别用polyval和polyvalm计算该多项式的值。,(2)矩阵多项式求值,26,4. 多项式求根,n次多项式具有n个根，当然这些根可能是实根，也可能含有若干对共轭复根。MATLAB提供的roots函数用于求多项式的全部根，其调用格式为：,x=roots(P),其中P为多项式的系数向量，求得的根赋给向量x，即x(1),x(2),x(n)分别代表多项式的n个根。,4. 多项式求根,27,例6.16 求多项式x,4,+8x,3,-10的根。,命令如下：,A=1,8,0,0,-10;,x=roots(A),若已知多项式的全部根，则可以用poly函数建立起该多项式，其调用格式为：,P=poly(x),若x为具有n个元素的向量，则poly(x)建立以x为其根的多项式，且将该多项式的系数赋给向量P。,例6.16 求多项式x4+8x3-10的根。,28,例6.17 已知 f(x),(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),例6.17 已知 f(x),29,6.2 数值微积分,6.2.1 数值微分,1. 数值差分与差商,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，按行计算差分。,6.2 数值微积分6.2.1 数值微分,30,例6.18 设x由0,2间均匀分布的10个点组成，求sin,x,的13阶差分。,命令如下,：,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),例6.18 设x由0,2间均匀分布的10个点组成，求,31,例6.19 用不同的方法求函数f(x)的数值导数，并在同一个坐标系中做出f(x)的图像。,程序如下：,f=inline(sqrt(x.3+2*x.2-x+12)+(x+5).(1/6)+5*x+2);,g=inline(3*x.2+4*x-1)./sqrt(x.3+2*x.2-x+12)/2+1/6./(x+5).(5/6)+5);,x=-3:0.01:3;,p=polyfit(x,f(x),5); %用5次多项式p拟合f(x),dp=polyder(p); %对拟合多项式p求导数dp,dpx=polyval(dp,x); %求dp在假设点的函数值,dx=diff(f(x,3.01)/0.01; %直接对f(x)求数值导数,gx=g(x); %求函数f的导函数g在假设点的导数,plot(x,dpx,x,dx,.,x,gx,-); %作图,例6.19 用不同的方法求函数f(x)的数值导数，并在同一,32,6.2.2 数值积分,1. 数值积分基本原理,求解定积分的数值方法多种多样，如简单的梯形法、辛普生(Simpson)法、牛顿柯特斯(Newton-Cotes)法等都是经常采用的方法。它们的基本思想都是将整个积分区间a,b分成n个子区间x,i,x,i+1,，i=1,2,n，其中x,1,=a，x,n+1,=b。这样求定积分问题就分解为求和问题。,6.2.2 数值积分,33,2. 数值积分的实现,(1)被积函数是一个解析式,MATLAB提供了quad函数和quadl函数来求定积分。它们的调用格式为：,quad(filename,a,b,tol,trace),quadl(filename,a,b,tol,trace),2. 数值积分的实现,34,例6.20 用两种不同的方法求定积分。,先建立一个函数文件ex.m：,function ex=ex(x),ex=exp(-x.2);,然后在MATLAB命令窗口,，输入命令：,format long,I=quad(ex,0,1) %注意函数名应加字符引号,I =,0.74682418072642,I=quadl(ex,0,1),I =,0.74682413398845,也可不建立关于被积函数的函数文件，而使用语句函数(内联函数)求解，命令如下：,g=inline(exp(-x.2); %定义一个语句函数g(x)=exp(-x2),I=quadl(g,0,1) %注意函数名不加号,I =,0.74682413398845,format short,例6.20 用两种不同的方法求定积分。,35,(2)被积函数由一个表格定义,在科学实验和工程应用中，函数关系往往是不知道的，只有实验测定的一组样本点和样本值，这时，就无法使用quad函数计算其定积分。在MATLAB中，对由表格形式定义的函数关系的求定积分问题用trapz(X,Y)函数。其中向量X、Y定义函数关系Y=f(X)。X、Y是两个等长的向量：X=(x1,x2,xn)，Y=(y1,y2,yn)，并且x1x2 chol,Matrix must be positive definite,命令执行时，出现错误信息，说明A为非正定矩阵。,例6.27 用Cholesky分解求解例6.24中的线性方,54,6.4.2 迭代解法,迭代解法非常适合求解大型系数矩阵的方程组。在数值分析中，迭代解法主要包括 Jacobi迭代法、Gauss-Serdel迭代法、超松弛迭代法和两步迭代法。,1Jacobi迭代法,对于线性方程组Ax=b，如果A为非奇异方阵，即aii0(i=1,2,n)，则可将A分解为A=D-L-U，其中D为对角阵，其元素为A的对角元素，L与U为A的下三角阵和上三角阵，于是Ax=b化为：,x=D,-1,(L+U)x+D,-1,b,与之对应的迭代公式为：,x(k+1)=D,-1,(L+U)x(k)+D,-1,b,这就是Jacobi迭代公式。如果序列x(k+1)收敛于x，则x必是方程Ax=b的解。,6.4.2 迭代解法,55,Jacobi迭代法的MATLAB函数文件Jacobi.m如下：,function y,n=jacobi(A,b,x0,eps),if nargin=3,eps=1.0e-6;,elseif nargin=eps,x0=y;,y=B*x0+f;,n=n+1;,end,Jacobi迭代法的MATLAB函数文件Jacobi.m如下,56,例6.28 用Jacobi迭代法求解线性方程组。设迭代初值为0，迭代精度为10,-6,。,在命令中调用函数文件Jacobi.m，命令如下：,A=10,-1,0;-1,10,-2;0,-2,10;,b=9,7,6;,x,n=jacobi(A,b,0,0,0,1.0e-6),例6.28 用Jacobi迭代法求解线性方程组。设迭代初值,57,2Gauss-Serdel迭代法,在Jacobi迭代过程中，计算时，已经得到，不必再用，即原来的迭代公式Dx,(k+1),=(L+U)x,(k),+b可以改进为Dx,(k+1),=Lx,(k+1,)+Ux,(k),+b，于是得到：,x,(k+1),=(D-L),-1,Ux,(k),+(D-L),-1,b,该式即为Gauss-Serdel迭代公式。和Jacobi迭代相比，Gauss-Serdel迭代用新分量代替旧分量，精度会高些。,2Gauss-Serdel迭代法,58,Gauss-Serdel迭代法的MATLAB函数文件gauseidel.m如下：,function y,n=gauseidel(A,b,x0,eps),if nargin=3,eps=1.0e-6;,elseif nargin=eps,x0=y;,y=G*x0+f;,n=n+1;,end,Gauss-Serdel迭代法的MATLAB函数文件gaus,59,例6.29 用Gauss-Serdel迭代法求解下列线性方程组。设迭代初值为0，迭代精度为10,-6,。,在命令中调用函数文件gauseidel.m，命令如下：,A=10,-1,0;-1,10,-2;0,-2,10;,b=9,7,6;,x,n=gauseidel(A,b,0,0,0,1.0e-6),例6.29 用Gauss-Serdel迭代法求解下列线性方,60,例6.30 分别用Jacobi迭代和Gauss-Serdel迭代法求解下列线性方程组，看是否收敛。,命令如下：,a=1,2,-2;1,1,1;2,2,1;,b=9;7;6;,x,n=jacobi(a,b,0;0;0),x,n=gauseidel(a,b,0;0;0),例6.30 分别用Jacobi迭代和Gauss-Serde,61,6.4.3 求线性方程组的通解,线性方程组的求解分为两类：一类是求方程组的惟一解即特解，另一类是求方程组的无穷解即通解。这里对线性方程组 Ax=b的求解理论作一个归纳。,(1)当系数矩阵A是一个满秩方阵时，方程Ax=b称为恰定方程，方程有惟一解x=A-1b，这是最基本的一种情况。一般用x=Ab求解速度更快。,(2)当方程组右端向量b=0时，方程称为齐次方程组。齐次方程组总有零解，因此称解x=0为平凡解。当系数矩阵A的秩小于n(n为方程组中未知变量的个数)时，齐次方程组有无穷多个非平凡解，其通解中包含n-rank(A)个线性无关的解向量，用MATLAB的函数null(A,r)可求得基础解系。,6.4.3 求线性方程组的通解,62,(3)当方程组右端向量b0时，系数矩阵的秩rank(A)与其增广矩阵的秩rank(A,b)是判断其是否有解的基本条件：,当rank(A)=rank(A,b)=n时，方程组有惟一解：x=Ab 或 x=pinv(A)*b。,当rank(A)=rank(A,b)n时，方程组有无穷多个解，其通解=方程组的一个特解+对应的齐次方程组Ax=0的通解。可以用Ab求得方程组的一个特解，用null(A,r)求得该方程组所对应的齐次方程组的基础解系，基础解系中包含n-rank(A)个线性无关的解向量。,当rank(A)rank(A,b)时，方程组无解。,(3)当方程组右端向量b0时，系数矩阵的秩rank(A)与,63,有了上面这些讨论，可以设计一个求解线性方程组的函数文件line_solution.m。在例中可以调用line_solution.m文件来解线性方程组。,有了上面这些讨论，可以设计一个求解线性方程组的函数文件,64,6.5 非线性方程与最优化问题求解,6.5.1 非线性方程数值求解,1. 单变量非线性方程求解,在MATLAB中提供了一个fzero函数，可以用来求单变量非线性方程的根。该函数的调用格式为：,z=fzero(fname,x0,tol,trace),其中fname是待求根的函数文件名，x0为搜索的起点。一个函数可能有多个根，但fzero函数只给出离x0最近的那个根。tol控制结果的相对精度，缺省时取tol=eps，trace指定迭代信息是否在运算中显示，为1时显示，为0时不显示，缺省时取trace=0。,6.5 非线性方程与最优化问题求解,65,例6.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.1926,fzero(fz,1) %,以,1,作为迭代初值,ans =,0.1926,例6.33 求f(x)在x0=-5和x0=1作为迭代初值时,66,2. 非线性方程组的求解,对于非线性方程组F(X)=0，用fsolve函数求其数值解。fsolve函数的调用格式为：,X=fsolve(fun,X0,option),其中X为返回的解，fun是用于定义需求解的非线性方程组的函数文件名，X0是求根过程的初值，option为最优化工具箱的选项设定。最优化工具箱提供了20多个选项，用户可以使用optimset命令将它们显示出来。如果想改变其中某个选项，则可以调用optimset()函数来完成。例如，Display选项决定函数调用时中间结果的显示方式，其中off为不显示，iter表示每步都显示，final只显示最终结果。optimset(Display,off)将设定Display选项为off。,2. 非线性方程组的求解,67,例6.34 求下列方程组在(1，1，1)附近的解并对结果进行验证。,首先建立函数文件myfun.m。,function F=myfun (X),x=X(1);,y=X(2);,z=X(3);,F(1)=sin(x)+y+z2*exp(x);,F(2)=x+y+z;,F(3)=x*y*z;,在给定的初值x0=1，y0=1，z0=1下，调用fsolve函数求方程的根。,X=fsolve(myfun,1,1,1,optimset(Display, off),X =,0.0224 -0.0224 -0.0000,例6.34 求下列方程组在(1，1，1)附近的解并对结果进,68,6.5.2 无约束最优化问题求解,在实际应用中，许多科学研究和工程计算问题都可以归结为一个最小化问题，如能量最小、时间最短等。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)上的最大值。,6.5.2 无约束最优化问题求解,69,例6.36 求函数在区间(-10，1)和(1，10)上的最小值点。,首先建立函数文件fx.m：,function f=f(x),f=x-1/x+5;,上述函数文件也可用一个语句函数代替：,f=inline(x-1/x+5),再在MATLAB命令窗口，输入命令：,fminbnd(fx,-10,-1) %求函数在(-10，-1)内的最小值点和最小值,fminbnd(f,1,10) %求函数在(1，10)内的最小值点。注意函数名f不用加,例6.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命令窗口，输入命令：,U,fmin=fminsearch(fxyz,0.5,0.5,0.5) %求函数的最小值点和最小值,例6.36 求函数在区间(-10，1)和(1，10)上的最,70,6.5.3 有约束最优化问题求解,MATLAB最优化工具箱提供了一个fmincon函数，专门用于求解各种约束下的最优化问题。该函数的调用格式为：,x,fval=fmincon(filename,x0,A,b, Aeq,beq,Lbnd,Ubnd, NonF,option),其中x、fval、filename、x0和option的含义与求最小值函数相同。其余参数为约束条件，参数NonF为非线性约束函数的M文件名。如果某个约束不存在，则用空矩阵来表示。,例6.38 求解有约束最优化问题。,首先编写目标函数M文件fop.m。,function f=fop(x),f=0.4*x(2)+x(1)2+x(2)2-x(1)*x(2)+1/30*x(1)3;,再设定约束条件，并调用fmincon函数求解此约束最优化问题。,x0=0.5;0.5;,A=-1,-0.5;-0.5,-1;,b=-0.4;-0.5;,lb=0;0;,option=optimset; option.LargeScale=off; option.Display =off;,x,f=fmincon(fop ,x0,A,b,lb,option),6.5.3 有约束最优化问题求解,71,6.6 常微分方程的数值求解,6.6.1 龙格库塔法简介,6.6.2 龙格库塔法的实现,基于龙格库塔法，MATLAB提供了求常微分方程数值解的函数，一般调用格式为：,t,y=ode23(fname,tspan,y0),t,y=ode45(fname,tspan,y0),其中fname是定义f(t,y)的函数文件名，该函数文件必须返回一个列向量。tspan形式为t0,tf,表示求解区间。y0是初始状态列向量。t和y分别给出时间向量和相应的状态向量。,6.6 常微分方程的数值求解,72,例6.39 设有初值问题，试求其数值解，并与精确解相比较(精确解为y(t)=)。,(1) 建立函数文件funt.m。,function yp=funt(t,y),yp=(y2-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为精确值，显然两者近似。,例6.39 设有初值问题，试求其数值解，并与精确解相比较(,73,例,6.40 已知一个二阶线性系统的微分方程，绘制系统的时间响应曲线和相平面图。,函数ode23和ode45是对一阶常微分方程组设计的，因此对高阶常微分方程，需先将它转化为一阶常微分方程组，即状态方程。,建立一个函数文件sys.m：,function xdot=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,axis equal,例6.40 已知一个二阶线性系统的微分方程，绘制系统的时间,74,6.7 稀疏矩阵6.7.1 矩阵存储方式MATLAB的矩阵有两种存储方式：完全存储方式和稀疏存储方式。1完全存储方式完全存储方式是将矩阵的全部元素按列存储。以前讲到的矩阵的存储方式都是按这个方式存储的，此存储方式对稀疏矩阵也适用。,6.7 稀疏矩阵6.7.1 矩阵存储方式MATLAB,75,2稀疏存储方式稀疏存储方式仅存储矩阵所有的非零元素的值及其位置，即行号和列号。在MATLAB中，稀疏存储方式也是按列存储的。注意，在讲稀疏矩阵时，有两个不同的概念，一是指矩阵的0元素较多，该矩阵是一个具有稀疏特征的矩阵，二是指采用稀疏方式存储的矩阵。,2稀疏存储方式稀疏存储方式仅存储矩阵所有的非零元素的值及,76,6.7.2 稀疏存储方式的产生1将完全存储方式转化为稀疏存储方式函数A=sparse(S)将矩阵S转化为稀疏存储方式的矩阵A。当矩阵S是稀疏存储方式时，则函数调用相当于A=S。sparse函数还有其他一些调用格式：sparse(m,n)：生成一个mn的所有元素都是0的稀疏矩阵。sparse(u,v,S)：u,v,S是3个等长的向量。S是要建立的稀疏矩阵的非0元素，u(i)、v(i)分别是S(i)的行和列下标，该函数建立一个max(u)行、max(v)列并以S为稀疏元素的稀疏矩阵。此外，还有一些和稀疏矩阵操作有关的函数。例如u,v,S=find(A)：返回矩阵A中非0元素的下标和元素。这里产生的u,v,S可作为sparse(u,v,S)的参数。full(A)：返回和稀疏存储矩阵A对应的完全存储方式矩阵。,6.7.2 稀疏存储方式的产生1将完全存储方式转化为稀,77,2产生稀疏存储矩阵只把要建立的稀疏矩阵的非0元素及其所在行和列的位置表示出来后由MATLAB自己产生其稀疏存储，这需要使用spconvert函数。调用格式为：B=spconvert(A)其中A为一个m3或m4的矩阵，其每行表示一个非0元素，m是非0元素的个数，A每个元素的意义是：(i,1) 第i个非0元素所在的行。(i,2) 第i个非0元素所在的列。(i,3) 第i个非0元素值的实部。(i,4) 第i个非0元素值的虚部，若矩阵的全部元素都是实数，则无须第四列。该函数将A所描述的一个稀疏矩阵转化为一个稀疏存储矩阵。,2产生稀疏存储矩阵只把要建立的稀疏矩阵的非0元素及其所在,78,例6.42 根据表示稀疏矩阵的矩阵A，产生一个稀疏存储方式矩阵B。命令如下：A=2,2,1;3,1,-1;4,3,3;5,3,8;6,6,12;B=spconvert(A),例6.42 根据表示稀疏矩阵的矩阵A，产生一个稀疏存储方式,79,3带状稀疏存储矩阵用spdiags函数产生带状稀疏矩阵的稀疏存储，调用格式是：A=spdiags(B,d,m,n)其中，参数m,n为原带状矩阵的行数与列数。B为rp阶矩阵，这里r=min(m,n)，p为原带状矩阵所有非零对角线的条数，矩阵B的第i列即为原带状矩阵的第i条非零对角线。,3带状稀疏存储矩阵用spdiags函数产生带状稀疏矩阵的,80,4单位矩阵的稀疏存储单位矩阵只有对角线元素为1，其他元素都为0，是一种具有稀疏特征的矩阵。函数eye产生一个完全存储方式的单位矩阵。MATLAB还有一个产生稀疏存储方式的单位矩阵的函数，这就是speye。函数speye(m,n)返回一个mn的稀疏存储单位矩阵。,4单位矩阵的稀疏存储单位矩阵只有对角线元素为1，其他元素,81,6.7.3 稀疏矩阵应用举例稀疏存储矩阵只是矩阵的存储方式不同，它的运算规则与普通矩阵是一样的。所以，在运算过程中，稀疏存储矩阵可以直接参与运算。当参与运算的对象不全是稀疏存储矩阵时，所得结果一般是完全存储形式。,6.7.3 稀疏矩阵应用举例稀疏存储矩阵只是矩阵的存储方,82,第6章MATLAB数值计算61-数据处理与多项式课件,83,