1、第6章数据分析与多项式计算【本章学习目标】掌握数据统计和分析的方法。掌握数值插值与曲线拟合的方法及其应用。了解快速傅立叶变换的应用。掌握多项式的常用运算6.1 数据统计处理1 1求最大值和最小值求最大值和最小值求一个数据序列最大值的函数max，其调用格式如下。Y =max(X)：若X X是向量，则返回向量X X的最大值；若X X是矩阵，则返回一个包含每一列最大值的行向量。Y,U=max(X,dim)：沿维度dim返回最大值。求矩阵最小值的函数是min，其用法和max相同。6.1.1 最大值和最小值6.1 数据统计处理【例6.1】分别求矩阵 中各列和各行元素中的最大值。x=54,86,453,4
2、5;90,32,64,54;-23,12,71,18;y=max(x)%求矩阵x中各列元素的最大值y=90 86 453 54 y=max(x,2)%求矩阵x中各行元素的最大值y=453 90 7154864534590326454231271186.1 数据统计处理2 2两个向量或矩阵对应元素的比较两个向量或矩阵对应元素的比较函数max和min还能对两个同型的向量或矩阵进行比较。max函数调用格式如下。max(X,Y)：返回值是与X X、Y Y同型的向量或矩阵，其中的每个元素等于X X、Y Y对应元素的较大者。max(X,n)：n是一个标量，返回值是与X X同型的向量或矩阵，其中的每个元素等
3、于X X对应元素和n中的较大者。6.1.1 最大值和最小值6.1 数据统计处理【例6.2】已知 ，求矩阵x、y所有同一位置上的较大元素构成的新矩阵p。x=443,45,43;67,34,-43;y=65,73,34;61,84,326;p=max(x,y)p=443 73 43 67 84 326 f=45;p=max(x,f)p=443 45 45 67 45 454434543673443x6573346184326y6.1 数据统计处理sum函数是用于对数据序列求和。调用格式如下。S=sum(X)：如果X X是一个向量，则返回向量各元素的和。如果X X是一个矩阵，则返回一个行向量，其第i
4、个元素是X X的第i列的元素和。S=sum(X,dim)：当dim为1（默认值）时，该函数等同于sum(X)；当dim为2时，返回一个列向量，其第i个元素是X X的第i行的各元素之和。数据序列求积的函数是prod6.1.2 求和与求积6.1 数据统计处理【例6.3】已知求矩阵A的每行元素之和和全部元素之和。A=9,10,11,12;100,200,300,400;50,60,50,60;S=sum(A,2)%求A的每行元素之和S=42 1000 220 P=sum(S)%求A的全部元素之和P=12626.1 数据统计处理求数据序列平均值的函数是mean，调用格式如下。M=mean(X)：如果X
5、是一个向量，则返回向量的算术平均值。如果X是一个矩阵，则返回一个行向量，其第i个元素是X的第i列的算术平均值。M=mean(X,dim)：当dim为1（默认值）时，该函数等同于mean(X)；当dim为2时，返回一个列向量，其第i个元素是X的第i行的算术平均值。求数据序列中值的函数是median6.1.3 平均值和中值6.1 数据统计处理例如，求向量x=-8，2，4，7，9与y=-8，2，4，7，9，15的平均值和中值。x=-8,2,4,7,9;%奇数个元素 mx=mean(x),median(x)mx=2.8000 4.0000 y=-8,2,4,7,9,15;%偶数个元素 my=mean(
6、y),median(y)my=4.8333 5.50006.1 数据统计处理求累加和的cumsum函数与求累乘积的cumprod函数。cumsum函数的调用格式如下。B=cumsum(X)：如果X X是一个向量，则返回累加和向量。如果X X是一个矩阵，返回的矩阵的第i列是X X的第i列的累加和向量。B=cumsum(X,dim)：返回多维数组的累加和。若X是矩阵，当dim为1（默认值）时，返回的矩阵的第i列是X X的第i列的累加和向量；当dim为2时，返回的矩阵的第i行是X X的第i行的累加和向量。累乘积的cumprod函数用法和cumsum函数相同。6.1.4 累加和与累乘积6.1 数据统计
7、处理 y=cumsum(1:10)y=1 3 6 10 15 21 28 36 45 55 s=sum(y)s=2206.1 数据统计处理1标准差对于向量X X，std(X)返回一个标量。对于矩阵X X，std(X)返回一个行向量，它的各个元素便是矩阵X X各列或各行的标准差。std函数的调用格式为s=std(X,w,dim)选项w用于指定标准差的计算方法，当w 为0（默认值）时，按1所列公式计算标准差；当w为1时，按2所列公式计算标准差。选项dim 为1（默认值）时，求各列元素的标准差；dim 为 2时，求各行元素的标准差。6.1.5 统计描述函数6.1 数据统计处理【例6.56.5】某次射
8、击选拔比赛中小明与小华的10次射击成绩（单位：环）如表6.1所示，试比较两人的成绩。hitmark=7,4,9,8,10,7,8,7,8,7;7,6,10,5,9,8,10,9,5,6;mean(hitmark,2)ans=7.5000 7.5000 std(hitmark,2)ans=1.5811 1.95796.1 数据统计处理2方差var函数的调用格式为V=var(X,w,dim)对于向量X X，var(X)返回一个标量。对于矩阵X X，var(X)返回一个行向量，各个元素是矩阵X X各列或各行的方差。其中，选项w用于指定权重方案。若w为0（默认值），则按观测值数量-1（即k=n-1）实
9、现归一化；若w为1，则按观测值数量（即k=n）实现归一化。选项dim指定沿维度dim计算方差。对于矩阵，当dim 为1（默认值）时，求各列元素的方差；当dim 为 2时，求各行元素的方差。6.1.5 统计描述函数6.1 数据统计处理【例6.6】考察一台机器的产品质量，判定机器工作是否正常。根据该行业通用法则：如果一个样本中的14个数据项的方差大于0.005，则该机器必须关闭待修。假设搜集的数据如表6.2所示，问此时的机器是否必须关闭？samples=3.43,3.45,3.43,3.48,3.52,3.50,3.39,3.48,3.41,3.38,3.49,.3.45,3.51,3.50;va
10、r_samples=var(samples)var_samples=0.00216.1 数据统计处理3相关系数corrcoef函数计算数据的相关系数。调用格式如下。R,P=corrcoef(X,Y)：返回相关系数矩阵和p值矩阵。如果得到的p值矩阵的非对角线元素小于显著性水平（即90%置信区间，默认为 0.05），则R中的相应相关性被视为显著。R,P=corrcoef(X)：返回矩阵X X各列的相关系数，计算时把矩阵X X的每列作为一个观测变量，然后求各列的相关系数。6.1.5 统计描述函数6.1 数据统计处理【例6.7】随机抽取15名健康成人，测定血液的凝血酶浓度及凝血时间，数据如表6.3所示
11、。分析凝血酶浓度与凝血时间之间的相关性。density=1.1,1.2,1.0,0.9,1.2,1.1,0.9,0.6,1.0,0.9,1.1,0.9,1.1,1,0.7;cruortime=14,13,15,15,13,14,16,17,14,16,15,16,14,15,17;R=corrcoef(density,cruortime)R=1.0000 -0.9265 -0.9265 1.00006.1 数据统计处理4协方差cov函数用于计算数据序列的协方差。调用格式如下。C=cov(x)：若x是向量，则返回x 的方差；若x是矩阵，则x的每一行代表一个样本，每一列代表一个观测变量。C=cov
12、(x,y)：返回变量x和y的协方差，x与y同样大小。如果x和y是矩阵，则cov(x,y)将x和y视为列向量，等价于 cov(x(:),y(:)。如果两个变量的协方差是正值，说明两者是正相关的，即两个变量的变化趋势一致；如果协方差为负值，则说明两者是负相关的，即两个变量的变化趋势相反；如果协方差为0，说明两者之间没有关系。6.1.5 统计描述函数6.1 数据统计处理【例6.7】随机抽取15名健康成人，测定血液的凝血酶浓度及凝血时间，数据如表6.3所示。分析凝血酶浓度与凝血时间之间的相关性。计算例6.7中两组数据的协方差，命令如下：density=1.1,1.2,1.0,0.9,1.2,1.1,0
13、.9,0.6,1.0,0.9,1.1,0.9,1.1,1,0.7;cruortime=14,13,15,15,13,14,16,17,14,16,15,16,14,15,17;C=cov(density,cruortime)C=0.0289 -0.2014 -0.2014 1.63816.1 数据统计处理对数组元素进行排序的函数sort(X)，调用格式为Y,I=sort(X,dim,mode)其中Y Y是排序后的矩阵，而I I记录Y Y中的元素在X X中的位置。选项dim指定排序的维度，若dim 为1（默认值），则按列排序；若dim=2，则按行排序。选项mode指明排序的方法，ascend（默
14、认值）为升序，descend为降序。6.1.6 排序6.1 数据统计处理【例6.8】对二维矩阵 做各种排序。A=1,-8,5;4,12,6;13,7,-13;Y=sort(A,2,descend)%对A的每行按降序排序Y=5 1 -8 12 6 4 13 7 -13 X,I=sort(A)%对A的每列按升序排序，矩阵I存储X各元素在A对应列中的行号X=1 -8 -13 4 7 5 13 12 6I=1 1 3 2 3 1 3 2 21854126137136.2 多项式计算1多项式的加减运算在MATLAB中，多项式的加减运算就是其所对应的系数向量的加减运算。对于次数相同的两个多项式，可直接对多
15、项式系数向量进行加减运算。如果多项式的次数不同，则应该把低次多项式系数不足的高次项用0补足，使参与运算的各多项式具有相同的次数。6.2.1 多项式的四则运算6.2 多项式计算例如，计算(x3 2x2+5x+3)+(6x 1)，对于和式的后一个多项式6x 1，它仅为1次多项式，而前面的是3次。为确保两者次数相同，应把后者的系数向量处理成0,0,6,1。命令如下：a=1,-2,5,3;b=0,0,6,-1;c=a+bc=1 -2 11 26.2 多项式计算2 2多项式的乘除运算多项式的乘除运算计算多项式乘积的conv函数和对多项式作除法运算的deconv函数。函数的用法如下。w=conv(P1,P
16、2)Q,r=deconv(P1,P2)P P1、P P2是两个多项式的系数向量w是两个多项式相乘所得多项式的系数向量Q Q是商式的系数向量r r是余式的系数向量。6.2.1 多项式的四则运算6.2 多项式计算【例6.9】求 和 A=1,8,0,0,-10;B=2,-1,3;C=conv(A,B)C=2 15 -5 24 -20 10 -30 P,r=deconv(A,B)P=0.5000 4.2500 1.3750r=0 0 0 -11.3750 -14.1250以下命令验证deconv和conv是互逆的。conv(B,P)+rans=1 8 0 0 -10432(810)(23)xxxx43
17、281023xxxx6.2 多项式计算k=k=polyderpolyder(P)(P)：求多项式：求多项式P P的导数，即的导数，即 k=k=polyderpolyder(P,Q)(P,Q)：求：求PQPQ的导数，即的导数，即 q,dq,d=polyderpolyder(P,Q)(P,Q)：求：求P/QP/Q的导数，的导数，即即6.2.2 多项式的求导)()(xPdxdxk)()()(xQxPdxdxk)()()()(qxQxPdxdxdx6.2 多项式计算【例6.10】求有理分式 的导数。P=1;Q=1,0,5;p,q=polyder(P,Q)p=-2 0q=1 0 10 0 2521()5
18、f xx6.2 多项式计算1 1代数多项式求值代数多项式求值polyvalpolyval函数用来求代数多项式的值，其调用格式为函数用来求代数多项式的值，其调用格式为y=y=polyvalpolyval(p,xp,x)p p是多项式系数向量是多项式系数向量。若若x x为标量，得到多项式在该点的值；若为标量，得到多项式在该点的值；若x x为向量或矩阵，则为向量或矩阵，则对向量或矩阵中的每个元素求多项式的值。对向量或矩阵中的每个元素求多项式的值。6.2.3 多项式的求值6.2 多项式计算【例6.116.11】已知多项式x4+8x3 10，分别取x=1.2和一个2 4矩阵为自变量计算该多项式的值。A=
19、1,8,0,0,-10;%4次多项式系数 x=1.2;%取自变量为一数值 y1=polyval(A,x)y1=5.8976 x=randi(9,2,4)x=8 2 6 3 9 9 1 5 y2=polyval(A,x)%分别计算矩阵x中各元素为自变量的多项式之值y2=8182 70 3014 287 12383 12383 -1 16156.2 多项式计算2矩阵多项式求值polyvalm函数用来求矩阵多项式的值，其调用格式与polyval相同，但含义不同。polyvalm函数要求x x为方阵。设A A为方阵，P代表多项式x3 5x2+8，polyvalm(P,A)的含义为A*A*A5*A*A+
20、8*eye(size(A)6.2.3 多项式的求值6.2 多项式计算【例6.126.12】以多项式x4+8x3 10为例，取一个2 2矩阵为自变量分别用polyval和polyvalm计算该多项式的值。A=1,8,0,0,-10;%多项式系数 x=-1,1.2;2,-1.8;%给出一个矩阵x y1=polyval(A,x)%计算代数多项式的值y1=-17.0000 5.8976 70.0000 -46.1584 y2=polyvalm(A,x)%计算矩阵多项式的值y2=-60.5840 50.6496 84.4160 -94.35046.2 多项式计算roots函数用于求多项式的全部根，其调用
21、格式为x=roots(P)P P为多项式的系数向量，求得的根赋给向量x，即x(1)，x(2)，x(n)分别代表多项式的n个根。6.2.4 多项式的求根6.2 多项式计算【例6.13】已知（1）计算f(x)=0的全部根。（2）由方程f(x)=0的根构造一个多项式g(x)，并与f(x)进行对比。P=2,-12,3,0,5;X=roots(P)%求方程f(x)=0的根X=5.7246+0.0000i0.8997+0.0000i -0.3122+0.6229i -0.3122-0.6229i G=poly(X)%求多项式g(x)G=1.0000 -6.0000 1.5000 -0.0000 2.500
22、0432()21235f xxxx6.2 多项式计算residue函数实现部份分式分解形式和多项式除式之间的相互转换。函数的用法如下。r,p,k=residue(b,a)b,a=residue(r,p,k)a、b 分别为分式的分母多项式、分子多项式的系数向量，r是分数多项式的商式的系数向量，p为分数多项式的极点，k为分数多项式的余式的系数向量6.2.5 多项式的除法变换6.2 多项式计算【例6.146.14】已知（1）将f(x)进行分式分解。（2）由分解的分式合成g(x)，并与f(x)进行对比。a=-4 0 8 3;r,p,k=residue(b,a)r=-1.4167 -0.6653 1.3
23、320p=1.5737 -1.1644 -0.4093k=-1.2500 b,a=residue(r,p,k)b=-1.2500 -0.7500 0.5000 -1.7500a=1.0000 -0.0000 -2.0000 -0.75006.3 数据插值实现一维插值的函数是interp1，其调用格式为vq=interp1(x,v,xq,method,extrapolation)x x、v v是两个等长的已知向量，分别存储采样点和采样值。若同一个采样点有多种采样值，则v v可以为矩阵，v v的每一列对应一种采样值。xqxq存储插值点，vqvq是一个列的长度与xqxq相同、宽度与v相同的矩阵。选项
24、methodmethod用于指定插值方法,可取值如下。linear（默认值）、pchip、spline、nearest、next、previous6.3.1 一维数据插值6.3 数据插值【例6.156.15】表6.4所示为我国06个月婴儿的体重、身长参考标准，用3次样条插值分别求得婴儿出生后半个月到5个半月每隔1个月的身长、体重参考值。tp=0:1:6;bb=50.6,3.27;56.5,4.97;59.6,5.95;62.3,6.73;.64.6,7.32;65.9,7.70;68.1,8.22;interbp=0.5:1:5.5;interbv=interp1(tp,bb,interbp,
25、spline)%用3次样条插值计算interbv=54.0847 4.2505 58.2153 5.5095 60.9541 6.3565 63.5682 7.0558 65.2981 7.5201 66.7269 7.91496.3 数据插值【例6.16】分别用pchip和spline方法对用以下函数生成的数据进行插值，绘制图形，对比两种插值方法的精度和运行效率。x1=1:7;subplot(1,2,1)y1=x1;y1(x15)=5;xq1=1:0.1:7;p1=interp1(x1,y1,xq1,pchip);s1=interp1(x1,y1,xq1,spline);plot(x1,y1
26、,ko,xq1,p1,r-,xq1,s1,b-.)subplot(1,2,2)x2=1:0.2:2*pi;y2=cos(5*x)./sqrt(x);xq2=1:0.1:2*pi;p2=interp1(x2,y2,xq2,pchip);s2=interp1(x2,y2,xq2,spline);plot(x2,y2,ko,xq2,p2,r-,xq2,s2,b-.)legend(Sample Points,pchip,spline)6.3 数据插值1.二维数据插值解决二维插值问题的函数interp2，其调用格式为Zq=interp2(X,Y,V,Xq,Yq,method,extrapval)X X、
27、Y Y分别存储采样点的平面坐标，V V存储采样点采样值。X Xq、Y Yq存储插值点的平面坐标，Z Zq是根据相应的插值方法得到的插值点的值。选项method的取值与一维插值函数相同，extrapval指定域外点的返回值。6.3.2 网格数据插值6.3 数据插值【例6.176.17】表6.5所示为某企业从19682008年、工龄为10年、20年和30年的职工的月均工资数据。试用线性插值求出19731993年每隔5年、工龄为15年和25年的职工月平均工资。x=1968:10:2008;h=10:10:30.;W=57,79,172,950,2496;69,95,239,1537,3703;87,
28、123,328,2267,4982;xi=1973:5:2003;hi=15;25;WI=interp2(x,h,W,xi,hi)WI=1.0e+03*0.0750 0.0870 0.1462 0.2055 0.7245 1.2435 2.1715 0.0935 0.1090 0.1963 0.2835 1.0928 1.9020 3.12236.3 数据插值2.2.多维数据插值多维数据插值MATLAB提供了3维、N维插值函数interp3、interpn，用法与interp2 一致。Vq=interp3(X,Y,Z,V,Xq,Yq,Zq,method)Vq=interpn(X1,X2,Xn,
29、V,Xq1,Xq2,Xqn,method)interp3函数的输入参数X、Y、Z以及interpn函数的输入参数 X1、X2、X3、.、Xn必须是网格格式。6.3.2 网格数据插值6.3 数据插值实现散乱数据插值的函数是griddata，调用格式为vq=griddata(x,y,v,xq,yq,method)vq=griddata(x,y,z,v,xq,yq,zq,method)x x、y y、z z存储采样点的坐标，v v是与采样点的采样值xqxq、yqyq、zqzq存储插值点的坐标，vq是根据相应的插值方法得到的插值结果。选项method指定插值方法，可取值如下。linear（默认值）、n
30、earest、natural、cubic、v46.3.3 散乱数据插值6.3 数据插值【例6.186.18】随机生成包含100个散点的数据集，绘制散点数据图和插值得到的网格数据图，观察插值结果。xy=rand(100,3)*10-5;x=xy(:,1);y=xy(:,2);z=xy(:,3);xq,yq=meshgrid(-4.9:0.08:4.9,-4.9:0.08:4.9);zq=griddata(x,y,z,xq,yq);mesh(xq,yq,zq)hold onplot3(x,y,z,rp)6.4 曲线拟合polyfit函数来求得最小二乘拟合多项式的系数。polyfit函数的调用格式为
31、p=polyfit(x,y,n)p,S=polyfit(x,y,n)p,S,mu=polyfit(X,Y,n)x x、y y是两个等长的向量，存储采样点x x和采样值y y产生一个n次多项式的系数向量p及其在采样点的误差向量S S。p是一个长度为n+1的向量，p的元素为多项式p1xn+p2xn1+.+pnx+pn+1的系数。mu是一个二元列向量，mu(1)是mean(x)，mu(2)是std(x)。6.4 曲线拟合【例6.196.19】某研究所为了研究氮肥的施肥量对土豆产量的影响，做了十次实验，实验数据如表6.6所示。试分析氮肥的施肥量与土豆产量之间的关系。data=0,15.18;34,21
32、.36;67,25.72;101,32.29;135,34.03;.202,39.45;259,43.15;336,43.46;404,40.83;471,30.75;x=data(:,1);y=data(:,2);f=polyfit(x,y,2);yi=polyval(f,x);plot(x,y,rp,x,yi)6.5 非线性方程和非线性方程组的数值求解fzero函数用来求非线性方程的根。该函数的调用格式为x=fzero(fun,x0,options)x,fval,exitflag,output=fzero(fun,x0,options)第2种格式可以在函数寻根失败时返回寻根过程的错误和信息
33、。fun是待求根的函数名，x0为搜索的起点。一个函数可能有多个根，但fzero函数只给出离x0最近的那个根。options为结构体变量，用于指定求解过程的优化参数，缺省时，使用默认值优化。返回方程的根，fval返回目标函数在解x处的值，exitflag返回求解过程终止原因，output返回寻根过程最优化的信息。6.5.1 非线性方程求解6.5 非线性方程和非线性方程组的数值求解fzero的优化参数通常调用optimset函数设置，optimset函数的调用方法如下。options=optimset(优化参数1,值1,优化参数2,值2,.)6.5.1 非线性方程求解6.5 非线性方程和非线性方程
34、组的数值求解【例6.206.20】求 e 2x x=0在x0=0附近的根。（1）建立函数文件funx.m。function fx=funx(x)fx=exp(-2*x)-x;（2）调用fzero函数求根。z=fzero(funx,0.0)z=0.4263如果要观测函数求根过程，可先设置优化参数，然后求解，命令如下。options=optimset(Display,iter);%设定显示迭代求解的中间结果 z=fzero(funx,0.0,options);6.5 非线性方程和非线性方程组的数值求解fsolve函数求非线性方程组F(X)=0的数值解。调用格式为：X=fsolve(fun,X0)X
35、,fval,exitflag,output=fsolve(fun,X0,options)第2种格式可以在函数寻根失败时返回寻根过程的错误和信息。fun是待求根的函数文件名，X0为搜索的起点。X返回方程组的解，fval返回目标函数在解X处的值，exitflag返回求解过程终止原因，output返回寻根过程最优化的信息。options为结构体变量，用于指定求解过程的优化参数，缺省时，使用默认值优化。6.5.2 非线性方程组的求解6.5 非线性方程和非线性方程组的数值求解fsolve的优化参数通常调用optimoptions函数设置，optimoptions函数的调用方法如下。options=opt
36、imoptions(SolverName,优化参数1,值1,优化参数2,值2,.)optimoptions函数用于设置除fzero、fminbnd、fminsearch和lsqnonneg这四个函数之外的其它优化函数的选项。SolverName是优化函数名，可以用字符串或函数句柄。6.5.2 非线性方程组的求解6.5 非线性方程和非线性方程组的数值求解【例6.216.21】求下列非线性方程组在（0.5，0.5）附近的数值解。（1）建立函数文件myfun.m。function q=myfun(p)x1=p(1);x2=p(2);q(1)=x12+x1-x22-1;q(2)=x2-sin(x12);（2）在给定的初值（0.5，0.5）下，调用fsolve函数求方程的根。x0=0.5;0.5;options=optimoptions(fsolve,Display,off);%不显示中间结果x=fsolve(myfun,x0,options)x=0.7260 0.5029221122211sin0 xxxxx
