MATLAB语言及应用3_zeng

MATLAB语言及应用──PART3数值计算曾然杭州电子科技大学通信工程学院数值计算内容:数值计算包括:多项式运算,拟合与插值，数值微积分，线性方程组求解,数据分析,信号处理和系统分析.1多项式运算多项式运算是数学中最基本的运算之一.多项式一般可以表示为:f(x)=a0xn+a1xn-1+a2xn-2+…+an-1x+an对于这种表示形式，很容易用一个行向量表示，即：T=[a0,a1,a2,…an-1,an]在MATLAB中,多项式正是用这样一个行向量表示的,系数是按降序排列的．幂多项式构造多项式可以直接用向量表示,因此,构造多项式最简单的方法是直接输入向量.【例】直接输入向量构造f(x)=2x5+5x4+4x2+x+4T=[2,5,0,4,1,4];fx=poly2sym(T)函数poly2sym是符号工具箱中的函数,在用此种方式构造多项式时,必须把多项式各项的系数写完整,而不管此项的系数是否为0fx=2*x^5+5*x^4+4*x^2+x+4r=[1,2,3,4];T1=poly(r);y=poly2sym(T1)y_class=class(y)%显示出y为符号变量【例】用多项式的根构造多项式,根为r=[1,2,3,4]T1=1-1035-5024y=x^4-10*x^3+35*x^2-50*x+24y_class=sym多项式的运算多项式的运算主要包括多项式的四则运算,导数运算,估值运算,求根运算以及多项式的拟合、插值等．1多项式的四则运算多项式四则运算主要是多项式的加,减,乘,除.其中,多项式的加减运算要求两个相加,减的多项式的大小必须相等,此时加,减才有效.当两个相加,减的多项式阶次不同时必须通过补0使其相同.加/减---+/-乘/除---conv,deconv[T,r]=deconv(T1,T3)商多项式余式T3为分母T1=[2,5,0,4,1,4];T2=[0,0,5,1,3,2];T3=[5,1,3,2];%重写T2，因除法运算中分母多项式第一个系数不能为0（为了与T1做加法运算的T2前两个系数补了0，不能再作除法运算）T=T1+T2;%必须是同维的才能相加T_add=poly2sym(T)T=T1-T2;T_sub=poly2sym(T)%注意：加减运算只要对系数操作即可，但乘除运算需要用conv和deconv指令T=conv(T1,T2);%乘法不要求同维T_mul=poly2sym(T)[A_coe,A_r]=deconv(T1,T3);T_coe=poly2sym(A_coe)T_rem=poly2sym(A_r)【例】多项式的加减乘除运算f1(x)=2x5+5x4+4x2+x+4,f2(x)=5x3+x2+3x+2【例】多项式求值,求上式f1(x)在x=0.5处的函数值T1=[2,5,0,4,1,4];x=0.5;y=polyval(T1,x)y=5.8750Polyval(p,X)求多项式的值指令：当X为输入数组时，则对应求出每个元素作为多项式自变量求出的各个函数值作为输出数组3多项式的导数运算---polyder【例】求多项式f1(x)=2x5+5x4+4x2+x+4的导数T1=[2,5,0,4,1,4];h=polyder(T1);poly2sym(h)2多项式求根---roots【例】求多项式f1(x)=2x5+5x4+4x2+x+4的根T1=[2,5,0,4,1,4];root=roots(T1)root=-2.77090.5611+0.7840i0.5611-0.7840i-0.4257+0.7716i-0.4257-0.7716i10*x^4+20*x^3+8*x+14拟合和插值---polyfit,interp1p=polyfit(x0,y0,n);[p,s]=polyfit(x0,y0,n);x0,y0---给定数据对n----拟合出的多项式次数p---多项式向量s---偏差信息yi=interp1(x0,y0,xi,'cubic');xi,yi---得到的新的插值点cubic---插值方式插值方式:linear---线性插值nearest---最近插值cubic---三次多项式插值spline---样条插值拟合和插值的区别:曲线拟合是研究如何寻找“平滑”曲线最好地表现带噪声的“测量数据”,但并不要求拟合曲线通过这些“测试数据”点。插值是在认定所给数据完全正确的情况下,研究如何平滑地估算出“基准数据”之间其他点的函数值,因此插值一定通过“基准数据”。【拟合的例子】由下列数据对(x0,y0)求三次和六次拟合多项式并绘图.x0=0:0.1:1;y0=[-0.447,1.978,3.11,5.25,5.02,4.66,4.01,4.58,3.45,5.35,9.22];通俗来讲,拟合就是由已知点得到一条曲线,使该曲线最能反映点所代表的规律.比如做欧姆定理的实验的时候,由于实验中存在误差,最后拟合得到的曲线是一条直线,而且肯定只有部分点落在拟合的直线上,但此时该直线和测试点的方差最小.由拟合直线的斜率就可以知道电阻的阻值.拟合是探测事物变化规律的办法.插值就是根据函数上某些已知点(或实验数据),按一定规律(插值方法)寻求未知的点,比如已知一个常用对数y=log(x)表,是按照x=0.1:0.1:10制表的,如果按已知数据求y=log(2.897)就可以用插值得到.表制得越密,插值越准确．x0=0:0.1:1;y0=[-0.447,1.978,3.11,5.25,5.02,4.66,4.01,4.58,3.45,5.35,9.22];n=3;%设定拟合次数为3[p,s]=polyfit(x0,y0,n);%得到拟合多项式向量和相关偏差信息xx=0:0.01:1;yy=polyval(p,xx);%按拟合曲线计算采样值，以便画图n1=6;%设定拟合次数为6[p1,s1]=polyfit(x0,y0,n1);yy1=polyval(p1,xx);plot(xx,yy,'-b',xx,yy1,'-m',x0,y0,'.r','MarkerSize',20);y3=polyval(p,0.5);y6=polyval(p1,0.5);text(0.5,y3-0.3,'n=3');text(0.5,y6+0.2,'n=6');注意：拟合只能在给定数据所限定的区间内使用,不要任意向外拓展．y=51.48x3-77.74x2+35.06x-0.20y=348.21x6-1060.23x5+1297.68x4-758.90x3+181.00x2+1.00x+0.48【例】按上例所给数据研究插值,并观察插值和拟合的区别.x0=0:0.1:1;y0=[-0.447,1.978,3.11,5.25,5.02,4.66,4.01,4.58,3.45,5.35,9.22];xi=0:0.02:1;yi=interp1(x0,y0,xi,'cubic');subplot(1,2,1);plot(xi,yi,'-b',x0,y0,'.r','MarkerSize',20);xlabel('x'),ylabel('y');[p,s]=polyfit(x0,y0,3);xx=0:0.01:1;yy=polyval(p,xx);subplot(1,2,2);plot(xx,yy,'-r',x0,y0,'.r','MarkerSize',20);xlabel('x');ylabel('y');插值拟合2数据分析函数(1)若输入宗量x是向量,那么不管是行向量还是列向量,运算是对整个向量进行的．(2)若输入宗量x是二维数组，那么指令运算是按列进行的．见下列例子MATLAB在进行数据分析时的约定:(1)随机数发生器和统计分析指令函数含义rand(m,n)产生mxn维的[0,1]区间均匀分布随机数组．randn(m,n)产生mxn维的N(0,1)正态分布随机数组(均值为0标准差为1)（设置随机数发生器参数，见前章）min(x)对mxn矩阵各列分别求最小值max(x)对mxn矩阵各列分别求最大值median(x)对mxn矩阵各列分别求中位数mean(x)对mxn矩阵各列分别求均值std(x)对mxn矩阵各列分别求标准差var(x)对矩阵各列分别求方差（标准差的平方）11mjijixxm122111mjijjisxxm【例】基本统计示例randn('state',0);A=randn(1000,4);AMAX=max(A),AMIN=min(A)%在+3,-3附近AMED=median(A),AMEAN=mean(A)%接近0ASTD=std(A)%接近1已知A=[1,2,3;4,5,6;7,8,9],演示max,min,median,mean,std,var各列最大,最小值：max(A),min(A)各行最大,最小值：max(A'),min(A')各列中间值：median(A)各列平均值：mean(A)各列方差,标准差：var(A),std(A)A=[1,2,3;4,5,6;7,8,9]789(2)差分和累计指令函数含义diff(X,m,n)沿第n维求m次差分(二维：n=1;按列,n=2,按行)[DZx,DZy]=gradient(z)对Z求x,y方向(行、列)梯度prod(X,n)沿第n维求乘积sum(X,n)沿第n维求和(二维：n=1;按列,n=2,按行)trapz(x,y)用梯形法沿列方向求函数y关于自变量x的积分cumprod(X,n)沿第n维求累计乘积cumsum(X,n)沿第n维求累计和cumtrapz(x,Y,n)用梯形法沿第n维求函数Y关于x的累计积分[XS,KK]=sort(X,n)沿第n维对X元素按模增大排列(2)差分和累计指令注：diff:当X是向量时，dx=X(2:n)-X(1:n-1)；当X是矩阵时，dx=X(2:n,:)-(1:n-1,:)即按列求差分。注意diff(A)的长度比X的长度短少一个元素。例：A=[142;369;325],diff(A),diff(A,1,1),diff(A,1,2)gradient:输入z为数组时，DZx,DZy给出与z同样大小的数组，DZx的每行给出z相应行元素之间的“梯度”，即如第一行DZx(1,1)=Z(1,2)-Z(1,1),DZx(1,end)=Z(1,end)-Z(1,end-1),DZx(1,2:end)=1/2*[Z(1,3:end)-Z(1,1:end-2)]；DZy的每列给出z相应列元素之间的“梯度”。例：Z=[1254;3245;31056;4863],[DZx,DZy]=gradient(Z)cumsum与cumtrapz的区别（矩形法、梯形法）trapz(x,y)给出采样点(x,y)所连接折线下的积分面积（为单值）；S=cumtrapz(x,y)为累计积分，即输出为一个与y大小相同的数组，每个元素S(k)为对应该位置的自变量区间[x(1)~x(k)]上的函数y积分值diff,trapz,cumtrapz指令单参数输入的演示a=[1,2,3,4,5];b=[a',a'+1,a'+2,a'+3]diff(a)=diff(a,1)=[1,1,1,1]diff(b)=[1,1,1,1;1,1,1,1;1,1,1,1;1,1,1,1]trapz(b)=[12162024]cumtrapz(b)00001.502.503.504.504.006.008.0010.007.5010.5013.5016.5012.0016.0020.0024.0012342345345645675678【例】(非累计的指令sum和prod:给出单值而非数组)求1+2+3+…+100以及50!x=1:1:100;sum_x=sum(x);sum_x=5050a=1:1:50;prod_a=prod(a)prod_a=3.0414e+064cumsum,cumprod的用法分别是求向量的累计和与累计乘积.如果a=