第4章 数值运算基础MATLAB的计算包括数值计算 和符号计算两类，本章将带大家学 习数值计算部分。其中将主要学习 与我们专业密切相关的多项式、方 程组求解、数据分析和数字信号处 理的快速傅里叶变换。1物理与电气工程学院第1节 多项式polynomial MATLAB用行向量表示多项式。将 多项式的系数按降幂次序存放在行向量 中。如：的系数行向量为 注意：多项式中缺少的幂次系数一 定要用“0”补齐。 2物理与电气工程学院一、创建多项式 1、系数矢量直接输入法在命令窗口直接输入多式的系数向量【例4-1】输入系数矢量，创建多项式x3-4*x2+3*x+2。 p=1 -4 3 2 poly2sym(p)%将矢量P表示为多项式的手写形式Polynomial coefficient vector to symbolic polynominal3物理与电气工程学院2、方阵特征多项式输入法 p=poly(A) 若A为nn的矩阵，则返回值P将是一个 含有n+1个元素的行向量，也就是该矩阵特 征多项式的系数 【例4-2】求矩阵1 2 3;4 5 6;7 8 0的特征多项 式系数，并转换为多项式形式。 a=1 2 3;4 5 6;7 8 0; p=poly(a); poly2sym(p) %将矢量P表示为多项式的手写形式 d1=roots(p) %由特征多项式求得的特征值 d2=eig(a) %由特征值函数求得的特征值43、根矢量创建p=poly(A) A为待求多项式的根 矢量，则返回值将 是对应多项式的系 数行矢量，该多项 式的根为矢量A。此 时p=poly(A)与 A=roots(p)互逆。系 统定义P0=1。Ax1 x2 x3p=1 p1 p2 p3 5【例4-3】根据根矢量-0.5 -0.3+0.4i -0.3-0.4i，创建多项式 r=-0.5 -0.3+0.4i -0.3-0.4i; p=poly(r) pr=real(p) ppr=poly2sym(pr) 二、多项式运算 1、求多项式的值 MATLAB的多项式求值方式有两种，按数组运算 规则和按矩阵运算规则计算多项式的值。y=polyval(p,x) 按数组规则运算，计算多项式 在x处的值，p是多项式的系数矢量；x是指定自 变量的值，可以是标量、向量或矩阵。如果x是 向量或矩阵，则函数的返回值是与x对应的同维 向量或矩阵。 6【例4-4】求多项式3x2+2x+1在5、7和9处的值 。 p = 3 2 1; polyval(p,5 7 9)y=polyvalm(p,x) 将矩阵整体（而不是矩阵元素 ）作为自变量进行计算。p是多项式的系数向 量；相当于用矩阵x代替多项式的变量对矩阵 而不是对数组进行计算，x必须是方阵【例4-5】求多项式3x2+2x+1对于矩阵2 5;7 9的值p= 3 2 1; polyvalm(p,2 5;7 9)7物理与电气工程学院2、求多项式的根格式：C=roots(p) p为多项式的系数矢量，C为函数返回多项式的根矢量 如果C为复数，则必成对出现。【例4-6】分别用两种方法求多项式x5-5x4+3x3-6x2+4x-10的根。 a=1 -5 3 -6 4 -10; r=roots(a)8物理与电气工程学院3、多项式的乘除运算 多项式的乘法conv格式： c=conv(a,b) 多项式的乘法运算,也是矢量的卷积运算 向量a长度为m，向量b长度为n，a和b的卷积 定义为：运算结果矢量c为长度k=m+n-1 【例4-7】计算两多项式x4-5x3+3x2-4x+2和 x3+2x2-5x+3的乘法a=1 -5 3 -4 2;b=1 2 -5 3;c=conv(a,b)poly2sym(c)9物理与电气工程学院 多项式的除法deconv格式： q, r=deconv( c,a)多项式的除法运算，也是矢量的解卷积运算 过程。向量a对向量c进行解卷积得到“商”向量q和“ 余”向量r。q和r分别是商多项式和余多项式；c和a分别 是被除多项式和除多项式，使得：c=conv(a,q)+r【例4-8】计算例4-7中求得的乘积被x3+2x2-5x+3除 所得结果c=1 -3 -12 30 -36 33 -22 6;b=1 2 -5 3;q ,r=deconv(c,b)104、多项式微积分 polyder(p) 返回多项式系数向量p 的导数【例4-9】计算多项式3x4-5x3+2x2-6x+10的微分 。p=3 -5 2 -6 10;dp=polyder(p)poly2sym(dp) polyint(p) 返回多项式系数p 的积分【例4-10】计算多项式12x3-15x2+4x-6的积分 。p=12 -15 4 -6;ip=polyint(p)poly2sym(ip)11物理与电气工程学院5、多项式的部分分式展开 MATLAB提供了residue命令来执行部分分式 展开或多项式系数之间的转换。该命令通 常用于信号与控制领域中。格式如下： r,p,k=residue(b,a) 该命令是求多项式之比b(s)/a(s)的部分分 式展开，返回留数r、极点p和直项向量k。 a和b分别是分母和分子多项式的系数向量 b,a=residue(r,p,k) 上一条语句的逆过程，只是分母多项式系 数以归一化形式：最高次项系数a(1)为1。12物理与电气工程学院 如果分母多项式a(s)不含重根，则两个多项式 可写成如下形式：其中，pi称为极点，ri称为留数，k(s)是直项 。 如果b的次数低于a的次数，则直项为空。 如果分母多项式a(s)含m重根p(j)=p(j+m-1),则 这m项应展开为13b(x) 5x3+3x2-2x+7 【例3-11】两多项式的比为 = ,a(x) -4x3+8x+3求部分分式展开。 a = -4 0 8 3; b = 5 3 -2 7; r, p, k = residue(b,a) b1,a1 = residue(r,p,k) %分母最高次项归1 r2,p2,k2 = residue(1 1,1 -2 1) %出现重根笔算结果=？14物理与电气工程学院6、多项式拟合对于给定的一组数据(xi ,yi)，i=1,2,n,如果 要采用多项式模型对数据组进行描述，形成如多 项式y(x)=f(x,p)=p1 xn+ p2 xn-1+ p3 xn-2+ pn+1的形式，求取参数p使得量值2(p)的值最小的过程 ，称为对数据组进行多项式拟合，其中MATLAB系统设计polyfit函数采用最小二乘 法原理对给定的数组(xi,yi)，i=1,2,n进行 多项式的曲线拟合，最后给出拟合的多项式 系数。15物理与电气工程学院p=polyfit(x,y,n) 其中，x,y分别表示横、纵 坐标向量；n是给定的拟合多项式的最高阶数 ，返回一个多项式系数向量p。如n=3,若p=1 0.5 1 2,则 y=1*x3+0.5*x2+1*x1+2*x0p,S=polyfit(x,y,n) 返回n阶多项式系数p, S为误差估计结构 p,S,mu=polyfit(x,y,n) 对归一化以后的数据 进行多项式拟合，用x x=(x-u1)/u2替代x，其中 mu=u1 u2,u1=mean(x),u2=std(x)16物理与电气工程学院【例3-13】 求误差函数的6阶拟合多项式。 x = (0: 0.1: 2.5); % 生成0至2.5间隔为0.1的自变量 y = erf(x); % 计算误差函数 p = polyfit(x,y,6) % 求6阶拟合多项式 x = (0: 0.1: 5); % 生成0至5间隔为0.1的自变量 y1 = erf(x); % 计算误差函数 f = polyval(p,x); % 计算拟合函数的值 plot(x,y1,o,x,f,-) % 绘图函数 p0,s0,mu0 = polyfit(x,y,6) %x=(x-mean(x)/std(x) p1,s1,mu1 = polyfit(x-mu(1)/mu(2),y,6)可以看出，拟合区间0 2.5内拟合曲线拟近原 曲线，而在区间以外的曲线误差较大17物理与电气工程学院第2节 线性代数 给定两个矩阵A和B,求X的解,使得:AX=B XA=B在MATLAB中，求解线性方程组时，主要采 用前面章节介绍的除法运算符“”和“/”求解： X=AB 是AX=B 的解 X=B/A 是XA=B 的解对于方程AX=B，要求矩阵A和B有相同的行数； X和B有相同的列数，且它们的行数等于A的列数。18物理与电气工程学院根据矩阵A的结构(m, n），可以将方 程分为以下3类： m=n 方阵系统，可偿试求精确解 mn 超定系统，可偿试求最小二乘解 mn,则方程没有精确解,此时方程组为超定 方程组。一般采用最小二乘法。 左除法：x=Ab 建立在奇异值分解基 础之上 广义逆法：x=pinv(A)*b 速度较快， 可靠性较差一些实验中数据的曲线拟合就是就可以解超定方程组 的方法来解。一般情况下需要将非线性问题转换 为线性问题来解决 22物理与电气工程学院一组实验数据，时间t和测量数据y，如下表所示：itexp(-ti)yi10.010.8220.30.7400.7230.80.4490.6341.10.3320.6051.60.2020.5562.20.1110.5认为x和y有上式的关系式 ，则由6组实验数据就可形 成超定方程组，就可求解 出C1和C2。使得将t代入函 数得到的值与实际y值之间 差值的最小平方和最小如何以矩阵的 形式表示？23物理与电气工程学院【例4-16】求表中数据的最小二乘解。 t=0 0.3 0.8 1.1 1.6 2.3; y=0.82 0.72 0.63 0.60 0.55 0.50; e=ones(size(t) exp(-t) %6组实验数据变换得到的 c=ey cp=pinv(e)*y t1=0:0.1:2.5; y1=ones(size(t1),exp(-t1)*c; yp=ones(size(t1),exp(-t1)*cp; plot(t,y,t1,y1,ro,t1,yp,b*)24物理与电气工程学院三 、欠定系统程组Ax=b，A为mn矩阵,如果mc=0.612 c=norminv(0.612,3,2) p=normcdf(c,3,2) %用于验证33物理与电气工程学院二、协方差矩阵和相关阵 C=cov(x) 如果x是向量，则返回C是该向量的方差 ,C是一个标量； C=cov(x)如果x是矩阵，则每一列相当于一个变量，返回C是 该矩阵的列与列之间的协方差矩阵：C是一个列数与x相同的 对称方阵，主对角线元素是x对应列矢量的方差 c(i,j)=meanx(:,i)-mean(x(:,i)x(:,j)-mean(x(:,j) C=cov(x,y) x和y是同维向量或矩阵，相当于cov(x(:) y(:), 既先将x和y按展开成列向量,再求这两个列向量的协方差diag(cov(x) 如果x是矩阵，则C就是该矩阵每个列向量的方 差34S=corrcoef(x) 计算x的相关系数矩阵S， 将x的每一列视为一个变量，S(i , j)是i列 与j列之间的相关