MATLAB5-1数值计算.ppt

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