matlab数值运算基础PPT课件

-,1,第3章数值运算基础,内容提要：,3.1多项式3.2线性代数3.3数据分析3.4插值3.5数字信号处理初步,-,2,3.1多项式,-,3,3.1.1创建多项式,Matlab多项式的表达,注：系数中的零不能省！,-,4,创建多项式的方法,系数矢量直接输入法特征多项式输入法由根矢量创建多项式,-,5,系数矢量直接输入法,方法:函数poly2sym+系数矢量,例:poly2sym(1234)ans=x3+2*x2+3*x+4,说明:poly2str以习惯方式显示多项式poly2sym双精度多项式系数转为符号多项式,例:poly2str(1234,y)ans=y3+2y2+3y+4,适用于:已知系数表达式,-,6,特征多项式输入法,方法:函数poly+矩阵,poly求矩阵特征多项式、根集对应的多项式,？,适用于:已知矩阵特征多项式,-,7,例:a=123;456;780;p=poly(a)p=1.0000-6.0000-72.0000-27.0000pstr=poly2str(p,x)pstr=x3-6x2-72x-27,-,8,说明:1.N阶方阵特征多项式系数矢量一定是n+1阶的2.特征多项式系数矢量的第一个元素必须为1。,-,9,由根矢量创建多项式,方法:函数poly+根矢量,poly求矩阵特征多项式、根集对应的多项式,适用于:已知根(因式分解形式)多项式,例:已知:(x-1)(x+2),求:系数形式的多项式,p=poly(1-2);poly2sym(p)ans=x2+x-2,-,10,说明:1.实系数多项式的根为实数或者共轭复数。2.当根是复数时,要使用real命令滤除根的虚部对多项式系数的影响(可能使系数为复数)。,-,11,例:r=-0.50.5+0.9i-0.3-0.4i;p=poly(r);pr=real(p);ppr1=poly2sym(p)ppr2=poly2sym(pr)ppr1=x3+3/10*x2-1/2*i*x2+11/100*x-18/25*i*x+21/200-47/200*ippr2=x3+3/10*x2+11/100*x+21/200,-,12,3.1.2多项式运算,求多项式的值求多项式的根多项式的乘除运算多项式的微积分多项式的部分分式展开多项式拟合,-,13,求多项式的值,方法:函数polyval:按数组运算规则求值函数polyvalm:按矩阵运算规则求值,格式:y=polyval(p,x)p为多项式,x可为标量/数组/矩阵y=polyvalm(p,x)x可为标量/方阵,-,14,例:p=111;x=0123;xm=01;23;y1=polyval(p,x)y2=polyval(p,xm)ym=polyvalm(p,xm),y1=13713y2=13713ym=34815,-,15,说明:矩阵运算时多项式常数c为c*单位阵,-,16,求多项式的根,方法:一、直接调用求根函数函数roots（p)二、多项式转化为伴随矩阵，然后求其特征值S=compan(p);r=eig(S),说明:根用列矢量表示,适用于:已知多项式根(因式分解形式),-,17,例:a=1-53-64-10;r=roots(a),r=4.61300.7621+0.9789i0.7621-0.9789i-0.5685+1.0419i-0.5685-1.0419i,-,18,a=1-53-64-10;s=compan(a)r=eig(s),s=5-36-41010000010000010000010r=4.61300.7621+0.9789i0.7621-0.9789i-0.5685+1.0419i-0.5685-1.0419i,-,19,多项式的乘除运算,方法:函数conv（a,b）:多项式乘法、矢量的卷积函数deconv（c,a）:多项式除法法、矢量的卷积的逆函数。其中c为被除多项式。,-,20,a=1-53-42;b=12-53;c=conv(a,b)d=deconv(c,b),c=1-3-1230-3633-226d=1-53-42,-,21,多项式的微积分,方法:函数polyder（p）:多项式的微分函数polyint（p）:多项式的积分,differentialcoefficientintegral,-,22,例:p=3-52-610;pd=polyder(p)poly2sym(pd)pi=polyint(pd)poly2sym(pi),pd=12-154-6ans=12*x3-15*x2+4*x-6pi=3-52-60ans=3*x4-5*x3+2*x2-6*x,-,23,多项式的部分分式展开,方法:函数residue:多项式的部分分式展开及逆过程,-,24,多项式b(x)和不含重根的n阶多项式a(x),当a(x)有m重根p(j)=p(j+1)=p(j+m-1)时,-,25,-,26,a=-4083;b=53-27;r,p,k=residue(b,a),r=-1.4167-0.66531.3320p=1.5737-1.1644-0.4093k=-1.2500,例:求部分分式展开,-,27,b,a=residue(r,p,k),b=-1.2500-0.75000.5000-1.7500a=1.0000-0.0000-2.0000-0.7500,例:部分分式展开逆过程,注意:表达式是分母的归一化处理形式,-,28,多项式拟合,曲线拟合与插值,用一个解析函数描述一组数据(通常是测量值),exp3_13_2.m,-,29,连接数据点的实线描绘了线性内插虚线是数据的曲线拟合,exp3_13_2.m,-,30,多项式拟合,曲线拟合与插值,两种方法:插值法:在已知数据点之间计算估值。以某种方法描述数据之间的情况。曲线拟合:设法找出某条光滑曲线，它最佳地拟合数据，但不必要经过任何数据点。根据数据，求自变量x和因变量y之间的近似函数(又称经验公式或数学模型)Y=F(x)。,用一个解析函数描述一组数据(通常是测量值),-,31,最佳拟合:可以有许多不同的方法定义最佳拟合；最小二乘法:使误差平方和最小。,多项式拟合:近似函数为多项式的拟合;,-,32,多项式拟合,方法:函数polyfit(x,y,n),X：自变量数据数组y：因变量数据数组n：多项式阶数p：拟合多项式系数向量,应用最小二乘法求出n阶拟合多项式,-,33,例exp3_13_2.m:x=(0:0.1:1);y=-.4471.9783.286.167.087.347.669.569.489.3011.2;n=2p=polyfit(x,y,n)xi=linspace(0,1,100);z=polyval(p,xi);plot(x,y,o,x,y,xi,z,:)xlabel(x),ylabel(y=f(x),title(SecondOrderCurveFitting),-,34,-,35,3.2线性代数,-,36,线性方程组系数矩阵三种情况,m=n“恰定”方程mn“超定”方程mn超定方程组,x=inv(AA)AB求逆法x=AB最小二乘法解,两种实现方法:,线性无关方程个数大于未知数个数不存在唯一解,-,40,例:,x1+2x2=12x1+3x2=23x1+4x2=4,A=12;23;34;B=1;2;4;x=ABx=inv(A*A)*A*B,x=2.1667-0.6667x=2.1667-0.6667,=,Ax=B,-,41,mn欠定方程组,x=pinv(A)B具有最小长度或范数的解x=AB具有最多零元素的解,两种实现方法:,线性无关方程个数小于未知数个数无穷多个解,-,42,例:,x1+2x2+3x3=12x1+3x2+4x3=2,A=123;234;B=1;2;x=ABx=pinv(A)*B,x=100 x=0.83330.3333-0.1667,Ax=B,=,-,43,3.3数据分析,-,44,MATLAB约定：,若输入量X为矢量，不论是行还是列矢量，对整个矢量运算，结果为一个数值；若输入量X为矩阵，则函数按列运算，其结果为一个行矢量。,-,45,3.3.1基本统计命令,max求最大元素min求最小元素median求中间数值sort递增排序,mean求均值sum求和cumsum求累计和std求标准差prod求积cumprod求累积,-,46,sum求和cumsum求累计和,A=12-1-223S=sum(A)CS=cumsum(A),S=23CS=120023,Prod与cumprod的区别和上述类似,思考:求X=(1！，2！，3！，10！),X=cumprod(1:10),?,-,47,A=12-1-223,M=23MA=3,思考:求A中的最大值,A=12-1-223M=max(A)MA=max(max(A),-,48,sort递增排序,格式:sort(X),A:矩阵Dim:(=1按列;dim=2按行)Y:排序后的矩阵I:记录Y中的元素在A中位置A矩阵,Y,I=sort(A,dim),按升序排列;x可为数组/矩阵,见例:exp3_19_2.m,-,49,A=1-854126137-13AC=sort(A)AL=sort(A,2)X,I=sort(A)X2,I2=sort(A,2),排序举例:exp3_19_2.m,-,50,回顾:矩阵翻转,格式:fliplr(a)%flipleftrigtflipud(a)%flipupdown,思考:如何降序排列?,-,51,3.3.2协方差阵和相关阵,cov（x，y）求两随机变量的协方差矩阵corrcoef（x，y）求两随机变量的相关系数矩阵,了解,-,52,3.3.3微分、差分与梯度,微分、差分近似梯度,-,53,微分、差分,方法:函数diff,格式:DX=diff(x):计算X的向前差分DX(i)=X(i+1)-X(i)DX=diff(X,n)：计算X的n阶向前差分DX=diff(A,n,dim)：计算矩阵A的n阶差分dim=1时(缺省状态)，按列计算差分dim=2，按行计算差分,-,54,(课后自学),例:exp3_21.m,a=rand(3,3,2)diff(a)diff(a,2)diff(a,3)diff(a,4)diff(a,5)diff(a,6),-,55,近似梯度,方法:函数gradient,格式:FX=gradient(F)计算F一维数值梯度FX,FY=gradient(F)计算F二维数值梯度FX,FY,FZ=gradient(F)计算F三维数值梯度,(课后自学),例:exp3_22.m,-,56,思考:产生一个随机5*5的矩阵；选择12个与几何中心元素差值为最小的元素；用这12个像素的均值替换掉中心元素值。,-,57,3.4插值,-,58,插值分类,插值:在已知数据点之间计算估值。以某种方法描述数据之间的情况。,空间,方法,-,59,3.4.1一维插值,方法:函数yi=interp1(x,y,xi,method),X：采样点数组y：样本值数组xi：插值点数组yi：插值结果method:插值方法,(interpose),-,60,插值方法:nearest-邻近点插值将插点值设为最近数据点值linear-线性插值(默认方法)两采样点间连直线spline-三次样条插值通过数据点拟合出三次样条曲线pchip,cubic-立方插值分段立方插值法计算插值,-,61,选择插值法时考虑的因素:运算时间、占用内存、曲线的光滑程度,nearest：执行速度最快，输出结果为直角转折linear：在样本点上斜率变化很大spline：最花时间，但输出结果也最平滑cubic：最占内存，输出结果与spline差不多,-,62,例exp3_23c.m:一维插值函数插值方法的对比x=0:10;y=sin(x);xi=0:.25:10;strmod=nearest,linear,spline,cubicstrlb=(a)nearest,(b)linear,.(c)spline,(d)cubic;fori=1:4yi=interp1(x,y,xi,strmodi);subplot(2,2,i),plot(x,y,ro,xi,yi,b).,xlabel(strlb(i)end,-,63,-,64,MATLAB中有一个专门的三次样条插值函数y1=spline(x,y,xi)等同于yi=interp1(x,y,xi,spline),-,65,3.4.2二维插值,方法:函数yi=interp2(x,y,z,xi,yi,method),3.4.3三维插值,方法:函数yi=interp3(x,y,z,w,xi,yi,zi,method),插值方法同一维插值(自学),-,66,3.5数字信号处理初步,波形的产生离散傅里叶变换,说明：本部分需有数字信号处理知识基础,此处做大概讲解,-,67,3.5.1波形的产生,y=square（2*pi*t）周期为1、峰值为1的方波y=square（2*pi*t,duty）duty为占空比,1.方波,-,68,例exp3_27_2.m:产生周期为0.5的标准方波和占空比为80%的方波。,fs=10000;t=0:1/fs:1;x=square(2*pi*2*t);y=square(2*pi*2*t,80);figure(1);plot(t,x);axis(01-1.51.5);figure(2);plot(t,y);axis(01-1.51.5);,-,69,y=sawtooth（2*pi*t）周期为1、峰值为1的锯齿波y=sawtooth（2*pi*t,width）不同形状的三角波width指定y最大值出现的位置，是位置横坐标与周期的比值。,2.锯齿波或三角波,-,70,例exp3_27_3.m:产生周期为2的锯齿波和对称三角波,fs=10000;t=0:1/fs:4;x=sawtooth(2*pi*0.5*t);y=sawtooth(2*pi*0.5*t,0.5);figure(1);plot(t,x);axis(04-1.51.5);figure(2);plot(t,y);axis(04-1.51.5);,-,71,3.5.2快速傅立叶变换,y=fft(x)用快速傅里叶变换计算离散傅里叶变换y=fft(x,n)用快速傅里叶变换计算n点的离散傅里叶变换,-,72,例exp3_27_4.m:产生一个含有两个不同频率的正弦信号，进行快速傅里叶变换，并画出幅值图、相位图,tp=0:2048;%时域数据点数Nyt=0.02*sin(0.03*pi*tp)+0