第6章 MATLAB数值计算

1、第第6 6章章 MATLABMATLAB数值计算数值计算XIAN UNIVERSITYOF TECHNOLOGYPage 26.1.1 数据统计与分析数据统计与分析1. 求矩阵最大元素和最小元素求矩阵最大元素和最小元素（1）求向量的最大值和最小值）求向量的最大值和最小值 y=max(X)：返回向量：返回向量X的最大值存入的最大值存入y，如果，如果X中包含复数元中包含复数元素，则按模取最大值。素，则按模取最大值。 y,I=max(X)：返回向量：返回向量X的最大值存入的最大值存入y，最大值的序号存入，最大值的序号存入I，如果如果X中包含复数元素，则按模取最大值。中包含复数元素，则按模取最大值。求

2、向量求向量X的最小值的函数是的最小值的函数是min(X)，用法和，用法和max(X)完全相同。完全相同。6.1 数据处理与多项式计算数据处理与多项式计算 XIAN UNIVERSITYOF TECHNOLOGYPage 3例例 求向量求向量x的最大值。的最大值。命令如下：命令如下：x=-43,72,9,16,23,47;y=max(x) %求向量求向量x中的最大值中的最大值y,l=max(x) %求向量求向量x中的最大值及其该元素的位置中的最大值及其该元素的位置多输入函数多输入函数 例如：例如： x=max(z , y) 多输出函数多个输出值用多输出函数多个输出值用 括起来，且输出值之间用逗括

3、起来，且输出值之间用逗号隔开。号隔开。6.1 数据处理与多项式计算数据处理与多项式计算 XIAN UNIVERSITYOF TECHNOLOGYPage 4（2）求矩阵的最大值和最小值）求矩阵的最大值和最小值 max(A)：返回一个行向量，向量的第：返回一个行向量，向量的第i个元素是矩阵个元素是矩阵A的第的第i列列上的最大值。上的最大值。 Y,U=max(A)：返回行向量：返回行向量Y和和U，Y向量记录向量记录A的每列的最的每列的最大值，大值，U向量记录每列最大值的行号。向量记录每列最大值的行号。 max(A,dim)：dim取取1或或2。dim取取1时，该函数和时，该函数和max(A)完完全

4、相同；全相同；dim取取2时，该函数返回一个列向量，其第时，该函数返回一个列向量，其第i个元素是个元素是A矩阵的第矩阵的第i行上的最大值。行上的最大值。6.1 数据处理与多项式计算数据处理与多项式计算 XIAN UNIVERSITYOF TECHNOLOGYPage 56.1 数据处理与多项式计算数据处理与多项式计算 A=13,-56,78;25,63,-235;78,25,563;1,0,-1; a1 = max(A,2); a2 = min(A,2); a3 = max(A); a4 = min(A); a5 = max(max(A); a6 = min(A(:);XIAN UNIVERS

5、ITYOF TECHNOLOGYPage 6（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

6、中的较大者。中的较大者。 例例 求两个求两个23矩阵矩阵x, y所有同一位置上的较大元素构成的新所有同一位置上的较大元素构成的新矩阵矩阵p。6.1 数据处理与多项式计算数据处理与多项式计算 x=4,5,6;1,4,8x = 4 5 6 1 4 8 y=1,7,5;4,5,7y = 1 7 5 4 5 7 p=max(x,y)p = 4 7 6 4 5 8 p2=max(y,5)p2 = 5 7 5 5 5 7XIAN UNIVERSITYOF TECHNOLOGYPage 72. 求矩阵的平均值和中值求矩阵的平均值和中值mean(X) 返回向量返回向量X的算术平均值的算术平均值median(X

7、) 返回向量返回向量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的

8、第的第i行的中值行的中值6.1 数据处理与多项式计算数据处理与多项式计算 XIAN UNIVERSITYOF TECHNOLOGYPage 86.1 数据处理与多项式计算数据处理与多项式计算 3. 矩阵元素求和与求积矩阵元素求和与求积sum(X) prod(X)7. 排序排序sort(X)XIAN UNIVERSITYOF TECHNOLOGYPage 96.1 数据处理与多项式计算数据处理与多项式计算 A = 2 5 4 ;2 3 6; 9 6 1A = 2 5 4 2 3 6 9 6 1 mean(A)ans = 4.3333 4.6667 3.6667 median(A)ans = 2

9、5 4 sort(A)ans = 2 3 1 2 5 4 9 6 6 sort(A,2)ans = 2 4 5 2 3 6 1 6 9 sum(A)ans = 13 14 11 prod(A)ans = 36 90 24XIAN UNIVERSITYOF TECHNOLOGYPage 106.1.2 数据插值数据插值 一维数据插值一维数据插值插值：是在认定所给插值：是在认定所给“基准数据基准数据”完全完全正确的情况下，研正确的情况下，研 究如何究如何“平滑平滑”的的估算出估算出“基准数据基准数据”之间其它点的函数之间其它点的函数值。因此，插值所值。因此，插值所得曲线一定得曲线一定穿过穿过“基准

10、数据基准数据”。6.1 数据处理与多项式计算数据处理与多项式计算 XIAN UNIVERSITYOF TECHNOLOGYPage 11 函数根据函数根据x,y的值，计算函数在的值，计算函数在xs处的值。处的值。x,y是两个等长是两个等长的已知向量，分别描述采样点和样本值，的已知向量，分别描述采样点和样本值，xs是一个向量或是一个向量或标量，描述欲插值的点，标量，描述欲插值的点，ys是一个与是一个与xs等长的插值结果。等长的插值结果。 其中其中:(1) x,y是量测数据对；是量测数据对； (2) xs是需要内插的点所构成的向量。是需要内插的点所构成的向量。 (3) method是指所使用的内插

11、方法。是指所使用的内插方法。ys=interp1(x,y,xs,method);说明：说明：interp1仅是插仅是插值指令的一种，还有值指令的一种，还有interp2 、interp3等。等。 插值算法：插值算法：nearest，linear，spline，cubic注意：注意：xs的取值范围不能超出的取值范围不能超出x的给定范围，否则，会给出的给定范围，否则，会给出“NaN”错误。错误。6.1 数据处理与多项式计算数据处理与多项式计算 XIAN UNIVERSITYOF TECHNOLOGYPage 12%给定数据对x0=0:0.1:1;y0=-.44,1.97,3.11,5.25,5.0

12、2,4.66,4.01,4.58,3.45,5.35,9.22; %采用三次多项式进行插值xi=0:0.02:1;yi=interp1(x0,y0,xi,cubic); %绘图plot(xi,yi,-b,x0,y0,.r,MarkerSize,20),xlabel(x) 6.1 数据处理与多项式计算数据处理与多项式计算 XIAN UNIVERSITYOF TECHNOLOGYPage 13例例:某观测站测得某日某观测站测得某日6:00时至时至18:00时之间每隔时之间每隔2小时的室小时的室内外温度内外温度()，用，用3次样条插值分别求得该日室内外次样条插值分别求得该日室内外6:00至至18:0

13、0时之间每隔时之间每隔0.5小时各点的近似温度小时各点的近似温度()。 设时间变量设时间变量h为一行向量，温度变量为一行向量，温度变量t为一个两列矩阵，为一个两列矩阵，其中第一列存放室内温度，第二列储存室外温度。命令如其中第一列存放室内温度，第二列储存室外温度。命令如下：下： h =6:2:18; t=18,20,22,25,30,28,24;15,19,24,28,34,32,30; XI =6:0.5:18; YI=interp1(h,t,XI, spline); %用用3次样条插值计算次样条插值计算 plot(XI,YI, -,h,t, .,MarkerSize,20),xlabel(时

14、间时间)6.1 数据处理与多项式计算数据处理与多项式计算 XIAN UNIVERSITYOF TECHNOLOGYPage 141）最邻近插值方法（nearest）插值点的值与其最邻近的点的函数值相等。*6.1 数据处理与多项式计算数据处理与多项式计算 XIAN UNIVERSITYOF TECHNOLOGYPage 15 2）线性插值方法（ linear ）插值点的值在前，后两个数据点所构成的直线上。*6.1 数据处理与多项式计算数据处理与多项式计算 XIAN UNIVERSITYOF TECHNOLOGYPage 16 3）三次样条插值方法(spline)利用一系列样条函数获得内插数据点，

15、从而确定已有数据点之间的函数。4）三次曲线内插方法（cubic）构造三次曲线函数来拟合已知数据x、y，从而确定内插点的值。6.1 数据处理与多项式计算数据处理与多项式计算 XIAN UNIVERSITYOF TECHNOLOGYPage 17说明:1、四种插值方法中，x中的数据是单调但不一定距均匀的。 2、若已知x为均匀的，则在method前加*，可使执行速度加快。3、按nearest linear cubic spline的顺序，对内存要求 从小到大，执行速度由快到慢，平滑度由差到好。 6.1 数据处理与多项式计算数据处理与多项式计算 XIAN UNIVERSITYOF TECHNOLOGY

16、Page 186.1.3 曲线拟合曲线拟合 p=polyfit(x,y,n); 1121.nnnnya xa xa xa多项式拟合多项式拟合6.1 数据处理与多项式计算数据处理与多项式计算 XIAN UNIVERSITYOF TECHNOLOGYPage 196.1 数据处理与多项式计算数据处理与多项式计算 XIAN UNIVERSITYOF TECHNOLOGYPage 20 在在MATLAB中，用中，用polyfit函数来求得最小二乘函数来求得最小二乘拟合多项式的系数，再用拟合多项式的系数，再用polyval函数按所得的多项函数按所得的多项式计算所给出的点上的函数近似值。式计算所给出的点上

17、的函数近似值。 p=polyfit(x,y,n); 函数根据采样点函数根据采样点x和采样点函数值和采样点函数值y，产生一个，产生一个n次多次多项式项式p。其中。其中x,y是两个等长的向量，是两个等长的向量，p是一个长度为是一个长度为n+1的的向量，向量，p的元素为多项式系数。的元素为多项式系数。 polyval函数的功能是按多项式的系数计算函数的功能是按多项式的系数计算x点多项式点多项式的值的值.6.1 数据处理与多项式计算数据处理与多项式计算 XIAN UNIVERSITYOF TECHNOLOGYPage 216.1 数据处理与多项式计算数据处理与多项式计算 %给定数据对x0=0:0.1:

18、1;y0=-.44,1.97,3.11,5.25,5.02,4.66,4.01,4.58,3.45,5.35,9.22;%求拟合多项式系数n=3;P=polyfit(x0,y0,n)%图示拟合情况 xx=0:0.01:1;yy=polyval(P,xx);plot(xx,yy,-b,x0,y0,.r,MarkerSize,20),xlabel(x) 拟合：寻找一条“平滑”曲线来最好的表现带噪声的“测量数据”。但并不要求拟合曲线穿过这些“测量数据”点。 XIAN UNIVERSITYOF TECHNOLOGYPage 22例例6.11 用一个用一个3次多项式在区间次多项式在区间0，2内逼近函数。

19、内逼近函数。命令如下命令如下：X=linspace(0,2*pi,50);Y=sin(X);P=polyfit(X,Y,3) %得到得到3次多项式的系数和误差次多项式的系数和误差Y1 = polyval(P,X);plot(X,Y,:o,X,Y1,-*);legend(原始曲线原始曲线,拟合曲线拟合曲线);6.1 数据处理与多项式计算数据处理与多项式计算 XIAN UNIVERSITYOF TECHNOLOGYPage 236.1.4 多项式计算多项式计算1、多项式的表示：、多项式的表示： MATALB中，用一个向量来表示多项式。这个向量中按照降幂的顺序排列多项式的各项系数。52)(3xxxp

20、在MATLAB中为：p=1,0,-2,5;说明：如果多项式中缺某幂次项，则认为该项的系数为零。6.1 数据处理与多项式计算数据处理与多项式计算 XIAN UNIVERSITYOF TECHNOLOGYPage 241. 多项式的四则运算多项式的四则运算（1）多项式的加减运算）多项式的加减运算（2）多项式乘法运算）多项式乘法运算因为多项式是用其系数构成的离散序列表示的，因此多项式因为多项式是用其系数构成的离散序列表示的，因此多项式的乘法和除法就对应于卷积和去卷积的操作。的乘法和除法就对应于卷积和去卷积的操作。6.1 数据处理与多项式计算数据处理与多项式计算 32( )(253)(61)p xxx

21、xx a=1,-2,5,3; b=0,0,6,-1; c=a+bc = 1 -2 11 2XIAN UNIVERSITYOF TECHNOLOGYPage 251）多项式乘法）多项式乘法卷积：卷积： c=conv(a,b);例如: (x-1)(x-2)=x2-3x+2求多项式求多项式a和和b的乘积的乘积/卷积。卷积。a、b是两个多项式系数向是两个多项式系数向量。量。 a=1,-1; b=1,-2; c=conv(a,b)c = 1 -3 26.1 数据处理与多项式计算数据处理与多项式计算 XIAN UNIVERSITYOF TECHNOLOGYPage 26)()()()()(sbsrsqsb

22、sc例如：例如：q,r=deconv(c,b); (1)(2)33(1)22xxxxx c=1,-3,5; b=1,-2; q,r=deconv(c,b)q = 1 -1r = 0 0 3 conv(q,b)+rans = 1 -3 5其中其中q返回多项式返回多项式c除以除以b的商式，的商式，r返回余式。返回余式。q和和r仍是多项式系仍是多项式系数向量。数向量。deconv是是conv的逆函数，的逆函数，即有即有c=conv(b,q)+r。（3）多项式除法）多项式除法解卷积解卷积6.1 数据处理与多项式计算数据处理与多项式计算 XIAN UNIVERSITYOF TECHNOLOGYPage

23、272. 多项式的导函数多项式的导函数)()(spdsdsq q=polyder(p);具体用法:(1)q=polyder(p)可以得到多项式p的导数。(2) q=polyder(a,b)可以求出多项式a,b之积的导数。(3)q,d=polyder(b,a)可以求出多项式之比b(s)/a(s)的导数。)()()()(sasbdsdsdsq)()()(sbsadsdsq6.1 数据处理与多项式计算数据处理与多项式计算 c=1,-3,5; polyder(c)ans = 2 -3XIAN UNIVERSITYOF TECHNOLOGYPage 283. 多项式求值多项式求值 y=polyval(p

24、,x); y=polyvalm(p,x);功能：按数组运算规则计算多项式的值。其中x可以是标量和数组。功能：按矩阵运算规则计算多项式的值。其中x必须为方阵。6.1 数据处理与多项式计算数据处理与多项式计算 XIAN UNIVERSITYOF TECHNOLOGYPage 29P= 2 3 1 ;X= 1 2 3 ; 3 2 1 ; 2 1 3 ;y=polyval ( P , X ) y=polyvalm ( P , X ) y=2*X.*X+3*X+1y=2*X*X+3*X+eye(3)6.1 数据处理与多项式计算数据处理与多项式计算 XIAN UNIVERSITYOF TECHNOLOGY

25、Page 304. 多项式求根多项式求根R=roots(P)0.1121nnnnpxpxpxp p=1,-3,2; R = roots(p)R = 2 1p1=poly(R)功能：计算多 项式P的根。6.1 数据处理与多项式计算数据处理与多项式计算 若已知多项式的全部根，则可以用若已知多项式的全部根，则可以用poly函数建函数建立起该多项式：立起该多项式：P=poly(x)当当x为具有为具有n个元素的向量，则个元素的向量，则poly(x)建立以建立以x为其根的多项式，且将该多项式的系数赋给为其根的多项式，且将该多项式的系数赋给向量向量P。XIAN UNIVERSITYOF TECHNOLOGY

26、Page 316.2 数值微积分数值微积分6.2.1 数值微分数值微分1. 数值差分与差商数值差分与差商2. 数值微分的实现数值微分的实现在在MATLAB中，没有直接提供求数值导数的函数，只有计算向中，没有直接提供求数值导数的函数，只有计算向前差分的函数前差分的函数diff，其调用格式为：，其调用格式为：DX=diff(X)：计算向量：计算向量X的向前差分的向前差分DX=diff(X,n)：计算：计算X的的n阶向前差分。阶向前差分。XIAN UNIVERSITYOF TECHNOLOGYPage 32clear;clc;p = 1 2 1;x=-1:0.01:1;dp = polyder(p)

27、;dpx = polyval(dp,x);y = x.2 + 2*x +1;dx = diff(y)/0.01;x1 = x(1:end-1);gx = 2*x+2;plot(x,dpx,r,x1,dx,g,x,gx,b);figure(2)subplot(1,3,1); plot(x,dpx,r);subplot(1,3,2); plot(x1,dx,g);subplot(1,3,3); plot(x,gx,b);例：例：f(x) = x2+2x+1用不同方法求其倒数，用不同方法求其倒数，并画出图象并画出图象XIAN UNIVERSITYOF TECHNOLOGYPage 33(1)被积函数

28、是一个解析式被积函数是一个解析式例：求以下的定积分。其精确值为0.7468204102dxeIx%（1）使用字符串表示被处理函数 fun=exp(-x.*x);%（2）调用积分指令求积分值 E1 = quad (fun,0,1) E2 = quadl (fun,0,1)区别：求数值积分是采取不同的方法6.2 数值微积分数值微积分6.2.2 数值积分数值积分2. 数值积分的实现数值积分的实现XIAN UNIVERSITYOF TECHNOLOGYPage 346.4.1 直接解法直接解法1利用左除运算符的直接解法利用左除运算符的直接解法对于线性方程组对于线性方程组Ax=b，可以利用左除运算符，可

29、以利用左除运算符“”求解：求解： x=Ab6.4 线性方程组求解线性方程组求解XIAN UNIVERSITYOF TECHNOLOGYPage 356.4.3 求线性方程组的通解求线性方程组的通解 线性方程组的求解分为两类：一类是求方程组的惟一线性方程组的求解分为两类：一类是求方程组的惟一解即特解，另一类是求方程组的无穷解即通解。这里对线解即特解，另一类是求方程组的无穷解即通解。这里对线性方程组性方程组 Ax=b的求解理论作一个归纳：的求解理论作一个归纳：p1726.4 线性方程组求解线性方程组求解XIAN UNIVERSITYOF TECHNOLOGYPage 366.5.1 非线性方程数值

30、求解非线性方程数值求解1. 单变量非线性方程求解单变量非线性方程求解z=fzero(fname,x0,tol,trace) 求单变量非线性方程的根求单变量非线性方程的根其中其中fname是待求根的函数文件名，是待求根的函数文件名，x0为搜索的起点。一个为搜索的起点。一个函数可能有多个根，但函数可能有多个根，但fzero函数只给出离函数只给出离x0最近的那个根。最近的那个根。tol控制结果的相对精度，缺省时取控制结果的相对精度，缺省时取tol=eps，trace 指定指定迭代信息是否在运算中显示，为迭代信息是否在运算中显示，为1时显示，为时显示，为0时不显示，时不显示，缺省时取缺省时取trace

31、=0。6.5 非线性方程与最优化问题求解非线性方程与最优化问题求解XIAN UNIVERSITYOF TECHNOLOGYPage 37例例6.33 求求f(x)=x-1/x+5在在x0=-5和和x0=1作为迭代初值时的零点。作为迭代初值时的零点。先建立函数文件先建立函数文件fz.m：function f=fz(x)f=x-1/x+5;然后调用然后调用fzero函数求根。：函数求根。：fzero(fz,-5) %以以-5作为迭代初值作为迭代初值fzero(fz,1) %以以1作为迭代初值作为迭代初值x=-6:0.1:2;f=x-1./x+5; plot(x,f,b-,x,f,r.)hold onplot(x,zeros(1,length(x),k-)6.5 非线性方程与最优化问题求解非线性方程与最优化问题求解XIAN UNIVERSITYOF TECHNOLOGYPage