1、MATLAB SDU,1,第五章 MATLAB的数值计算,1.多项式计算 2. 方程组求解 3. 函数及其数值分析 4.数据处理 5.数据插值 6.曲线拟合 7.傅里叶变换,MATLAB SDU,2,1.多项式 运算,多项式表达方式的约定 创建多项式的方法 多项式运算函数,MATLAB SDU,3,多项式表达方式的约定,MATLAB约定降幂多项式P(x)=a0xn+ a1xn-1 +an-1x+an用以下系数矢量（系数行向量）表示：p= a0, a1 ,an-1,an,即把多项式的各项系数依降幂次序排放在行向量的元素位置上。 注意I：假如多项式中缺某幂次项，则应认为该幂次项的系数为零。注意II
2、：多项式加减时低阶多项式必须首零填补使得与高阶多项式阶次相同。,MATLAB SDU,4,创建多项式的方法,系数矢量的直接输入法 在命令窗直接输入多项式的系数矢量，由根矢量创建多项式 由给定的根矢量创建多项式，由函数poly实现。特征多项式输入法 由矩阵的特征多项式取得，由函数poly实现。多项式的字符标示：poly2str,MATLAB SDU,5,多项式求根,n次多项式具有n个根，当然这些根可能是实根，也可能含有若干对共轭复根。MATLAB提供的roots函数用于求多项式的全部根，其调用格式为： x=roots(P) 其中P为多项式的系数向量，求得的根赋给向量x，即x(1),x(2),x(
3、n)分别代表多项式的n个根。,MATLAB SDU,6,例: 求多项式x4+8x3-10的根。 命令如下： A=1,8,0,0,-10; x=roots(A) 若已知多项式的全部根，则可以用poly函数建立起该多项式，其调用格式为： P=poly(x) 若x为具有n个元素的向量，则poly(x)建立以x为其根的多项式，且将该多项式的系数赋给向量P。,MATLAB SDU,7,多项式运算：乘运算,例:a(x)=x2+2x+3; b(x)=4x2+5x+6; c = (x2+2x+3)(4x2+5x+6) a=1 2 3;b=4 5 6; c=conv(a,b)=conv(1 2 3,4 5 6)
4、 c = 4.00 13.00 28.00 27.00 18.00 p=poly2str(c,x) p = 4 x4 + 13 x3 + 28 x2 + 27 x + 18,多项式加减法：补零使阶数相等 多项式乘法运算 conv 函数conv(P1,P2)用于求多项式P1和P2的乘积。这里，P1、P2是两个多项式系数向量。,多项式乘除对应向量卷积和解卷积,MATLAB SDU,8,多项式除运算deconv,函数Q,r=deconv(P1,P2)用于对多项式P1和P2作除法运算。其中Q返回多项式P1除以P2的商式，r返回P1除以P2的余式。这里，Q和r仍是多项式系数向量。 deconv是conv
5、的逆函数，即有P1=conv(P2,Q)+r。例：求多项式x4+8x3-10除以多项式2x2-x+3的结果。,MATLAB SDU,9,多项式微分,matlab提供了polyder函数多项式的微分。 命令格式： polyder(p): 求p的微分 polyder(a,b): 求多项式a,b乘积的微分 p,q=polyder(a,b): 求多项式a,b商的微分 例：a=1 2 3 4 5; poly2str(a,x) ans = x4 + 2 x3 + 3 x2 + 4 x + 5 b=polyder(a) b = 4 6 6 4 poly2str(b,x) ans =4 x3 + 6 x2 +
6、 6 x + 4,MATLAB SDU,10,多项式积分,多项式积分命令函数polyint 1 polyint(p) 多项式p的积分，常数项为0 2 polyint(p,k)多项式p的积分，常数项为k,MATLAB SDU,11,有理多项式的运算函数,多项式之比的分式展开 r,p,k=residue(a,b) 从展开分式得到多项式a,b=residue(r,p,k)注意residue函数的可逆性,MATLAB SDU,12,多项式的求值,MATLAB提供了两种求多项式值的函数：polyval与polyvalm，它们的输入参数均为多项式系数向量P和自变量x。两者的区别在于前者是代数多项式求值，而
7、后者是矩阵多项式求值。,MATLAB SDU,13,代数多项式求值,polyval函数用来求代数多项式的值，其调用格式为： Y=polyval(P,x) 若x为一数值，则求多项式在该点的值；若x为向量或矩阵，则对向量或矩阵中的每个元素求其多项式的值。例 已知多项式x4+8x3-10，分别取x=1.2和一个23矩阵为自变量计算该多项式的值。 例: 画出上例x在-5到5间的曲线,MATLAB SDU,14,矩阵多项式求值【了解】,polyvalm函数用来求矩阵多项式的值，其调用格式与polyval相同，但含义不同。polyvalm函数要求x为方阵，它以方阵为自变量求多项式的值。设A为方阵，P代表多
8、项式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)例 仍以多项式x4+8x3-10为例，取一个22矩阵为自变量分别用polyval和polyvalm计算该多项式的值。,MATLAB SDU,15,2. 方程组求解,MATLAB SDU,16,代数方程组求解,方法1 矩阵的除运算 （教材）左除求逆（了解：各种方程类型 恰定，超定，欠定）方法2 符号方程函数求解（符号计算章节）,MATLAB SDU,17,3 函数及其数值分析,函数定义方
9、法(1) m文件生成函数文件(2) 函数句柄教材：匿名函数handle = functionnamehandle = (arglist)anonymous_function(3) inline函数来定义 g = inline(expr) g = inline(expr,arg1,arg2,.) g = inline(expr,n),MATLAB SDU,18,函数曲线绘制,绘制函数曲线的专用函数fplot的调用 FPLOT(FUN,LIMS)特点：绘图数据由函数在指定范围内自适应产生，根据函数曲线的平滑程度自动调整数据点的密度,绘制函数曲线的一般方法，计算出函数在某一区间值，然后根据两组数据值
10、绘制出函数曲线，但是如果函数在某些区间是平坦无激励的，某些区间却是失控的，传统方法无法表达函数的真正特性,MATLAB SDU,19,函数的图形绘制,绘制函数fplotfplot(function,limits,LineSpec)FUN can be specified using , an inline object, or an expression: (注意引号)举例：subplot(2,2,1), fplot(humps,0 1)f = (x)abs(exp(-j*x*(0:9)*ones(10,1);subplot(2,2,2), fplot(f,0 2*pi)subplot(2,2
11、,3), fplot(tan(x),sin(x),cos(x),2*pi*-1 1 -1 1)subplot(2,2,4), fplot(sin(1 ./ x), 0.01 0.1),MATLAB SDU,20,函数极值,MATLAB中只存在处理极小值命令的函数fminbnd，极大值的处理等价于-f(x)的极小值 局域极值的函数调用：x = fminbnd(fun,x1,x2)：一元函数的x1,x2范围内极小值时x的取值x,fval = fminbnd(fun,x1,x2): 同时返回fval的值例题: 0,2内x3-2x-5的极值x = fminsearch(fun,x0): 单纯形法求函数
12、在x0点附近极值，x0为标量或者向量(多元函数)x,fval = fminsearch(fun,x0): 同时返回fval的值 X=fminunc(fun,X0,options): 拟牛顿法多元函数极值点,MATLAB SDU,21,函数零点,Matlab中用fzero来寻找单变量函数值为零的自变量的值，调用格式：x = fzero(fun,x0)x,fval=fzero(fun,x0)x0指定搜索的点注意： fzero并不一定能找到零点搜索方法：先猜测一个初时零点所在的区间；然后通过一些计算，使得猜测值不断精确，或者使得猜测区间不断收缩，直至达到预先指定的精度，终止计算。help fzero
13、,MATLAB SDU,22,数据统计处理,最大值和最小值 MATLAB提供的求数据序列的最大值和最小值的函数分别为max和min，两个函数的调用格式和操作过程类似。 1求向量的最大值和最小值 求一个向量X的最大值的函数有两种调用格式，分别是： (1) y=max(X)：返回向量X的最大值存入y，如果X中包含复数元素，则按模取最大值。,MATLAB SDU,23,(2) y,I=max(X)：返回向量X的最大值存入y，最大值的序号存入I，如果X中包含复数元素，则按模取最大值。 求向量X的最小值的函数是min(X)，用法和max(X)完全相同。例 求向量x的最大值。 x=-43,72,9,16,
14、23,47; y=max(x) %求向量x中的最大值 y,l=max(x) %求向量x中的最大值及其该元素的位置,MATLAB SDU,24,求矩阵的最大值和最小值,求矩阵A的最大值的函数有3种调用格式，分别是： (1) max(A)：返回一个行向量，向量的第i个元素是矩阵A的第i列上的最大值。 (2) Y,U=max(A)：返回行向量Y和U，Y向量记录A的每列的最大值，U向量记录每列最大值的行号。,MATLAB SDU,25,(3) max(A,dim)：dim取1或2。dim取1时，该函数和max(A)完全相同；dim取2时，该函数返回一个列向量，其第i个元素是A矩阵的第i行上的最大值。
15、求最小值的函数是min，其用法和max完全相同。例 分别求34矩阵x中各列和各行元素中的最大值，并求整个矩阵的最大值和最小值。,MATLAB SDU,26,两个向量或矩阵对应元素的比较,函数max和min还能对两个同型的向量或矩阵进行比较，调用格式为： (1) U=max(A,B)：A,B是两个同型的向量或矩阵，结果U是与A,B同型的向量或矩阵，U的每个元素等于A,B对应元素的较大者。 (2) U=max(A,n)：n是一个标量，结果U是与A同型的向量或矩阵，U的每个元素等于A对应元素和n中的较大者。 min函数的用法和max完全相同。例 求两个23矩阵x, y所有同一位置上的较大元素构成的新
16、矩阵p。,MATLAB SDU,27,求和与求积,数据序列求和与求积的函数是sum和prod，其使用方法类似。设X是一个向量，A是一个矩阵，函数的调用格式为： sum(X)：返回向量X各元素的和。 prod(X)：返回向量X各元素的乘积。 sum(A)：返回一个行向量，其第i个元素是A的第i列的元素和。,MATLAB SDU,28,prod(A)：返回一个行向量，其第i个元素是A的第i列的元素乘积。 sum(A,dim)：当dim为1时，该函数等同于sum(A)；当dim为2时，返回一个列向量，其第i个元素是A的第i行的各元素之和。 prod(A,dim)：当dim为1时，该函数等同于prod
17、(A)；当dim为2时，返回一个列向量，其第i个元素是A的第i行的各元素乘积。例 求矩阵A的每行元素的乘积和全部元素的乘积。,MATLAB SDU,29,平均值和中值,求数据序列平均值的函数是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个元素是
18、A的第i行的算术平均值。 median(A,dim)：当dim为1时，该函数等同于median(A)；当dim为2时，返回一个列向量，其第i个元素是A的第i行的中值。例 分别求向量x与y的平均值和中值。,MATLAB SDU,30,累加和与累乘积 在MATLAB中，使用cumsum和cumprod函数能方便地求得向量和矩阵元素的累加和与累乘积向量，函数的调用格式为： cumsum(X)：返回向量X累加和向量。 cumprod(X)：返回向量X累乘积向量。 cumsum(A)：返回一个矩阵，其第i列是A的第i列的累加和向量。 cumprod(A)：返回一个矩阵，其第i列是A的第i列的累乘积向量。
19、 cumsum(A,dim)：当dim为1时，该函数等同于cumsum(A)；当dim为2时，返回一个矩阵，其第i行是A的第i行的累加和向量。 cumprod(A,dim)：当dim为1时，该函数等同于cumprod(A)；当dim为2时，返回一个向量，其第i行是A的第i行的累乘积向量。 例,MATLAB SDU,31,标准方差与相关系数,求标准方差 在MATLAB中，提供了计算数据序列的标准方差的函数std。对于向量X，std(X)返回一个标准方差。对于矩阵A，std(A)返回一个行向量，它的各个元素便是矩阵A各列或各行的标准方差。std函数的一般调用格式为： Y=std(A,flag,di
20、m) 其中dim取1或2。当dim=1时，求各列元素的标准方差；当dim=2时，则求各行元素的标准方差。flag取0或1，当flag=0时，按公式1(1/n-1)计算标准方差，当flag=1时，按1/n公式计算标准方差。缺省flag=0，dim=1。例对二维矩阵x，从不同维方向求出其标准方差。,MATLAB SDU,32,排序,MATLAB中对向量X是排序函数是sort(X)，函数返回一个对X中的元素按升序排列的新向量。 sort函数也可以对矩阵A的各列或各行重新排序，其调用格式为： Y,I=sort(A,dim) 其中dim指明对A的列还是行进行排序。若dim=1，则按列排；若dim=2，则
21、按行排。Y是排序后的矩阵，而I记录Y中的元素在A中位置。例：help sort,MATLAB SDU,33,协方差与相关系数,MATLAB提供了cov协方差函数赖体现向量中各个元素的分散程度（距阵的各列的相关程度）corrcoef函数，可以求出数据的相关系数矩阵，为相关性的归一化表示。corrcoef函数的调用格式为： corrcoef(X)：返回从矩阵X形成的一个相关系数矩阵。此相关系数矩阵的大小与矩阵X一样。它把矩阵X的每列作为一个变量，然后求它们的相关系数。 corrcoef(X,Y)：在这里，X,Y是向量，它们与corrcoef(X,Y)的作用一样。,MATLAB SDU,34,【补充
22、】曲线拟合,在MATLAB中，用polyfit函数来求得最小二乘拟合多项式的系数，再用polyval函数按所得的多项式计算所给出的点上的函数近似值。 polyfit函数的调用格式为：P,S=polyfit(X,Y,m) 函数根据采样点X和采样点函数值Y，产生一个m次多项式P及其在采样点的误差向量S。其中X,Y是两个等长的向量，P是一个长度为m+1的向量，P的元素为多项式系数。注意：polyval函数的功能是按多项式的系数计算x点多项式的值,MATLAB SDU,35,多项式拟合,多项式拟合又称为曲线拟合，其目的就是在众多的样本点中进行拟合，找出满足样本点分布的多项式 Polyfit(x,y,n
23、) 例：x=1 2 3 4 5 y=3.2 46.7 138.3 278.6 467.2polyfit(x,y,2)polyfit(x,y,3)polyfit(x,y,4),MATLAB SDU,36,【补充】数据插值,一维数据插值在MATLAB中，实现这些插值的函数是interp1，其调用格式为： Y1=interp1(X,Y,X1,method) 函数根据X,Y的值，计算函数在X1处的值。X,Y是两个等长的已知向量，分别描述采样点和样本值，X1是一个向量或标量，描述欲插值的点，Y1是一个与X1等长的插值结果。method是插值方法，允许的取值有linear、nearest、cubic、sp
24、line。,注意：X1的取值范围不能超出X的给定范围，否则，会给出“NaN”错误 Interp1(X,Y,X1,spline)=spline(X,Y,X1),MATLAB SDU,37,例 某观测站测得某日6:00时至18:00时之间每隔2小时的室内外温度()，用3次样条插值分别求得该日室内外6:30至17:30时之间每隔2小时各点的近似温度()。 设时间变量h为一行向量，温度变量t为一个两列矩阵，其中第一列存放室内温度，第二列储存室外温度。命令如下： h =6:2:18; t=18,20,22,25,30,28,24;15,19,24,28,34,32,30; XI =6.5:2:17.5
25、YI=interp1(h,t,XI,spline) %用3次样条插值计算,MATLAB SDU,38,二维数据插值,在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”错误。,MATLAB SDU,39,例
26、 某实验对一根长10米的钢轨进行热源的温度传播测试。用x表示测量点0:2.5:10(米)，用h表示测量时间0:30:60(秒)，用T表示测试所得各点的温度()。试用线性插值求出在一分钟内每隔20秒、钢轨每隔1米处的温度TI。 命令如下： x=0:2.5:10; h=0:30:60; T=95,14,0,0,0;88,48,32,12,6;67,64,54,48,41; xi=0:10; hi=0:20:60; TI=interp2(x,h,T,xi,hi),MATLAB SDU,40,拟合和插值区别,拟合：寻找平滑曲线以最好地表现带噪声的“测量数据”，不要求拟合曲线穿过这些数据点； 插值：研究
27、如何平滑估算出基准数据之间其它点的函数值，所以插值所得曲线必定穿过基准数据。,MATLAB SDU,41,数值差分,MATLAB SDU,42,数值微分,数值微分的实现 在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，按行计算差分。,MAT
28、LAB SDU,43,数值积分,数值积分基本原理求解定积分的数值方法多种多样，如简单的梯形法、辛普生(Simpson)法、牛顿柯特斯(Newton-Cotes)法等都是经常采用的方法。它们的基本思想都是将整个积分区间a,b分成n个子区间xi,xi+1，i=1,2,n，其中x1=a，xn+1=b。这样求定积分问题就分解为求和问题。,MATLAB SDU,44,数值积分的实现方法,低阶法-自适应递推辛普生法 基于变步长辛普生法，MATLAB给出了quad函数来求定积分。该函数的调用格式为：I=quad(fname,a,b,tol,trace)I,n=quad(fname,a,b,tol,trace
29、)其中fname是被积函数名。a和b分别是定积分的下限和上限。tol用来控制积分精度，缺省时取tol=0.001。trace控制是否展现积分过程，若取非0则展现积分过程，取0则不展现，缺省时取trace=0。返回参数I即定积分值，n为被积函数的调用次数。,MATLAB SDU,45,例 求0.3pi定积分 f=exp(-0.5*x)*sin(x+pi/6); 。调用数值积分函数quad求定积分。 S,n=quad(exp(-0.5*x).*sin(x+pi/6),0,3*pi) S =0.9008 n =77,MATLAB SDU,46,2高阶法：自适应牛顿柯特斯法 基于牛顿柯特斯法，MATL
30、AB给出了quadl函数来求定积分。该函数的调用格式为： I,n=quadl(fname,a,b,tol,trace) 其中参数的含义和quad函数相似，只是用高阶自适应递推法，该函数可以更精确地求出定积分的值，且一般情况下函数调用的步数明显小于quad函数，从而保证能以更高的效率求出所需的定积分值。,MATLAB SDU,47,例：前一例子 分别用quad函数和quadl函数求定积分的近似值，并在相同的积分精度下，比较函数的调用次数例： 求0,pi 定积分f=x*sin(x)/(1+cos(x)*cos(x)调用函数quadl求定积分。 I=quadl(x.*sin(x)./(1+cos(x
31、).*cos(x),0,pi) I =2.4674,MATLAB SDU,48,3Trapz : 计算梯形面积的和来计算定积分 在MATLAB中，对由表格形式定义的函数关系的求定积分问题用trapz(X,Y)函数。其中向量X,Y定义函数关系Y=f(X)。例 用trapz函数计算定积分。 命令如下： X=1:0.01:2.5; Y=exp(-X); %生成函数关系数据向量 trapz(X,Y) ans =0.28579682416393,MATLAB SDU,49,7 离散傅里叶变换,离散傅里叶变换是信号分析中的一种重要工具,它将时域内的问题转化为频域内的问题,在很多情况下大大地简化了问题的求解
32、过程;另外计算时域信号的频谱也主要依靠离散傅里叶变换来完成. 离散序列的离散傅里叶变换(DFT)定义为: 在时间片段等距抽取N个抽样时间tm处的样本值为f(tm),m=0,1,N,称向量F(k)为f(m)的一个离散傅里叶变换 由于MATLAB中序列的下标不允许出现0,而是从1开始, 离散序列x(n)的离散傅里叶逆变换定义为:,MATLAB SDU,50,7 离散傅里叶变换（cont.),MATLAB中，提供了对向量(或直接对矩阵的行或列)进行离散傅立叶变换的函数，其调用格式是： Y=fft(X,n,dim) (1)当X是一个向量时，返回对X的离散傅立叶变换。 (2)当X是一个矩阵时，返回一个矩
33、阵并送Y，其列(行)是对X的列(行)的离散傅立叶变换。,MATLAB SDU,51,7 离散傅里叶变换,MATLAB中，提供了对向量(或直接对矩阵的行或列)进行离散傅立叶变换的函数，其调用格式是： Y=fft(X,n,dim) (1)当X是一个向量时，返回对X的离散傅立叶变换。 (2)当X是一个矩阵时，返回一个矩阵并送Y，其列(行)是对X的列(行)的离散傅立叶变换。,MATLAB SDU,52,例求X=(1，0，-3，5，2)的离散傅立叶逆变换。 在MATLAB命令窗口，输入命令： X=1,0,-3,5,2; Y=fft(X) %对X进行变换3. 离散傅立叶变换的逆变换 MATLAB中，对向量
34、(或直接对矩阵的行或列)进行离散傅立叶逆变换的函数的调用方法是： Y=ifft(X,n,dim) 函数对X进行离散傅立叶逆变换。其中X、n、dim的意义及用法和离散傅立叶变换函数fft完全相同。,MATLAB SDU,53,例: 对矩阵A的列向量、行向量分别进行离散傅立叶变换、并对变换结果进行逆变换。 命令如下： A=3,2,1,1;-5,1,0,1;3,2,1,5; fftA=fft(A) %求A的列向量的傅立叶变换 fftA2=fft(A,4,2) %求A的行向量的傅立叶变换 ifft(fftA) %对矩阵fftA的列向量进行傅立叶逆变换，结果应等于A ifft(fftA2,4,2) %对
35、矩阵fftA2的行向量进行傅立叶逆变换，其结果应等于A,MATLAB SDU,54,例：已知模拟信号f(t)= 2*sin(6*pi*t)+7*cos(8*pi*t); 求N=128点DFT的幅度谱和相位谱. MATLAB命令窗口输入如下: N=128; n=0:127; t=0:0.01:127; %时域数据点数 q=n*2*pi/N; %频域数据点数 X=2*sin(6*pi*t)+7*cos(8*pi*t); Y=fft(X,N); subplot(2,1,1); plot(q,abs(Y) subplot(2,1,2); plot(q,angle(Y),MATLAB SDU,55,第5章小结,多项式: 表达、求值、求根、加减乘除微分积分运算 代数方程求解 函数创建，函数图形绘制，函数的极值，零点，积分 数值统计处理：最大最小值，平均值中值，求和求积，累加和累乘，标准方差，排序， 曲线拟合，一维二维数值插值,MATLAB SDU,56,思考题,1 已知 f(x)3x6+4x4-5x3-7.2x2+5(1) 计算f(x)=0 的全部根。(2) 由方程f(x)=0的根构造一个多项式g(x)，并与f(x)进行对比。 2 寻找一个三次多项式在0,2pi内逼近函数sinx,
