第5章 CHAPTER 5 MATLAB数值计算 本章要点：  多项式运算；  数据插值；  多项式拟合；  数据统计；  数值计算。 5.1多项式 多项式在代数中占有重要的地位，广泛用于数据插值、数据拟合和信号与系统等应用领域。MATLAB提供了各种多项式的创建以及运算方法，应用起来简单方便。 5.1.1多项式的创建 一个多项式按降幂排列为 p(x)=anxn+an-1xn-1+…+a1x+a0(51) 在MATLAB中多项式的各项系数用一个行向量表示，使用长度为n+1的行向量按降幂排列，多项式中某次幂的缺项用0表示，则表示为 p=［an,an-1,…,a1,a0］(52) 例如，多项式p1(x)=x3-2x2+4x+6，在MATLAB可以表示为p1=［1,-2,4,6］； p2(x)=x3+3x+6可表示为p2=［1,0,3,6］。 在MATLAB中，创建一个多项式，可以用poly2str、poly2sym函数实现，其调用格式如下： f =poly2str(p,'x')%p为多项式的系数，x为多项式的变量 f =poly2sym(p)%p为多项式的系数 其中，f=poly2str(p,'x')表示创建一个系数为p，变量为x的字符串型多项式； f=poly2sym(p)表示创建一个系数为p，默认变量为x的符号型多项式。两者在命令窗口显示形式类似，但数据类型是不一样的，一个是字符串型，另一个是符号型。 微课视频 【例51】已知多项式系数为p=［1,-2,4,6］， 分别用poly2str(p,'x')和poly2sym(p)创建多项式，比较它们有什么不同。 程序代码如下： >> p=［1 -2 4 6］ p = 1-246 >> f1=poly2str(p,'x') f1 = 'x^3 - 2 x^2 + 4 x + 6' >> f2=poly2sym(p) f2 = x^3 - 2*x^2 + 4*x + 6 显然，两种函数创建的多项式f1和f2显示形式类似，但数据类型和大小都不一样，如图51所示。 图51两种多项式的比较 5.1.2多项式的值和根 1. 多项式的值 在MATLAB里，求多项式的值可以用polyval和polyvalm函数。它们的输入参数都是多项式系数和自变量，二者的区别是前者为代数多项式求值，后者为矩阵多项式求值。 (1) 代数多项式求值 polyval函数可以求解代数多项式的值，其调用格式如下： y=polyval(p,x) 其中，p为多项式的系数，x为自变量。当x为一个数值时，求解的是多项式在该点的值； 若x为向量或矩阵，则是对向量或矩阵每个元素求多项式的值。 微课视频 【例52】已知多项式f(x)=x3-2x2+4x+6，分别求 x1=2和x=［0,2,4,6,8,10］向量的多项式的值。 在文件编辑窗口编写命令文件，保存为exam_5_2.m脚本文件。程序代码如下： x1=2; x=［0:2:10］; p=［1 -2 4 6］; y1=polyval(p,x1) y=polyval(p,x) 在命令空间输入文件名exam_5_2.m，就能直接运行该脚本文件。结果如下： >> exam_5_2 y1 = 14 y = 61454174422846 (2) 矩阵多项式求值 polyvalm函数是以矩阵为自变量求解多项式的值，其调用格式如下： Y=polyvalm(p,X) 其中，p为多项式系数，X为自变量，自变量要求为方阵。 因为运算规则不一样，所以MATLAB用polyvalm和polyval函数求解多项式的值是不一样的。例如，假设A为方阵，p为多项式x2-5x+6的系数，则polyvalm(p,A)表示A*A-5*A+6*eye(size(A))，而polyval(p,A)表示A.*A-5*A+6*ones(size(A))。 微课视频 【例53】已知多项式为f(x)=x2-3x+2，分别用polyvalm和polyval函数，求X=12 34的多项式值。 在文件编辑窗口编写命令文件，保存为exam_5_3.m脚本文件。程序代码如下： X=［1 2;3 4］; p=［1 -3 2］; Y=polyvalm(p,X) Y1=polyval(p,X) 程序运行结果： >> exam_5_3 Y = 64 612 Y1 = 00 26 2. 多项式的根 一个n次多项式有n个根，这些根有可能是实根，也有可能包含若干对共轭复根。MATLAB提供了roots函数用于求解多项式的全部根，其调用格式为 r=roots(p) 其中，p为多项式的系数向量，r为多项式的根向量，r(1),r(2),…,r(n)分别表示多项式的n个根。 MATLAB还提供一个由多项式的根求多项式的系数的函数poly，其调用格式如下： p=poly(r) 其中，r为多项式的根向量，p为由根r构造的多项式系数向量。 微课视频 【例54】已知多项式为f(x)=x4+4x3-3x+2。 (1) 用roots函数求该多项式的根r。 (2) 用poly函数求根为r的多项式系数。 在文件编辑窗口编写命令文件，保存为exam_5_4.m脚本文件。程序代码如下： p=［1 4 0 -3 2］; r=roots(p) p1=poly(r) 程序运行结果： >> exam_5_4 r = -3.7485 + 0.0000i -1.2962 + 0.0000i 0.5224 + 0.3725i 0.5224 - 0.3725i p1 = 1.00004.0000-0.0000-3.00002.0000 显然，roots和poly函数的功能正好相反。 5.1.3多项式的四则运算 多项式之间可以进行四则运算，其结果仍为多项式。在MATLAB中，用多项式系数向量进行四则运算，得到的结果也显示为多项式系数向量。 1. 多项式的加减运算 MATLAB没有提供多项式加减运算的函数。事实上多项式的加减运算，是合并同类项，可以用多项式系数向量进行加减运算。如果多项式阶次不同，则把低次多项式系数缺少的高次项用0补足，使得多项式系数矩阵具有相同维度，以便实现加减运算。 2. 多项式乘法运算 在MATLAB中，两个多项式的乘积可以用conv函数实现。其调用格式为 p=conv(p1,p2) 其中，p1和p2是两个多项式的系数向量； p是两个多项式乘积的系数向量。 3. 多项式除法运算 MATLAB用deconv函数实现两个多项式的除法运算。其调用格式为 ［q,r］=deconv(p1,p2) 其中，q为多项式p1除以p2的商式； r为多项式p1除以p2的余式。q和r都是多项式系数向量。 deconv是conv的逆函数，因此满足p1=conv(p2，q)+r。 微课视频 【例55】已知两个多项式为 f(x)=x4+4x3-3x+2，g(x)=x3-2x2+x。 (1) 求两个多项式相加f(x)+g(x),两个多项式相减f(x)-g(x)的结果。 (2) 求两个多项式相乘f(x)*g(x),两个多项式相除f(x)/g(x)的结果。 在文件编辑窗口编写命令文件，保存为exam_5_5.m脚本文件。程序代码如下： p1=［1 4 0 -3 2］; p2=［0 1 -2 1 0］; p3=［1 -2 1 0］; p=p1+p2%f(x)+g(x) poly2sym(p) p=p1-p2%f(x)-g(x) poly2sym(p) p=conv(p1,p2)%f(x)*g(x) poly2sym(p) ［q,r］=deconv(p1,p3)%f(x)/g(x) p4=conv(q,p3)+r%验证deconv是conv的逆函数 程序运行结果： >> exam_5_5 p = 15-2-22 ans = x^4 + 5*x^3 - 2*x^2 - 2*x + 2 p = 132-42 ans = x^4 + 3*x^3 + 2*x^2 - 4*x + 2 p = 012-718-720 ans = x^7 + 2*x^6 - 7*x^5 + x^4 + 8*x^3 - 7*x^2 + 2*x q = 16 r = 0011-92 p4 = 140-32 5.1.4多项式的微积分运算 1. 多项式的微分 对n阶多项式p(x)=anxn+an-1xn-1+…+a1x+a0求导，其导数为n-1阶多项式dp(x)=nanxn-1+(n-1)an-1xn-2+…+a1。原多项式及其导数多项式的系数分别为p=［an,an-1,…,a1,a0］, d=［nan,(n-1)an-1,…,a1］。 在MATLAB中，可以用polyder函数来求多项式的微分运算，polyder函数可以对单个多项式求导，也可以对两个多项式的乘积或商求导，其调用格式如下： p=polyder(p1)%求多项式p1的导数 p=polyder(p1,p2)%求多项式p1×p2的积的导数 ［p,q］=polyder(p1,p2)%p1÷p2的导数，p为导数的分子多项式系数，q为导数的分母多 %项式系数 微课视频 【例56】 已知两个多项式为f(x)=x4+4x3-3x+2，g(x)=x3-2x2+x。 (1) 求多项式f(x)的导数。 (2) 求两个多项式乘积f(x)*g(x)的导数。 (3) 求两个多项式相除g(x)/f(x)的导数。 在文件编辑窗口编写命令文件，保存为exam_5_6.m脚本文件。程序代码如下： p1=［1 4 0 -3 2］; p2=［1 -2 1 0］; p=polyder(p1) poly2sym(p) p=polyder(p1,p2) poly2sym(p) ［p,q］=polyder(p2,p1) 程序运行结果： >> exam_5_6 p = 4120-3 ans = 4*x^3 + 12*x^2 - 3 p = 712-35424-142 ans = 7*x^6 + 12*x^5 - 35*x^4 + 4*x^3 + 24*x^2 - 14*x + 2 p = -145-1412-82 q = 1816-6-20169-124 2. 多项式的积分 对于n阶多项式p(x)=anxn+an-1xn-1+…+a1x+a0，其不定积分为n+1阶多项式i(x)=1n+1anxn+1+1nan-1xn+…+12a1x2+a0x+k，其中k为常数项。原多项式和积分多项式分别可以表示为系数向量p=［an,an-1,…,a1,a0］, I=1n+1an,1nan-1,…,12a1,k。 在MATLAB中，提供了polyint函数用于多项式的积分。其调用格式为 I=polyint(p,k)%求以p为系数的多项式的积分，k为积分常数项 I=polyint(p)%求以p为系数的多项式的积分，积分常数项为默认值0 显然polyint是polyer的逆函数，因此有p=polyder(I)。 微课视频 【例57】求多项式的积分I=∫(x4+4x3-3x+2)dx。 在文件编辑窗口编写命令文件，保存为exam_5_7.m脚本文件。程序代码如下： p=［1 4 0 -3 2］; I=polyint(p)%求多项式的积分，常数项为默认的0 poly2sym(I)%显示多项式积分的多项式 p=polyder(I)%验证polyint是polyder的逆函数 syms k%定义常数项k I1=polyint(p,k)%求多项式的积分，常数项为k poly2sym(I1) 程序运行结果： >> exam_5_7 I = 0.20001.00000-1.50002.00000 ans = x^5/5 + x^4 - (3*x^2)/2 + 2*x p = 140-32 I1 = ［ 1/5, 1, 0, -3/2, 2, k］ ans = x^5/5 + x^4 - (3*x^2)/2 + 2*x + k 5.1.5多项式的部分分式展开 由分子多项式B(s)和分母多项式A(s)构成的分式表达式进行多项式的部分分式展开，表达式如下： B(s)A(s)=r1s-p1+r2s-p2+…+rns-pn+k(s)(53) MATLAB可以用residue函数实现多项式的部分分式展开，residue函数的调用格式如下： ［r,p,k］=residue(B,A) 其中，B为分子多项式系数行向量； A为分母多项式系数行向量； ［p1;p2;…;pn］为极点列向量； ［r1;r2;…;rn］为零点列向量； k为余式多项式行向量。 residue函数还可以将部分分式展开式转换为两个多项式的分式表达式，其调用格式为： ［B,A］=residue(r,p,k) 微课视频 【例58】已知分式表达式为f(s)=B(s)A(s)=3s3+1s2-5s+6。 (1) 求f(s)的部分分式展开式。 (2) 将部分分式展开式转换为分式表达式。 在文件编辑窗口编写命令文件，保存为exam_5_8.m脚本文件。程序代码如下： a=［1 -5 6］; b=［3 0 0 1］; ［r,p,k］=residue(b,a)%部分分式展开 ［b1,a1］=residue(r,p,k)%将部分分式展开转换为分式表达式 程序运行结果： >> exam_5_8 r = 82.0000 -25.0000 p = 3.0000 2.0000 k = 315 b1 = 3001 a1 = 1-56 5.2数据插值 在工程测量与科学实验中，得到的数据通常都是离散的。如果要得到这些离散数据点以外的其他数据值，就需要根据这些已知数据进行插值。假设已测量得到n个点数据，(x1,y1)、(x2y2)、…、(xnyn)，且这些测量值满足某一个未知的函数关系y=f(x)。数据插值的任务就是根据这n个测量数据，构造一个函数y=p(x)，使得yi=p(xi)(i=1,2,…,n)成立，称p(x)为f(x)关于点x1,x2,…,xn的插值函数。求插值函数p(x)的方法为插值法。插值函数p(x)一般可以用线性函数、多项式或样条函数实现。 根据插值函数的自变量个数，数据插值可以分为一维插值、二维插值和多维插值等； 根据不同的插值函数，又可以分为线性插值、多项式插值和样条函数插值等。MATLAB提供了一维插值interp1、二维插值interp2、三维插值interp3和N维插值interpn函数，以及三次样条插值spline函数等。 5.2.1一维插值 所谓一维插值是指被插值函数的自变量是一个单变量的函数。一维插值采用的方法一般有一维多项式插值、一维快速插值和三次样条插值。 1. 一维多项式插值 MATLAB中提供了interp1函数进行一维多项式插值。interp1函数使用了多项式函数，通过已知数据点，计算目标插值点的数据。interp1函数调用格式如下： yi=interp1(Y,xi) 其中，Y是在默认自变量x选为1： n的值。 yi=interp1(X,Y,xi) 其中X和Y是长度一样的已知向量数据，xi可以是一个标量，也可以是向量。 yi=interp1(X,Y,xi,'method') 其中，method是插值方法，其取值有下面几种： (1) linear线性插值： 这是默认插值方法，它将与插值点靠近的两个数据点直线连接，在直线上选取对应插值点的数据。这种插值方法兼顾速度和误差，插值函数具有连续性，但平滑性不好。 (2) nearest最邻近点插值： 根据插值点和最接近的已知数据点进行插值，这种插值方法速度快，占用内存小，但一般误差最大，插值结果最不平滑。 (3) next下一点插值： 根据插值点和下一点的已知数据点插值，这种插值方法的优缺点与最邻近点插值一样。 (4) previous前一点插值： 根据插值点和前一点的已知数据点插值，这种插值方法的优缺点与最邻近点插值一样。 (5) spline三次样条插值： 采用三次样条函数获得插值点数据，要求在各点处具有光滑条件。这种插值方法连续性好，插值结果最光滑，缺点为运行时间长。 (6) cubic 三次多项式插值： 根据已知数据求出一个3次多项式进行插值。这种插值方法连续性好，光滑性较好，缺点是占用内存多，速度较慢。 需要注意，xi的取值如果超出已知数据X的范围，就会返回NaN错误信息。 MATLAB还提供interp1q函数用于一维插值。它与interp1函数的主要区别是，当已知数据不是等间距分布时，interp1q插值速度比interp1快。需要注意，interp1q执行的插值数据x必须是单调递增的。 微课视频 【例59】某气象台对当地气温进行测量，实测数据见表51所示，用不同插值方法计算t=12时的气温。 表51某地不同时间的气温 测量时间t(小时)681014161820 温度T(度)1617.519.32221.219.518 在文件编辑窗口编写命令文件，保存为exam_5_9.m脚本文件。程序代码如下： t=［6 81014161820］;%测量时间t T=［16 17.5 19.3 2221.2 19.5 18］;%测量的温度T t1=12;%插值点时间t1 T1=interp1(t,T,t1,'nearest')%最接近点插值 T2=interp1(t,T,t1,'linear')%线性插值 T3=interp1(t,T,t1,'next')%下一点插值 T4=interp1(t,T,t1,'previous')%前一点插值 T5=interp1(t,T,t1,'pchip')%三次多项式插值 T6=interp1(t,T,t1,'spline')%三次样条插值 程序运行结果： >> exam_5_9 T1 = 22 T2 = 20.6500 T3 = 22 T4 = 19.3000 T5 = 21.0419 T6 = 21.1193 微课视频 【例510】假设测量的数据来自函数f(x)=e-0.5xsinx，试根据生成的数据，用不同方法进行插值，比较插值结果。 在文件编辑窗口编写命令文件，保存为exam_5_10.m脚本文件。程序代码如下： clear x=(0:0.4:2*pi)'; y=exp(-0.5*x).*sin(x);%生成测试数据 x1=(0:0.1:2*pi)';%插值点 y0=exp(-0.5*x1).*sin(x1);%插值点真实值 y1=interp1(x,y,x1,'nearest');%最接近点插值 disp('interp1函数插值时间');tic y2=interp1(x,y,x1); toc;%interp1插值时间 y3=interp1(x,y,x1,'spline') ;%三次样条插值 disp('interp1q函数插值时间');tic yq=interp1q(x,y,x1);toc;%interp1q插值时间 plot(x1,y1,'--',x1,y2,'-',x1,y3,'-.',x,y,'*',x1,y0,':') legend('nearest插值数据','linear插值数据','spline插值数据',... '样本数据点','插值点真实数据') max(abs(y0-y3)) 程序运行结果如下，插值效果如图52所示。 >> exam_5_10 interp1函数插值时间 历时 0.001270 秒。 interp1q函数插值时间 历时 0.001472 秒。 ans = 6.5673e-04 图52各种插值结果比较 由上面结果可知，最接近点拟合误差大，直线拟合得到的曲线不平滑； 三次样条插值的效果最好，曲线平滑，误差很小，基本逼近真实值。 2. 一维快速傅里叶插值 在MATLAB中，一维快速傅里叶插值可以用interpft函数实现。该函数利用傅里叶变换将输入数据变换到频率域，然后用更多点进行傅里叶逆变换，实现对数据的插值。函数调用格式为 y=interpft(x,n)%表示对x进行傅里叶变换，然后采用n点傅里叶逆变换，得到插值后的数据 y=interpft(x,n,dim)%表示在dim维上进行傅里叶插值 微课视频 【例511】假设测量的数据来自函数f(x)=sinx，试根据生成的数据，用一维快速傅里叶插值，并比较插值结果。 在文件编辑窗口编写命令文件，保存为exam_5_11.m脚本文件。程序代码如下： clear x=0:0.4:2*pi; y=sin(x);%原始数据 N=length(y); M=N*4; x1=0:0.1:2*pi; y1=interpft(y,M-1);%傅里叶插值 y2=sin(x1);%插值点真实数据 plot(x,y,'O',x1,y1,'*',x1,y2,'-') legend('原始数据','傅里叶插值数据','插值点真实数据') max(abs(y1-y2)) 程序运行结果如下，插值效果如图53所示。 >> exam_5_11 ans = 0.0980 图53一维快速傅里叶插值及比较 由以上结果可知，一维快速傅里叶插值实现插值的速度比较快，曲线平滑，误差很小，基本逼近真实值。 3. 三次样条插值 三次样条插值是利用多段多项式逼近插值，降低了插值多项式的阶数，使得曲线更为光滑。在MATLAB中，interp1插值函数的method选为spline样条插值选项，就可以实现三次样条插值。另外，MATLAB专门提供三次样条插值函数spline，其格式如下： yi=spline(x,y,xi)%利用初始值x，y，对插值点数据xi进行三次样条插值。采用这种调 %用方式，相当于yi=interp1(x,y,xi,'spline') 微课视频 【例512】已知数据x=［-5 -4 -3 -2 -1 0 1 2 3 4 5］,y=［25 16 9 4 1 0 1 4 9 16 25］,对xi=-5:0.5:5，用spline进行三次样条插值，并比较用interp1实现三次样条插值。 在文件编辑窗口编写命令文件，保存为exam_5_12.m脚本文件。程序代码如下： x=-5:5; y=x.*x; xi=-5:0.5:5; y0=xi.*xi; y1=spline(x,y,xi); y2=interp1(x,y,xi,'spline'); plot(x,y,'O',xi,y0,xi,y1,'+',xi,y2,'*') legend('原始数据','插值点真实数据','spline插值','interp1样条插值') max(abs(y1-y2)) 程序运行结果如下，插值效果如图54所示。 >> exam_5_12 ans = 0 由程序结果可知，三次样条插值spline函数实现插值的效果和interp1(x,y,xi,'spline')一样。 图54三次样条插值及比较 5.2.2二维插值 二维插值是指已知一个二元函数的若干个采样数据点x，y和z(x,y)，求插值点(x1，y1)处的z1的值。在MATLAB中，提供interp2函数用于实现二维插值，其调用格式为 Z1=interp2(X,Y,Z,X1,Y1,'method') 其中，X和Y是两个参数的采样点，一般是向量，Z是参数采样点对应的函数值。X1和Y1是插值点，可以是标量也可以是向量。Z1是根据选定的插值方法(method)得到的插值结果。插值方法method和一维插值函数相同，linear为线性插值(默认算法)，nearest为最近点插值，spline为三次样条插值，cubic为三次多项式插值。需要注意，X1和Y1不能超出X和Y的取值范围，否则会得到NaN错误信息。 微课视频 【例513】某实验对电脑主板的温度分布做测试。用x表示主板的宽度(cm)，y表示主板的深度(cm)，用T表示测得的各点温度(℃)，测量结果如表52所示。 (1) 分别用最近点二维插值和线性二维插值法求(12.6,7.2)点的温度。 (2) 用三次多项式插值求主板宽度每1cm，深度每1cm 处各点的温度，并用图形显示插值前后主板的温度分布图。 表52主板各点温度测量值 y x 0510152025 0303234333231 5333741383533 10353844433734 15323436353332 在文件编辑窗口编写命令文件，保存为exam_5_13.m脚本文件。程序代码如下： clear x=［0:5:25］; y=［0:5:15］'; T=［303234333231; 333741383533; 353844433734; 323436353332］; x1=12.6;y1=7.2;%插值点(12.6,7,2) T1=interp2(x,y,T,x1,y1,'nearest')%最近点二维插值 T2=interp2(x,y,T,x1,y1,'linear')%线性二维插值 xi=［0:1:25］; yi=［0:1:15］'; Ti=interp2(x,y,T,xi,yi,'cubic');%三次多项式二维插值 subplot(1,2,1) mesh(x,y,T) xlabel('板宽度(cm)');ylabel('板深度(cm)');zlabel('温度(摄氏度)') title('插值前主板温度分布图') subplot(1,2,2) mesh(xi,yi,Ti) xlabel('板宽度(cm)');ylabel('板深度(cm)');zlabel('温度(摄氏度)') title('插值后主板温度分布图') 运行程序，结果如下，图55是插值前后主板温度分布图。可知，用插值技术处理数据，可以使得温度分布图更加光滑。 >> exam_5_13 T1 = 38 T2 = 41.2176 图55插值前后主板温度分布图 5.2.3多维插值 1. 三维插值 在MATLAB中，还提供了三维插值的函数interp3，其调用格式为： U1=interp3(X,Y,Z,U，X1,Y1,Z1,'method') 其中，X、Y、 Z是三个参数的采样点，一般是向量，U是参数采样点对应的函数值。X1、Y1、Z1是插值点，可以是标量也可以是向量。U1是根据选定的插值方法(method)得到的插值结果。插值方法method和一维插值函数相同，linear为线性插值(默认算法)，nearest为最近点插值，spline为三次样条插值，cubic为三次多项式插值。需要注意，X1、Y1和Z1不能超出X、Y和Z的取值范围，否则会得到NaN错误信息。 2. n维插值 在MATLAB中，还可以实现更高维的插值，interpn函数用于实现n维插值。其调用格式为： U1=interpn(X1,X2,…,Xn,U,Y1,Y2,…,Yn,'method') 其中，X1,X2,…，Xn是n个参数的采用点，一般是向量，U是参数采样点对应的函数值。Y1,Y2,…，Yn是插值点，可以是标量也可以是向量。U1是根据选定的插值方法(method)得到的插值结果。插值方法method和一维插值函数相同，linear为线性插值(默认算法)，nearest为最近点插值，spline为三次样条插值，cubic为三次多项式插值。需要注意，Y1,Y2,…，Yn不能超出X1,X2,…，Xn的取值范围，否则会得到NaN错误信息。 5.3数据拟合 与数据插值类似，数据拟合的目的也是用一个较为简单的函数g(x)去逼近一个未知的函数f(x)。利用已知的测量数据(xi,yi)(i=1,2,…,n)，构造函数y=g(x)，使得误差δi=g(xi)-f(xi)(i=1,2,…,n)在某种意义上达到最小。 一般用得比较多的是多项式拟合，所谓多项式拟合是利用已知的测量数据(xi,yi)(i=1,2,…,n)，构造一个m(m> exam_5_14 p1 = -0.00240.0462-0.27820.47600.1505 s1 = 包含以下字段的 struct: R: ［5×5 double］ df: 25 normr: 0.4086 g1 = '-0.002378 x^4 + 0.04625 x^3 - 0.27815 x^2 + 0.476 x + 0.15048' p2 = 0.0007-0.01910.1856-0.75931.08260.0046 s2 = 包含以下字段的 struct: R: ［6×6 double］ df: 24 normr: 0.0909 g2 = '0.00071166 x^5 - 0.019146 x^4 + 0.18564 x^3 - 0.75929 x^2 + 1.0826 x + 0.0045771' 图564阶多项式和5阶多项式拟合f(x)函数 由上述例题结果可知，用高阶多项式拟合f(x)函数的效果更好，误差小，更加逼近实际函数f(x)。 5.4数据统计 在科学研究和生产实际中经常需要对数据进行统计， MATLAB语言提供了很多数据统计方面的函数。 5.4.1矩阵元素的最大值和最小值 1. 求向量的最大元素和最小元素 1) 求向量的最大元素 在MATLAB中，可以用函数max(X)求一个向量X的最大元素，其调用格式为 y=max(X)%返回向量X的最大元素给y，如果X中包括复数元素，则按模取最大元素 ［y,k］=max(X)%返回向量X的最大元素给y，最大元素所在的位置序号给k，如果X中包括复数 %元素，则按模取最大元素 例如，求向量X=［34,23,-23,6,76,56,14,35］的最大值。 >> X=［34,23,-23,6,76,56,14,35］; >> y=max(X) y = 76 >> ［y,k］=max(X) y = 76 k = 5 2) 求向量的最小元素 在MATLAB中，可以用函数min(X)求一个向量X的最小元素，其调用格式及用法与max(X)函数一样。 例如，求向量X=［34,10,-23,6,76,0,14,35］的最小值。 >> X=［34,10,-23,6,76,0,14,35］; >> y=min(X) y = -23 >> ［y,k］=min(X) y = -23 k = 3 2. 求矩阵的最大元素和最小元素 1) 求矩阵的最大元素 在MATLAB中，可以用函数max求一个矩阵A的最大元素，其调用格式为： Y=max(A)%返回矩阵A的每列上最大元素给Y，Y是一个行向量 ［Y, K］=max(A) %返回矩阵A的每列上最大元素给Y，K向量记录每列最大元素所在的行号 %如果X中包括复数元素，则按模取最大元素 ［Y, K］=max(A, ［］, dim) 其中dim为1时，该函数和max(A)完全相同。当dim为2时，该函数返回一个每行上最大元素的列向量。 2) 求矩阵的最小元素 在MATLAB中，可以用函数min求一个矩阵A的最小元素，其调用格式及用法和max函数一样。 微课视频 【例515】在MATLAB中，用max和min函数，求矩阵A的每行和每列的最大和最小元素，并求整个A的最大和最小元素。 A=121-624 -423120 2-3186 45316-7 程序代码如下： >> A=［12 1 -6 24;-4 23 12 0;2 -3 18 6;45 3 16 -7］; >> Y1=max(A,［］,2)%求每行最大元素 Y1 = 24 23 18 45 >> ［Y2,K］=min(A,［］,2)%求每行最小元素，及每行最小值的列数 Y2 = -6 -4 -3 -7 K = 3 1 2 4 >> Y3=max(A)%求每列的最大元素 Y3 = 45231824 >> ［Y4,K1］=min(A)%求每列的最小元素，及最小元素所在的行数 Y4 = -4-3-6-7 K1 = 2314 >> ymax=max(max(A))%求矩阵A的最大元素 ymax = 45 >> ymin=min(min(A))%求矩阵A的最小元素 ymin = -7 3. 两个维度一样的向量或矩阵对应元素比较 max和min函数还能对两个维度一样的向量或矩阵对应元素求较大值和较小值，其调用格式为： Y=max(A, B) 其中，A和B是同维度的向量或矩阵，Y的每个元素为A和B对应元素的较大者，与A和B同维。 min函数的用法和max一样。 例如，求A和B矩阵对应元素的较大元素Y1和较小元素Y2。 程序代码如下： >> A=［1 5 6;7 3 1;3 7 4］ A = 156 731 374 >> B=［2 9 4;9 1 3;-1 0 3］ B = 294 913 -103 >> Y1=max(A,B) Y1 = 296 933 374 >> Y2=min(A,B) Y2 = 154 711 -103 5.4.2矩阵元素的平均值和中值 数据序列的平均值指的是算术平均，中值是指数据序列中其值位于中间的元素，如果数据序列个数为偶数，中值等于中间两项的平均值。 在MATLAB中，求矩阵或向量元素的平均值用mean函数，求中值用median函数。它们的调用方法如下： (1) y=mean(X)%返回向量X的算术平均值 (2) Y=mean(A)%返回一个矩阵A每列的算术平均值的行向量 (3) y=median(X)%返回向量X的中值 (4) Y=median(A)%返回一个矩阵A每列的中值的行向量 (5) Y=mean(A，dim) %当dim为1时，等同于mean(A)； 当dim为2时，返回一个矩阵A %每行的算术平均值的列向量 (6) Y=median(A，dim)%当dim为1时，等同于median(A)； 当dim为2时，返回一个 %矩阵A每行的中值的列向量 例如，求向量X和矩阵A的平均值和中值。 程序代码如下： >> X=［1,12,23,7,9,-5,30］; >> y1=mean(X) y1 = 11 >> y2=median(X) y2 = 9 >> A=［0 9 2;7 3 3;-1 0 3］ A = 092 733 -103 >> Y1=mean(A) Y1 = 2.00004.00002.6667 >> Y2=median(A) Y2 = 033 >> Y3=mean(A,2) Y3 = 3.6667 4.3333 0.6667 >> Y4=median(A,2) Y4 = 2 3 0 5.4.3矩阵元素的排序 在MATLAB中，可以用函数sort实现数据序列的排序。对于向量X的排序，可以用函数sort(X)，函数返回一个对向量X的元素按升序排列的向量。 sort函数还可以对矩阵A的各行或各列的元素重新排序，其调用格式为 ［Y,I］=sort(A, dim, mode) 其中，当dim为1时，矩阵元素按列排序； 当dim为2时，矩阵元素按行排序。dim默认为1。当mode为'ascend'，则按升序排序； 当mode为'descend'，则按降序排序。mode默认取'ascend'。Y为排序后的矩阵，而I记录Y中元素在A中的位置。 例如，对一个向量X和一个矩阵A做各种排序。 程序代码如下： >> X=［1,12,23,7,9,-5,30］; >> Y=sort(X) Y = -5179122330 >> A=［0 9 2;7 3 1;-1 0 3］ A = 092 731 -103 >> Y1=sort(A) Y1 = -101 032 793 >> Y2=sort(A,1,'descend') Y2 = 793 032 -101 >> Y3=sort(A,2,'ascend') Y3 = 029 137 -103 >> ［Y4,I］=sort(A,2,'descend') Y4 = 920 731 30-1 I = 231 123 321 5.4.4矩阵元素求和与求积 在MATLAB中，向量和矩阵求和与求积的基本函数是sum和prod，它们的使用方法类似，调用格式为： (1) y=sum(X)%返回向量X各元素的和 (2) y=prod(X)%返回向量X各元素的乘积 (3) Y=sum(A)%返回一个矩阵A各列元素的和的行向量 (4) Y=prod(A)%返回一个矩阵A各列元素的乘积的行向量 (5) Y=sum(A，dim) %当dim为1时，该函数等同于sum(A)； 当dim为2时，返回一个 %矩阵A各行元素的和的列向量 (6) Y=prod(A，dim) %当dim为1时，该函数等同于prod(A)； 当dim为2时，返回一个 %矩阵A各行元素的乘积的列向量 例如，求一个向量X和一个矩阵A的各元素的和与乘积。 程序代码如下： >> X=［1,3,9,-2,7］; >> y=sum(X)%求向量X的各元素的和 y = 18 >> y=prod(X)%求向量X的各元素的乘积 y = -378 >> A=［1 9 2;7 3 1;-1 1 3］ A = 192 731 -113 >> Y1=sum(A)%求矩阵A的各列元素的和 Y1 = 7136 >> Y2=sum(A,2)%求矩阵A的各行元素的和 Y2 = 12 11 3 >> Y3=prod(A)%求矩阵A的各列元素的乘积 Y3 = -7276 >> Y4=prod(A,2)%求矩阵A的各行元素的乘积 Y4 = 18 21 -3 >> y5=sum(Y1)%求矩阵A所有元素的和 y5 = 26 >> y6=prod(Y3)%求矩阵A所有元素的乘积 y6 = -1134 5.4.5矩阵元素的累加和与累乘积 在MATLAB中，向量和矩阵的累加和与累乘积的基本函数是cumsum和cumprod，它们的使用方法类似，调用格式为 (1) y=cumsum(X)%返回向量X累加和向量 (2) y=cumprod(X)%返回向量X累乘积向量 (3) Y=cumsum(A)%返回一个矩阵A各列元素的累加和的矩阵 (4) Y=cumprod(A)%返回一个矩阵A各列元素的累乘积的矩阵 (5) Y=cumsum(A，dim) %当dim为1时，该函数等同于cumsum(A)； 当dim为2时， %返回一个矩阵A各行元素的累加和矩阵 (6) Y=cumprod(A，dim) %当dim为1时，该函数等同于cumprod(A)； 当dim为2时， %返回一个矩阵A各行元素的累乘积矩阵 例如，求一个向量X和一个矩阵A的各元素的累加和与累乘积。 程序代码如下： >> X=［1,3,9,-2,7］; >> Y=cumsum(X) Y = 14131118 >> Y=cumprod(X) Y = 1327-54-378 >> A=［1 9 2;7 3 1;-1 1 3］ A = 192 731 -113 >> Y1=cumsum(A) Y1 = 192 8123 7136 >> Y2=cumsum(A,2) Y2 = 11012 71011 -103 >> Y3=cumprod(A) Y3 = 192 7272 -7276 >> Y4=cumprod(A,2) Y4 = 1918 72121 -1-1-3 5.4.6标准方差和相关系数 1. 标准方差 对于具有N个元素的向量数据x1,x2,…,xN，有如下两种标准方差的公式： D1=1N-1∑Ni=1 (xi-x-)2(55) 或 D2=1N∑Ni=1(xi-x-)2(56) 其中， x-=1N∑Ni=1xi(57) 在MATLAB中，可以用函数std计算向量和矩阵的标准方差。对于向量X，std(X)返回一个标准方差； 对于矩阵A，std(A)返回一个矩阵A各列或者各行的标准方差向量。std函数的调用格式为： (1) d=std(X)%求向量X的标准方差 (2) D=std(A，flag，dim) 其中，当dim为1时，求矩阵A的各列元素的标准方差； 当dim为2时，则求矩阵A的各行元素的标准方差。当flag为0时，按公式D1计算标准方差； 当flag为1时，按D2计算标准方差。默认flag=0，dim=1。 例如，求一个向量X和一个矩阵A的标准方差。 程序代码如下： >> X=［1,3,9,-2,7］; >> d=std(X)%求向量X的标准方差 d = 4.4497 >> A=［1 9 2;7 3 1;-1 1 3］ A = 192 731 -113 >> D1=std(A,0,1)%按D1标准方差公式，求矩阵A的列元素标准方差 D1 = 4.16334.16331.0000 >> D2=std(A,0,2)%按D1标准方差公式，求矩阵A的行元素标准方差 D2 = 4.3589 3.0551 2.0000 >> D3=std(A,1,1)%按D2标准方差公式，求矩阵A的列元素标准方差 D3 = 3.39933.39930.8165 >> D4=std(A,1,2)%按D2标准方差公式，求矩阵A的行元素标准方差 D4 = 3.5590 2.4944 1.6330 2. 相关系数 对于两组数据序列xi,yi(i=1,2,…,N)，可以用下列式子定义两组数据的相关系数： ρ=∑Ni=1(xi-x-)(yi-y-)∑Ni=1(xi-x-)2∑Ni=1(yi-y-)2(58) 其中， x-=1N∑Ni=1xi，y-=1N∑Ni=1yi(59) 在MATLAB中，可以用corrcoef函数计算数据的相关系数。corrcoef函数的调用格式为： (1) R=corrcoef(X，Y)%返回相关系数，其中X和Y是长度相等的向量 (2) R=corrcoef(A)%返回矩阵A的每列之间计算相关系数形成的相关系数矩阵 例如，求两个向量X和Y的相关系数，并求正态分布随机矩阵A的均值、标准方差和相关系数。 程序代码如下： >> X=［1,3,9,-2,7］; >> Y=［2,3,7,0,6］; >> r=corrcoef(X,Y)%求X和Y向量的相关系数 r = 1.00000.9985 0.99851.0000 >> A=randn(1000,3);%产生一个均值为0，方差为1的正态分布随机矩阵 >> y=mean(A)%计算矩阵A的列均值 y = 0.02530.00420.0427 >> D=std(A)%计算矩阵A的列标准方差 D = 0.99020.99191.0014 >> R=corrcoef(A)%计算A矩阵列的相关系数 R = 1.00000.0023-0.0028 0.00231.00000.0454 -0.00280.04541.0000 由上述结果可知，每列的均值接近0，每列的标准方差接近1，验证了A为标准正态分布随机矩阵。 5.5数值计算 数值计算是指利用计算机求解数学问题(比如，函数的零点、极值、积分和微分以及微分方程)近似解的方法。常用的数值分析有求函数的最小值、求过零点、数值微分、数值积分和解微分方程等。 5.5.1函数极值 数学上利用计算函数的导数来确定函数的最大和最小值点，然而，很多函数很难找到导数为零的点。为此，可以通过数值分析来确定函数的极值点。MATLAB只有求极小值的函数，没有专门求极大值的函数，因为f(x)的极大值问题等价于-f(x)的极小值问题。MATLAB求函数的极小值使用fminbnd和fminsearch函数。 1. 一元函数的极值 fminbnd函数可以获得一元函数在给定区间内的最小值，函数调用格式如下： (1) x=fminbnd(fun,x1,x2) 其中，fun是函数的句柄或匿名函数； x1和x2是寻找函数最小值的区间范围为x1> exam_5_16 x = 4.5150 y1 = -0.3975 图57在区间［0,5π］函数曲线 由图57可知，函数在x=4.5附近出现极小值点，极小值约为-0.4，验证了用极小值fminbnd函数求的极小值点和极小值是正确的。 2. 多元函数的极值 fminsearch函数可以获得多元函数的最小值，使用该函数时需要指定开始的初始值，获得初始值附近的局部最小值。该函数调用格式如下： (1) x = fminsearch(fun, x0) (2) ［x,y］= fminsearch(fun, x0) 其中，fun是多元函数的句柄或匿名函数； x0是给定的初始值； x是最小值的取值点； y是返回的最小值，可以省略。 【例517】使用fminsearch函数获取 f(x,y)二元函数在初始值(0, 0)附近的极小值，已知f(x,y)=100(y-x2)2+(1-x)2。 微课视频 在文件编辑窗口编写命令文件，保存为exam_5_17.m脚本文件。程序代码如下： clear fun=@(x)(100*(x(2)-x(1)^2)^2+(1-x(1))^2);%创建句柄函数 x0=［0,0］; ［x,y1］=fminsearch(fun,x0)%计算局部函数的极小值 程序运行结果如下： >> exam_5_17 x = 1.00001.0000 y1 = 3.6862e-10 由结果可知，由函数fminsearch计算出局部最小值点是［1, 1］，最小值为y1=3.6862e-10，和理论上是一致的。 5.5.2函数零点 一元函数f(x)过零点的求解相当于求解f(x)=0方程的根，MATLAB可以用fzero函数实现，使用该函数时需要指定一个初始值，在初始值附近查找函数值变号时的过零点，也可以根据指定区间来求过零点。该函数的调用格式为： (1) x=fzero(fun, x0) (2) ［x,y］=fzero(fun, x0) 其中，x为过零点的位置，如果找不到，则返回NaN； y是指函数在零点处函数的值； fun是函数句柄或匿名函数； x0是一个初始值或初始值区间。 需要指出，fzero函数只能返回一个局部零点，不能找出所有的零点，因此需要设定零点的范围。 微课视频 【例518】使用fzero函数求 f(x)=x2-5x+4分别在初始值x0=0，x0=5附近的过零点，并求出过零点函数的值。 在文件编辑窗口编写命令文件，保存为exam_5_18m脚本文件。程序代码如下： clear fun=@(x)(x^2-5*x+4);%创建句柄函数 x0=0; ［x,y1］=fzero(fun,x0)%求初始值x0为0附近，函数的过零点 x0=5; ［x,y1］=fzero(fun,x0)%求初始值x0为5附近，函数的过零点 x0=［0,3］; ［x,y1］=fzero(fun,x0)%求初始值x0区间内，函数的过零点 程序运行结果如下： >> exam_5_18 x = 1 y1 = 0 x = 4.0000 y1 = -3.5527e-15 x = 1 y1 = 0 由结果可知，用fzero函数可以求在初始值x0附近的函数过零点。不同的零点，需要设置不同的初始值x0。 5.5.3数值差分 任意函数f(x)在x点的前向差分定义为 Δf(x)=f(x+h)-f(x)(510) 称Δf(x)为函数f(x)在x点处以h(h>0)为步长的向前差分。 在MATLAB中，没有直接求数值导数的函数，只有计算前向差分的函数diff，其调用格式为： (1) D=diff(X)%计算向量X的向前差分，即D=X(i+1)-X(i),i=1,2,…n-1 (2) D=diff(X, n) %计算向量X的n阶向前差分。即diff(X, n)=diff(diff(X,n-1)) (3) D=diff(A, n, dim) %计算矩阵A的n阶差分。当dim=1(默认)，按行计算矩阵A的差分； %当dim=2，按列计算矩阵的差分 例如，已知矩阵A=163 624 581，分别求矩阵A行和列的一阶和二阶前向差分。 >> A=［1 6 3;6 2 4;5 8 1］ A = 163 624 581 >> D=diff(A,1,1) D = 5-41 -16-3 >> D=diff(A,1,2) D = 5-3 -42 3-7 >> D=diff(A,2,1) D = -610-4 >> D=diff(A,2,2) D = -8 6 -10 5.5.4数值积分 数值积分是研究定积分的数值求解的方法。MATLAB提供了很多种求数值积分的函数，主要包括一重积分和二重积分两类函数。 1. 一重数值积分 MATLAB提供了quad函数和quadl函数求一重定积分。调用格式为： (1) q=quad(fun, a, b, tol, trace) 它是一种采用自适应的Simpson方法的一重数值积分，其中fun为被积函数，是函数句柄； a和b为定积分的下限和上限； tol为绝对误差容限值，默认是10-6； trace控制是否展现积分过程，当trace取非0时,则展现积分过程，默认取0。 (2) q=quadl(fun, a, b, tol, trace) 它是一种采用自适应的Lobatto方法的一重数值积分，参数定义和quad一样。 微课视频 【例519】分别使用quad函数和quadl函数求 q=∫3π0e-0.2xsin(x)dx的数值积分。 在文件编辑窗口编写命令文件，保存为exam_5_19.m脚本文件。程序代码如下： clear fun=@(x)(exp(-0.2*x).*sin(x));%定义一个函数句柄 a=0;b=3*pi; q1=quad(fun,a,b)%自适应Simpson方法的数值积分 q2=quadl(fun,a,b)%自适应Lobatto方法的数值积分 q3=quad(fun,a,b,1e-3,1)%定义积分精度和显示积分过程 程序运行结果如下： >> exam_5_19 q1 = 1.1075 q2 = 1.1075 90.00000000002.55958120e+001.3793949196 110.00000000001.27979060e+000.6053358622 131.27979059931.27979060e+000.7742537042 152.55958119864.30561556e+00-0.6459997048 172.55958119862.15280778e+00-0.3430614927 194.71238898042.15280778e+00-0.3052258622 216.86519676222.55958120e+000.3762543321 q3 = 1.1076 其中，迭代过程最后一列的和为数值积分q3的值。 2. 多重数值积分 MATLAB提供了dblquad函数和triplequad函数求二重积分和三重积分。调用格式如下： (1) q2=dblquad(fun, xmin, xmax, ymin,ymax, tol) (2) q3=triplequad(fun, xmin, xmax, ymin, ymax, zmin, zmax, tol) 函数的参数定义和一重积分一样。 例如，求二重数值积分q=∫3π0∫2π0sin(x) y+xsin(y)dxdy。 代码如下： >> q=dblquad('sin(x)*y+x*sin(y)',0,2*pi,0,3*pi) q = 39.4784 5.5.5常微分方程求解 MATLAB为解常微分方程提供了多种数值求解方法，包括ode45、ode23、ode113、ode15s、ode23s、ode23t和ode23tb等函数，用得最多的是4/5阶龙格库塔法ode45函数。该函数格式如下： ［t,y］=ode45(fun,ts,y0,options) 其中，fun是待解微分方程的函数句柄； ts是自变量范围，可以是范围［t0, tf］，也可以是向量［t0,…,tf］； y0是初始值，y0和y是具有相同长度的列向量； options是设定微分方程解法器的参数，可以省略，也可以由odeset函数获得。 需要注意，用ode45求解时，需要将高阶微分方程 y(n)=f(t,y,y′,…,y(n-1))，改写为一阶微分方程组，通常解法是，假设y1=y，从而y1=y,y2=y′,…,yn=y(n-1)，于是高 阶微分方程可以转换为下述常微分方程组求解： y′1=y2 y′2=y3  y′n=f(t,y,y′,…,y(n-1)) (511) 微课视频 【例520】已知二阶微分方程 d2ydt2-3y′+2y=1, y(0)=1,dy(0)dt=0，t∈［0,1］，试用ode45函数解微分方程，作出y~t的关系曲线图。 (1) 首先把二阶微分方程改写为一阶微分方程组。 令y1=y,y2=y′1，则 dy1dt dy2dt=y2 3y2-2y1+1,y1(0) y2(0)=1 0(512) (2) 在文件编辑窗口编写命令文件，保存为exam_5_20.m脚本文件。程序代码如下： clear t0=［0,1］;%求解的时间区域 y0=［1;0］;%初值条件 ［t,y］=ode45(@f05_20,t0,y0);%采用ode45函数解微分方程 plot(t,y(:,1)) xlabel('t'),ylabel('y') title('y(t)-t') grid on 定义f05_20函数文件 function y= f05_20(t,y) %f05_20定义微分方程的函数文件 y=［y(2);3*y(2)-2*y(1)+1］; end 程序运行结果如图58所示。 图58二阶微分方程的数值解 习题 1. 已知多项式p1(x)=x4-3x3+5x+1，p2(x)=x3+2x2-6，求： (1) p(x)=p1(x)+p2(x)； (2) p(x)=p1(x)-p2(x)； (3) p(x)=p1(x)×p2(x)； (4) p(x)=p1(x)/p2(x)。 2. 已知多项式为p(x)=x4-2x2+4x-6，分别求x=3和x=［0,2,4,6,8］向量的多项式的值。 3. 已知多项式为p(x)=x4-2x2+4x-6，试求： (1) 用roots函数求该多项式的根r； (2) 用poly函数求根为r的多项式系数。 4. 已知两个多项式为p1(x)=x4-3x3+x+2，p2(x)=x3-2x2+4 (1) 求多项式p1(x)的导数； (2) 求两个多项式乘积p1(x)*p2(x)的导数； (3) 求两个多项式相除p2(x)/p1(x)的导数。 5. 已知分式表达式为f(s)=B(s)A(s)=s+1s2-7s+12。 (1) 求f(s)的部分分式展开式； (2) 将部分分式展开式转换为分式表达式。 6. 某电路元件，测试两端的电压U与流过的电流I的关系，实测数据见表53所示。用不同的插值方法(最接近点法、线性法、三次样条法和三次多项式法)计算I=9A处的电压U。 表53实测数据 流过的电流I/A024681012 两端的电压U/V0258.2121621 7. 某实验对一幅灰度图像灰度分布做测试。用i表示图像的宽度(PPI)，j表示图像的深度(PPI)，用I表示测得的各点图像颜色的灰度，测量结果如表54所示。 (1) 分别用最近点二维插值、三次样条插值、线性二维插值法求(13,12)点的灰度值； (2) 用三次多项式插值求图像宽度每1PPI，深度每1PPI处各点的灰度值，并用图形显示插值前后图像的灰度分布图。 表54图像各点颜色灰度测量值 j i 0510152025 0130132134133132131 5133137141138135133 10135138144143137134 15132134136135133132 8. 用polyfit函数实现一个3阶和5阶多项式在区间［0,2］内逼近函数f(x)=e-0.5x+sinx，利用绘图的方法，比较拟合的5阶多项式、7阶多项式和f(x)的区别。 9. 已知矩阵A=1047 962 394,试求： (1) 用max和min函数，求每行和每列的最大和最小元素，并求整个A的最大和最小元素； (2) 求矩阵A的每行和每列的平均值和中值； (3) 对矩阵A做各种排序； (4) 对矩阵A的各列和各行求和与求乘积； (5) 求矩阵A的行和列的标准方差； (6) 求矩阵A列元素的相关系数。 10. 已知y=e-0.5xsin(2*x)，在0≤x≤π区间内，使用fminbnd函数获取y函数的极小值。 11. 使用fzero函数求f(x)=x2-8x+12分别在初始值x=0，x=7附近的过零点，并求出过零点函数的值。 12. 已知矩阵A=1047 962 394，分别求矩阵A行和列的一阶和二阶前向差分。 13. 分别使用quad函数和quadl函数求 q=∫2π0sin(x)x+cos2xdx的数值积分。 14. 求二重数值积分q=∫2π0∫2π0xcos(y)+ysin(x)dxdy。 15. 已知二阶微分方程d2ydt2-2y′+y=0,y(0)=1,dy(0)dt=0，t∈［0,2］，试用ode45函数解微分方程，作出y~t的关系曲线图。 16. 洛伦兹(Lorenz)模型的状态方程表示为： dx1(t)dt=-βx1(t)+x2(t)x3(t) dx2(t)dt=-δx2(t)+δx3(t) dx3(t)dt=-x2(t)x1(t)+ρx2(t)-x3(t) , x1(0)=0 x2(0)=0 x3(0)=10-10 取δ=10,ρ=28,β=8/3，解该微分方程，并绘制出x1(t)~t时间曲线和x1(t)~x2(t)相空间曲线。