Matlab常用命令{2}

工程技术中的一些问题，如振动和稳定性问题，常归结为求一个方阵的特征值和特征向量。 特征值与特征向量的求法 设A为n阶方阵，如果数“ ”和n维列向量x使得关系式 成立，则称 为方阵A的特征值，非零向量x称为A对应于特征值“ ”的特征向量。 求矩阵 的特征值和特征向量 >>A=[-2 1 1;0 2 0;-4 1 3]; >>[V,D]=eig(A) V = -0.7071 -0.2425 0.3015 0 0 0.9045 -0.7071 -0.9701 0.3015 D = -1 0 0 0 2 0 0 0 2 即：特征值-1对应特征向量(-0.7071 0 -0.7071)T 特征值2对应特征向量(-0.2425 0 -0.9701)T和(-0.3015 0.9045 -0.3015)T 提高特征值的计算精度 [T,B] = balance(A) %求相似变换矩阵T和平衡矩阵B，满足 。 B = balance(A) %求平衡矩阵B 复对角矩阵转化为实对角矩阵 [V,D] = cdf2rdf (v,d) %将复对角阵d变为实对角阵D，在对角线上，用2×2实数块代替共轭复数对。 例1-91 >> A=[1 2 3;0 4 5;0 -5 4]; >> [v,d]=eig(A) v = 1.0000 -0.0191 - 0.4002i -0.0191 + 0.4002i 0 0 - 0.6479i 0 + 0.6479i 0 0.6479 0.6479 d = 1.0000 0 0 0 4.0000 + 5.0000i 0 0 0 4.0000 - 5.0000i >> [V,D]=cdf2rdf(v,d) V = 1.0000 -0.0191 -0.4002 0 0 -0.6479 0 0.6479 0 D = 1.0000 0 0 0 4.0000 5.0000 0 -5.0000 4.0000 正交基 B=orth(A) %将矩阵A正交规范化，B的列与A的列具有相同的空间，B的列向量是正交向量，且满足：B'*B = eye(rank(A))。 将矩阵 正交规范化。 >>A=[4 0 0; 0 3 1; 0 1 3]; >>B=orth(A) >>Q=B'*B 则显示结果为 P = 1.0000 0 0 0 0.7071 -0.7071 0 0.7071 0.7071 Q = 1.0000 0 0 0 1.0000 0 0 0 1.0000 二次型 求一个正交变换X=PY，把二次型 化成 标准 excel标准偏差excel标准偏差函数exl标准差函数国标检验抽样标准表免费下载红头文件格式标准下载 形。 解：先写出二次型的实对称矩阵 在Matlab编辑器中建立M文件如下： A=[0 1 1 -1;1 0 -1 1;1 -1 0 1;-1 1 1 0]; [P,D]=schur(A) syms y1 y2 y3 y4 y=[y1;y2;y3;y4]； X=vpa(P,2)*y %vpa表示可变精度计算，这里取2位精度 f=[y1 y2 y3 y4]*D*y 运行后结果显示如下： P = 780/989 780/3691 1/2 -390/1351 780/3691 780/989 -1/2 390/1351 780/1351 -780/1351 -1/2 390/1351 0 0 1/2 1170/1351 D = 1 0 0 0 0 1 0 0 0 0 -3 0 0 0 0 1 X = [ .79*y1+.21*y2+.50*y3-.29*y4] [ .21*y1+.79*y2-.50*y3+.29*y4] [ .56*y1-.56*y2-.50*y3+.29*y4] [ .50*y3+.85*y4] f = y1^2+y2^2-3*y3^2+y4^2 即 f = 矩阵和向量组的秩以及向量组的线性相关性 矩阵A的秩是矩阵A中最高阶非零子式的阶数；向量组的秩通常由该向量组构成的矩阵来计算。 rank(A) %返回矩阵A的行（或列）向量中线性无关个数 rank(A,tol) %tol为给定误差\不明白为什么向量的秩会与误差有关系? 求向量组 ， ， ， ， 的秩，并判断其线性相关性。 >>A=[1 -2 2 3;-2 4 -1 3;-1 2 0 3;0 6 2 3;2 -6 3 4]; >>k=rank(A) k = 3 由于秩为3 < 向量个数，因此向量组线性相关。 求行阶梯矩阵及向量组的基 行阶梯使用初等行变换，矩阵的初等行变换有三条： 1．交换两行 （第i、第j两行交换） 2．第i行的K倍 3．第i行的K倍加到第j行上去 通过这三条变换可以将矩阵化成行最简形，从而找出列向量组的一个最大无关组，Matlab将矩阵化成行最简形的命令是rref或rrefmovie。 rref(A) %用高斯—约当消元法和行主元法求A的行最简行矩阵R [R,jb] = rref(A) %jb是一个向量，其含义为：r = length(jb)为A的秩；A(:, jb)为A的列向量基；jb中元素表示基向量所在的列。 [R,jb] = rref(A,tol) %tol为指定的精度 rrefmovie(A) %给出每一步化简的过程 求向量组a1=(1,-2,2,3)，a2=(-2,4,-1,3)，a3=(-1,2,0,3)，a4=(0,6,2,3)，a5=(2,-6,3,4)的一个最大无关组。 >> a1=[1 -2 2 3]'; >>a2=[-2 4 -1 3]'; >>a3=[-1 2 0 3]'; >>a4=[0 6 2 3]'; >>a5=[2 -6 3 4]'; A=[a1 a2 a3 a4 a5] A = 1 -2 -1 0 2 -2 4 2 6 -6 2 -1 0 2 3 3 3 3 3 4 >> [R,jb]=rref(A) R = 1.0000 0 0.3333 0 1.7778 0 1.0000 0.6667 0 -0.1111 0 0 0 1.0000 -0.3333 0 0 0 0 0 jb = 1 2 4 >> A(:,jb) ans = 1 -2 0 -2 4 6 2 -1 2 3 3 3 即：a1 a2 a4为向量组的一个基。 稀疏矩阵的创建 sparse(A) %将矩阵A转化为稀疏矩阵形式，即由A的非零元素和下标构成稀疏矩阵S。若A本身为稀疏矩阵，则返回A本身。 sparse(m,n) %生成一个m×n的所有元素都是0的稀疏矩阵 sparse(i,j,s) %生成一个由长度相同的向量i，j和s定义的稀疏矩阵S，其中i，j是整数向量，定义稀疏矩阵的元素位置(i,j)，s是一个标量或与i，j长度相同的向量，表示在(i,j)位置上的元素。 sparse(i,j,s,m,n) %生成一个m×n的稀疏矩阵，(i,j)对应位置元素为si，m = max(i)且n =max(j)。 S = sparse(i,j,s,m,n,nzmax) %生成一个m×n的含有nzmax个非零元素的稀疏矩阵S，nzmax的值必须大于或者等于向量i和j的长度。 >> S=sparse(1:10,1:10,1:10) S = (1,1) 1 (2,2) 2 (3,3) 3 (4,4) 4 (5,5) 5 (6,6) 6 (7,7) 7 (8,8) 8 (9,9) 9 (10,10) 10 将稀疏矩阵转化为满矩阵 full(S) %S为稀疏矩阵，A为满矩阵。 >> S=sparse(1:5,1:5,4:8) S = (1,1) 4 (2,2) 5 (3,3) 6 (4,4) 7 (5,5) 8 >> A=full(S) A = 4 0 0 0 0 0 5 0 0 0 0 0 6 0 0 0 0 0 7 0 0 0 0 0 8 稀疏矩阵非零元素的索引 k = find(x) %按行检索X中非零元素的点，若没有非零元素，将返回空矩阵。 [i,j] = find(X) %检索X中非零元素的行标i和列标j [i,j,v] = find(X) %检索X中非零元素的行标i和列标j以及对应的元素值v 上例中 >> [i,j,v]=find(S) I =[ 1 2 3 4 5]’ j =[ 1 2 3 4 5]’ v =[4 5 6 7 8]’ 外部数据转化为稀疏矩阵 S=spconvert(D) %D是只有3列或4列的矩阵 说明：先运用load函数把外部数据（.mat文件或.dat文件）装载于MATLAB内存空间中的变量T；T数组的行维为nnz或nnz+1，列维为3（对实数而言）或列维为4（对复数而言）；T数组的每一行（以[i,j,Sre,Sim]形式）指定一个稀疏矩阵元素。 >> D=[1 2 3;2 5 4;3 4 6;3 6 7] D = 1 2 3 2 5 4 3 4 6 3 6 7 >> S=spconvert(D) S = (1,2) 3 (3,4) 6 (2,5) 4 (3,6) 7 >> D=[1 2 3 4; 2 5 4 0;3 4 6 9;3 6 7 4]; D = 1 2 3 4 2 5 4 0 3 4 6 9 3 6 7 4 >> S=spconvert(D) S = (1,2) 3.0000 + 4.0000i (3,4) 6.0000 + 9.0000i (2,5) 4.0000 (3,6) 7.0000 + 4.0000i 基本稀疏矩阵 带状（对角）稀疏矩阵 [B,d] = spdiags(A) %从矩阵A中提取所有非零对角元素，这些元素保存在矩阵B中，向量d表示非零元素的对角线位置。 B = spdiags(A,d) %从A中提取由d指定的对角线元素，并存放在B中。 A = spdiags(B,d,A) %用B中的列替换A中由d指定的对角线元素，输出稀疏矩阵。 A = spdiags(B,d,m,n) %产生一个m×n稀疏矩阵A，其元素是B中的列元素放 在由d指定的对角线位置上。 >>A = [11 0 13 0 0 22 0 24 0 0 33 0 41 0 0 44 0 52 0 0 0 0 63 0 0 0 0 74]; >>[B,d] = spdiags(A) B = 41 11 0 52 22 0 63 33 13 74 44 24 d = -3 %表示B的第1列元素在A中主对角线下方第3条对角线上 0 %表示B的第2列在A的主对角线上 2 %表示B的第3列在A的主对角线上方第2条对角线上 >> B=[1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16]; >> d=[-2 0 1 3]; >> A=spdiags(B,d,4,4); >> full(A) ans = 2 7 0 16 0 6 11 0 1 0 10 15 0 5 0 14 speye(m,n) %生成m×n的单位稀疏矩阵 speye(n) %生成n×n的单位稀疏矩阵 sprand(S) %生成与S具有相同稀疏结构的均匀分布随机矩阵 R = sprand(m,n,density) %生成一个m×n的服从均匀分布的随机稀疏矩阵，非零元素的分布密度是density。 R = sprand(m,n,density,rc) %生成一个近似的条件数为1/rc、大小为m×n的均匀分布的随机稀疏矩阵。 sprandn(S) %生成与S具有相同稀疏结构的正态分布随机矩阵。 sprandn(m,n,density) %生成一个m×n的服从正态分布的随机稀疏矩阵，非零元素的分布密度是density。 sprandn(m,n,density,rc) %生成一个近似的条件数为1/rc、大小为m×n的均匀分布的随机稀疏矩阵。 sprandsym(S) %生成稀疏对称随机矩阵，其下三角和对角线与S具有相同的结构，其元素服从均值为0、方差为1的标准正态分布。 R = sprandsym(n,density) %生成n×n的稀疏对称随机矩阵，矩阵元素服从正态分布，分布密度为density。 R = sprandsym(n,density,rc) %生成近似条件数为1/rc的稀疏对称随机矩阵 R = sprandsym(n,density,rc,kind) %生成一个正定矩阵，参数kind取值为kind=1表示矩阵由一正定对角矩阵经随机Jacobi旋转得到，其条件数正好为1/rc；kind=2表示矩阵为外积的换位和，其条件数近似等于1/rc；kind=3表示生成一个与矩阵S结构相同的稀疏随机矩阵，条件数近似为1/rc ，density被忽略。 稀疏矩阵的运算 nnz(X) %返回矩阵X中非零元素的个数 nonzeros(A) %返回矩阵A中非零元素按列顺序构成的列向量 >>A =[ 2 7 0 16 0 6 11 0 1 0 10 15 0 5 0 14]; >> s=nonzeros(A) s =[ 2 1 7 6 5 11 10 16 15 14]’ nzmax(S) %返回非零元素分配的内存总数n spalloc(m,n,nzmax) %产生一个m×n阶只有nzmax个非零元素的稀疏矩阵，这样可以有效减少存贮空间和提高运算速度。 spfun('function',S) %用S中非零元素对函数'function'求值，如果'function'不是对稀疏矩阵定义的，同样可以求值。 S = (1,1) 1 (2,2) 2 (3,3) 3 (4,4) 4 f= spfun('exp',S) %即指数e的非零元素方幂。 f= (1,1) 2.7183 (2,2) 7.3891 (3,3) 20.0855 (4,4) 54.5982 spones(S) %将稀疏矩阵S中的非零元素全换为1 spy(S) %画出稀疏矩阵S中非零元素的分布图形。S也可以是满矩阵。 spy(S,markersize) % markersize为整数，指定点阵大小。 spy(S,'LineSpec') %'LineSpec'指定绘图标记和颜色 spy(S,'LineSpec',markersize) %参数与上面相同 >> load west0479 >> A=west0479; >> spy(A,'ro',3) 结果如图1-4所示。 矩阵变换 colamd(S) %返回稀疏矩阵S的列的近似最小度排序向量p colmmd(S) %返回稀疏矩阵S的列的最小度排序向量p，按p排列后的矩阵为S(: , p)。 比较稀疏矩阵S与排序后的矩阵S(: , p) >> load west0479; >> S=west0479; >> p=colmmd(S); >> subplot(2,2,1),spy(S) >> subplot(2,2,2),spy(S(:,p)) >> subplot(2,2,3),spy(lu(S)) >> subplot(2,2,4),spy(lu(S(:,p))) 图1-5 稀疏矩阵的排序图 colperm(S) %返回一个稀疏矩阵S的列变换的向量。列按非0元素升序排列。有时这是LU分解前有用的变换：LU(S(:,j))。如果S是一个对称矩阵，对行和列进行排序，有利于Cholesky分解：chol(S(j,j)) dmperm (A) %返回A的行排列向量p，这样，如果A满列秩，就使得A(p,:)是具有非0对角线元素的方阵。 [p,q,r] = dmperm(A) %A为方阵，p为行排列向量，q为列排列向量，使得A(p,q) 是上三角块形式，r为索引向量。 [p,q,r,s] = dmperm(A) %A不是方阵，p,q,r含义与前面相同，s也是索引向量。 >> A=[11 0 13 0;41 0 0 44;0 22 0 24;0 0 63 0] A = 11 0 13 0 41 0 0 44 0 22 0 24 0 0 63 0 >> [p,q,r]=dmperm(A) p = 3 2 1 4 q = 2 4 1 3 r = 1 2 3 4 5 >> A(p,q) ans = 22 24 0 0 0 44 41 0 0 0 11 13 0 0 0 63 randperm (n) %对正整数1，2，3，…，n的随机排列，可以用来创建随机变换矩阵。 >> p=randperm(6) p = 3 4 6 5 1 2 p = symamd(S) %S为对称正定矩阵，返回排列向量p。 r = symrcm (S) %返回S的对称逆Cuthill-McKee排序r，使S的非0元素集中在主对角线附近。 p = symmmd(S) %返回S的对称最小度排列向量p，S为对称正定矩阵。 >> B = bucky+4*speye(60); >>r = symrcm(B); >>p = symmmd(B); >>R = B(r,r); >>S = B(p,p); >>subplot(2,2,1), spy(R), title('B(r,r)') >>subplot(2,2,2), spy(S), title('B(s,s)') >>subplot(2,2,3), spy(chol(R)), title('chol(B(r,r))') >>subplot(2,2,4), spy(chol(S)), title('chol(B(s,s))') 图1-6 稀疏对称最小度排列图 稀疏矩阵的近似欧几里得范数和条件数 c = condest(A) %计算方阵A的1-范数中条件数的下界值c [c,v] = condest(A) %方阵A的1-范数中条件数的下界值c和向量v，使得 ，即： norm(A*v,1) = norm(A,1)*norm(v,1)/c。 nrm = normest(S) %返回矩阵S的2-范数的估计值，相对误差为10-6。 nrm = normest(S,tol) %tol为指定的相对误差，而不用默认误差10-6。 [nrm,count] = normest(…) %count为给出的计算范数迭代的次数 稀疏矩阵的分解 不完全的LU分解,[L,U] = luinc(X,'0') %X为稀疏方阵；L为下三角矩阵的置换形式；U为上三角矩阵；'0'是一种分解标准。 [L,U,P] = luinc(X,'0') %L为下三角阵，其主对角线元素为1；U为上三角阵，p为单位矩阵的置换形式。 [L,U] = luinc(X,options) %options取值为：droptol表示指定的舍入误差；milu表示改变分解以便于从上三角分解因子中抽取被去掉的列元素。ugiag为1表示用droptol值代替上三角因子中的对角线上的零元素，默认值为0。thresh为中心临界值。 [L,U] = luinc(X,droptol) %droptol表示指定不完全分解的舍入误差 [L,U,P] = luinc(X,options) [L,U,P] = luinc(X,droptol) >> S=[11 0 13 0;41 0 0 44;0 22 0 24;0 0 63 0] S = 11 0 13 0 41 0 0 44 0 22 0 24 0 0 63 0 >> S=sparse(S) S = (1,1) 11 (2,1) 41 (3,2) 22 (1,3) 13 (4,3) 63 (2,4) 44 (3,4) 24 >> luinc(S,'0') ans = (1,1) 41.0000 (4,1) 0.2683 (2,2) 22.0000 (3,3) 63.0000 (4,3) 0.2063 (1,4) 44.0000 (2,4) 24.0000 >> [L,U,p]=luinc(S,'0') L = (1,1) 1.0000 (4,1) 0.2683 (2,2) 1.0000 (3,3) 1.0000 (4,3) 0.2063 (4,4) 1.0000 U = (1,1) 41 (2,2) 22 (3,3) 63 (1,4) 44 (2,4) 24 p = (4,1) 1 (1,2) 1 (2,3) 1 (3,4) 1 稀疏矩阵的不完全Cholesky分解 R = (X,droptol) %稀疏矩阵X的不完全Cholesky分解，droptol为指定误差。 R = cholinc(X,options) %options取值为：droptol表示舍入误差；michol表示如果michol=1，则从对角线上抽取出被去掉的元素。rdiag表示用droptol值代替上三角分解因子中的对角线上的零元素。 R = cholinc(X,'0') %'0'是一种分解标准 [R,p] = cholinc(X,'0') %不产生任何出错信息，如果R存在，则p=0；如果R不存在，则p为正整数。 R = cholinc(X,'inf') %进行Cholesky无穷分解 稀疏矩阵的特征值分解 d = eigs(A) %求稀疏矩阵A的6个最大特征值d，d以向量形式存放。 d = eigs(A,B) %求稀疏矩阵的广义特征值问题。满足AV=BVD，其中D为特征值对角阵，V为特征向量矩阵，B必须是对称正定阵或Hermitian正定阵。 d = eigs(A,k) %返回k个最大特征值 d = eigs(A,B,k) %返回k个最大特征值 d = eigs(A,k,sigma) %sigma取值：'lm' 表示最大数量的特征值；'sm' 最小数量特征值；对实对称问题：'la'表示最大特征值；'sa'为最小特征值；对非对称和复数问题：'lr' 表示最大实部；'sr' 表示最小实部；'li' 表示最大虚部；'si'表示最小虚部。 d = eigs(A,B,k,sigma) %同上 d = eigs(A,k,sigma,options) % options为指定参数：参见eigs帮助文件。 d = eigs(A,B,k,sigma,options) %同上。以下的参数k、sigma、options相同。 d = eigs(Afun,n) %用函数Afun代替A，n为A的阶数，D为特征值。 d = eigs(Afun,n,B) d = eigs(Afun,n,k) d = eigs(Afun,n,B,k) d = eigs(Afun,n,k,sigma) d = eigs(Afun,n,B,k,sigma) d = eigs(Afun,n,k,sigma,options) d = eigs(Afun,n,B,k,sigma,options) [V,D] = eigs(A,…) %D为6个最大特征值对角阵，V的列向量为对应特征向量。 [V,D] = eigs(Afun,n,…) [V,D,flag] = eigs(A,…) %flag表示特征值的收敛性，若flag=0，则所有特征值都收敛，否则，不是所有都收敛。 [V,D,flag] = eigs(Afun,n,…) 基本数学函数: Nchoosek:二项式系数或所有的组合数。该命令只有对n<15时有用。 C = nchoosek(n,k) %参量n,k为非负整数，返回n! / ( (n-k)! k!)，即一次从n个物体中取出k个的组合数。 C = nchoosek(v,k) %参量v为n维向量，返回一矩阵，其行向量的分量为一次性从v个物体中取k个物体的组合数。矩阵 C包含 =n! / ( (n-k)! k!)行与k列。 >>C = nchoosek(2:2:10,4) C = 2 4 6 8 2 4 6 10 2 4 8 10 2 6 8 10 4 6 8 10 s = rand('state') %返回一有35元素的列向量s，其中包含均匀分布生成器的当前状态。该改变生成器的当前的状态，见表2-1。 表2-1 命 令 含 义 Rand(’state’,s) 设置状态为s Rand(’state’,0) 设置生成器为初始状态 Rand(’state’,k) 设置生成器第k个状态(k为整数) Rand(’state’,sum(100*clock)) 设置生成器在每次使用时的状态都不同（因为clock每次都不同） >>R1 = rand(4,5) >>a = 10; b = 50; >>R2 = a + (b-a) * rand(5) % 生成元素均匀分布于(10,50)上的矩阵 计算结果可能为： R1 = 0.6655 0.0563 0.2656 0.5371 0.6797 0.3278 0.4402 0.9293 0.5457 0.6129 0.6325 0.4412 0.9343 0.9394 0.3940 0.5395 0.6501 0.5648 0.7084 0.2206 R2 = 33.6835 19.8216 36.9436 49.6289 46.4679 18.5164 34.2597 15.3663 31.0549 49.0377 19.0026 37.1006 33.6046 39.5361 13.9336 12.4641 12.9804 35.5420 23.2916 46.8304 28.5238 48.7418 49.0843 13.0512 10.9265 Y = randn(n) %返回n*n阶的方阵Y，其元素服从正态分布N(0,1)。若n不是一标量，则显示一出错信息。 Y = randn(m,n)、Y = randn([m n]) %返回阶数为m*n的，元素均匀分布于区间(0,1)上矩阵Y。 Y = randn(m,n,p,…)、Y = randn([m n p…]) %生成阶数m*n*p*…的，元素服从正态分布的多维随机阵列Y。 Y = randn(size(A)) %生成一与阵列A同型的随机正态阵列Y randn %该命令在每次单独使用时，都返回一随机数（服从正态分布）。 s = randn('state') %返回一有2元素的向量s，其中包含正态分布生成器的当前状态。该改变生成器的当前状态，见表2-2。 表2-2 命 令 含 义 randn(’state’,s) 设置状态为s randn(’state’,0) 设置生成器为初始状态 rand(’state’,k) 设置生成器第k个状态(k为整数) rand(’state’,sum(100*clock)) 设置生成器在每次使用时的状态都不同（因为clock每次都不同） >>R1 = rand(4,5) >>R2 = 0.6 + sqrt(0.1) * randn(5) 计算结果可能为： R1 = 0.2778 0.2681 0.5552 0.5167 0.8821 0.2745 0.3710 0.1916 0.3385 0.5823 0.9124 0.5129 0.4164 0.2993 0.0550 0.4125 0.2697 0.1508 0.9370 0.5878 R2 = 0.4632 0.9766 0.5410 0.6360 0.6931 0.0733 0.9760 0.8295 0.9373 0.1775 0.6396 0.5881 0.4140 0.6187 0.8259 0.6910 0.7035 1.2904 0.5698 1.1134 0.2375 0.6552 0.5569 0.3368 0.3812 插值法是实用的数值方法，是函数逼近的重要方法。在生产和科学实验中，自变量x与因变量y的函数y = f(x)的关系式有时不能直接写出表达式，而只能得到函数在若干个点的函数值或导数值。当要求知道观测点之外的函数值时，需要估计函数值在该点的值。 如何根据观测点的值，构造一个比较简单的函数y=φ(x)，使函数在观测点的值等于已知的数值或导数值。用简单函数y=φ(x)在点x处的值来估计未知函数y=f(x)在x点的值。寻找这样的函数φ(x)，办法是很多的。φ(x)可以是一个代数多项式，或是三角多项式，也可以是有理分式；φ(x)可以是任意光滑（任意阶导数连续）的函数或是分段函数。函数类的不同，自然地有不同的逼近效果。在许多应用中，通常要用一个解析函数（一、二元函数）来描述观测数据。 根据测量数据的类型： 1．测量值是准确的，没有误差。 2．测量值与真实值有误差。 这时对应地有两种处理观测数据方法： 1．插值或曲线拟合。 2．回归分析（假定数据测量是精确时，一般用插值法，否则用曲线拟合）。 MATLAB中提供了众多的数据处理命令。有插值命令，有拟合命令，有查表命令。 插值命令 interp1 一维数据插值（表格查找）。该命令对数据点之间计算内插值。它找出一元函数f(x)在中间点的数值。其中函数f(x)由所给数据决定。各个参量之间的关系示意图为图2-14。 图2-14 数据点与插值点关系示意图 yi = interp1(x,Y,xi) %返回插值向量yi，每一元素对应于参量xi，同时由向量x与Y的内插值决定。参量x指定数据Y的点。若Y为一矩阵，则按Y的每列计算。yi是阶数为length(xi)*size(Y,2)的输出矩阵。 yi = interp1(Y,xi) %假定x=1:N，其中N为向量Y的长度，或者为矩阵Y的行数。 yi = interp1(x,Y,xi,method) %用指定的算法计算插值： ’nearest’：最近邻点插值，直接完成计算； ’linear’：线性插值（缺省方式），直接完成计算； ’spline’：三次样条函数插值。对于该方法，命令interp1调用函数spline、ppval、mkpp、umkpp。这些命令生成一系列用于分段多项式操作的函数。命令spline用它们执行三次样条函数插值； ’pchip’：分段三次Hermite插值。对于该方法，命令interp1调用函数pchip，用于对向量x与y执行分段三次内插值。该方法保留单调性与数据的外形； ’cubic’：与’pchip’操作相同； ’v5cubic’：在MATLAB 5.0中的三次插值。 对于超出x范围的xi的分量，使用方法’nearest’、’linear’、’v5cubic’的插值算法，相应地将返回NaN。对其他的方法，interp1将对超出的分量执行外插值算法。 yi = interp1(x,Y,xi,method,'extrap') %对于超出x范围的xi中的分量将执行特殊的外插值法extrap。 yi = interp1(x,Y,xi,method,extrapval) %确定超出x范围的xi中的分量的外插值extrapval，其值通常取NaN或0。 >>x = 0:10; y = x.*sin(x); >>xx = 0:.25:10; yy = interp1(x,y,xx); >>plot(x,y,'kd',xx,yy) 插值图形为图2-15。 >> year = 1900:10:2010; >> product = [75.995 91.972 105.711 123.203 131.669 150.697 179.323 203.212 226.505 249.633 256.344 267.893 ]; >>p1995 = interp1(year,product,1995) >>x = 1900:1:2010; >>y = interp1(year,product,x,'pchip'); >>plot(year,product,'o',x,y) p1995 = 252.9885 插值图形为图2-16。 图2-15 一元函数插值图形 图2-16 离散数据的一维插值图 interp2 ，二维数据内插值（表格查找） ZI = interp2(X,Y,Z,XI,YI) %返回矩阵ZI，其元素包含对应于参量XI与YI（可以是向量、或同型矩阵）的元素，即Zi(i,j)←[Xi(i,j),yi(i,j)]。用户可以输入行向量和列向量Xi与Yi，此时，输出向量Zi与矩阵meshgrid(xi,yi)是同型的。同时取决于由输入矩阵X、Y与Z确定的二维函数Z=f(X,Y)。参量X与Y必须是单调的，且相同的划分格式，就像由命令meshgrid生成的一样。若Xi与Yi中有在X与Y范围之外的点，则相应地返回nan（Not a Number）。 ZI = interp2(Z,XI,YI) %缺省地，X=1:n、Y=1:m，其中[m,n]=size(Z)。再按第一种情形进行计算。 ZI = interp2(Z,n) %作n次递归计算，在Z的每两个元素之间插入它们的二维插值，这样，Z的阶数将不断增加。interp2(Z)等价于interp2(z,1)。 ZI = interp2(X,Y,Z,XI,YI,method) %用指定的算法method计算二维插值： ’linear’：双线性插值算法（缺省算法）； ’nearest’：最临近插值； ’spline’：三次样条插值； ’cubic’：双三次插值。 >>[X,Y] = meshgrid(-3:.25:3); >>Z = peaks(X,Y); >>[XI,YI] = meshgrid(-3:.125:3); >>ZZ = interp2(X,Y,Z,XI,YI); >>surfl(X,Y,Z);hold on; >>surfl(XI,YI,ZZ+15) >>axis([-3 3 -3 3 -5 20]);shading flat >>hold off 插值图形为图2-17。 >>years = 1950:10:1990; >>service = 10:10:30; >>wage = [150.697 199.592 187.625 179.323 195.072 250.287 203.212 179.092 322.767 226.505 153.706 426.730 249.633 120.281 598.243]; >>w = interp2(service,years,wage,15,1975) w =190.6288 interp3,三维数据插值（查表） 格式 VI = interp3(X,Y,Z,V,XI,YI,ZI) %找出由参量X,Y,Z决定的三元函数V=V(X,Y,Z)在点（XI,YI,ZI）的值。参量XI,YI,ZI是同型阵列或向量。若向量参量XI,YI,ZI是不同长度，不同方向（行或列）的向量，这时输出参量VI与Y1,Y2,Y3为同型矩阵。其中Y1,Y2,Y3为用命令meshgrid(XI,YI,ZI)生成的同型阵列。若插值点(XI,YI,ZI)中有位于点(X,Y,Z)之外的点，则相应地返回特殊变量值NaN。 VI = interp3(V,XI,YI,ZI) %缺省地，X=1:N，Y=1:M，Z=1:P，其中，[M,N,P]=size(V)，再按上面的情形计算。 VI = interp3(V,n) %作n次递归计算，在V的每两个元素之间插入它们的三维插值。这样，V的阶数将不断增加。interp3(V)等价于interp3(V,1)。 VI = interp3(…,method) %用指定的算法method作插值计算： ‘linear’：线性插值（缺省算法）； ‘cubic’：三次插值； ‘spline’：三次样条插值； ‘nearest’：最邻近插值。 说明 在所有的算法中，都要求X,Y,Z是单调且有相同的格点形式。当X,Y,Z是等距且单调时，用算法’*linear’，’*cubic’，’*nearest’，可得到快速插值。 >>[x,y,z,v] = flow(20); >>[xx,yy,zz] = meshgrid(.1:.25:10, -3:.25:3, -3:.25:3); >>vv = interp3(x,y,z,v,xx,yy,zz); >>slice(xx,yy,zz,vv,[6 9.5],[1 2],[-2 .2]); shading interp;colormap cool 插值图形为图2-18。 图2-18 三维插值图 interpft, 用快速Fourier算法作一维插值 y = interpft(x,n) %返回包含周期函数x在重采样的n个等距的点的插值y。若length(x)=m，且x有采样间隔dx，则新的y的采样间隔dy=dx*m/n。注意的是必须n≥m。若x为一矩阵，则按x的列进行计算。返回的矩阵y有与x相同的列数，但有n行。 y = interpft(x,n,dim) %沿着指定的方向dim进行计算 griddata,数据格点 ZI = griddata(x,y,z,XI,YI) %用二元函数z=f(x,y)的曲面拟合有不规则的数据向量x,y,z。griddata将返回曲面z在点（XI,YI）处的插值。曲面总是经过这些数据点（x,y,z）的。输入参量（XI,YI）通常是规则的格点（像用命令meshgrid生成的一样）。XI可以是一行向量，这时XI指定一有常数列向量的矩阵。类似地，YI可以是一列向量，它指定一有常数行向量的矩阵。 [XI,YI,ZI] = griddata(x,y,z,xi,yi) %返回的矩阵ZI含义同上，同时，返回的矩阵XI,YI是由行向量xi与列向量yi用命令meshgrid生成的。 […] = griddata(…,method) %用指定的算法method计算： ‘linear’：基于三角形的线性插值（缺省算法）； ‘cubic’： 基于三角形的三次插值； ‘nearest’：最邻近插值法； ‘v4’：MATLAB 4中的griddata算法。 spline,三次样条数据插值 yy = spline(x,y,xx) %对于给定的离散的测量数据x,y（称为断点），要寻找一个三项多项式 ，以逼近每对数据(x,y)点间的曲线。过两点 和 只能确定一条直线，而通过一点的三次多项式曲线有无穷多条。为使通过中间断点的三次多项式曲线具有唯一性，要增加两个条件（因为三次多项式有4个系数）： 1．三次多项式在点 处有： ； 2．三次多项式在点 处有： ； 3．p(x)在点 处的斜率是连续的（为了使三次多项式具有良好的解析性，加上的条件）； 4．p(x)在点 处的曲率是连续的； 对于第一个和最后一个多项式，人为地规定如下条件： ①． ②． 上述两个条件称为非结点(not-a-knot)条件。综合上述 内容 财务内部控制制度的内容财务内部控制制度的内容人员招聘与配置的内容项目成本控制的内容消防安全演练内容 ，可知对数据拟合的三次样条函数p(x)是一个分段的三次多项式： ，其中每段 都是三次多项式。 该命令用三次样条插值计算出由向量x与y确定的一元函数y=f(x)在点xx处的值。若参量y是一矩阵，则以y的每一列和x配对，再分别计算由它们确定的函数在点xx处的值。则yy是一阶数为length(xx)*size(y,2)的矩阵。 pp = spline(x,y) %返回由向量x与y确定的分段样条多项式的系数矩阵pp，它可用于命令ppval、unmkpp的计算。 对离散地分布在y=exp(x)sin(x)函数曲线上的数据点进行样条插值计算： >>x = [0 2 4 5 8 12 12.8 17.2 19.9 20]; y = exp(x).*sin(x); >>xx = 0:.25:20; >>yy = spline(x,y,xx); >>plot(x,y,'o',xx,yy) 插值图形结果为图2-19。 图2-19 三次样条插值 interpn, n维数据插值（查表） VI = interpn(X1,X2,,…,Xn,V,Y1,Y2,…,Yn) %返回由参量X1,X2,…,Xn,V确定的n元函数V=V(X1,X2,…,Xn)在点（Y1,Y2,…,Yn）处的插值。参量Y1,Y2,…,Yn是同型的矩阵或向量。若Y1,Y2,…,Yn是向量，则可以是不同长度，不同方向（行或列）的向量。它们将通过命令ndgrid生成同型的矩阵，再作计算。若点(Y1,Y2,…,Yn)中有位于点（X1,X2,…,Xn）之外的点，则相应地返回特殊变量NaN。 VI = interpn(V,Y1,Y2,…,Yn) %缺省地，X1=1:size(V,1)，X2=1:size(V,2)，…，Xn=1:size(V,n)，再按上面的情形计算。 VI = interpn(V,ntimes) %作ntimes次递归计算，在V的每两个元素之间插入它们的n维插值。这样，V的阶数将不断增加。interpn(V)等价于interpn(V,1)。 VI = interpn(…,method) %用指定的算法method计算： ‘linear’：线性插值（缺省算法）； ‘cubic’：三次插值； ‘spline’：三次样条插值法； ‘nearest’：最邻近插值算法。 meshgrid, 生成用于画三维图形的矩阵数据。 [X,Y] = meshgrid(x,y) 将由向量x，y（可以是不同方向的）指定的区域[min(x)，max(x)，min(y)，max(y)]用直线x=x(i),y=y(j)（i=1,2,…,length(x) ，j=1,2,…,length(y)）进行划分。这样，得到了length(x)*length(y)个点，这些点的横坐标用矩阵X表示，X的每个行向量与向量x相同；这些点的纵坐标用矩阵Y表示，Y的每个列向量与向量y相同。其中X,Y可用于计算二元函数z=f(x,y)与三维图形中xy平面矩形定义域的划分或曲面作图。 [X,Y] = meshgrid(x) %等价于[X,Y]=meshgrid(x,x)。 [X,Y,Z] = meshgrid(x,y,z) %生成三维阵列X,Y,Z，用于计算三元函数v=f(x,y,z)或三维容积图。 [X,Y] = meshgrid(1:3,10:14) 计算结果为： X = 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 Y = 10 10 10 11 11 11 12 12 12 13 13 13 14 14 14 ndgrid, 生成用于多维函数计算或多维插值用的阵列 [X1,X2,…,Xn] = ndgrid(x1,x2,…,xn) %把通过向量x1,x2,x3…,xn指定的区域转换为数组x1,x2,x3,…,xn。这样，得到了 length(x1)* length(x2)*…*length(xn)个点，这些点的第一维坐标用矩阵X1表示，X1的每个第一维向量与向量x1相同；这些点的第二维坐标用矩阵X2表示，X2的每个第二维向量与向量x2相同；如此等等。其中X1,X2,…,Xn可用于计算多元函数y=f(x1,x2,…,xn)以及多维插值命令用到的阵列。 [X1,X2,…,Xn] = ndgrid(x) %等价于[X1,X2,…,Xn] = ndgrid(x,x,…,x) 查表命令 table1(TAB,X0) %返回用表格矩阵TAB中的行线性插值元素，对X0（TAB的第一列查找X0）进行线性插值得到的结果Y。矩阵TAB第一列包含关键值，而其他列包含数据的矩阵。X0中的每一元素将相应地返回一线性插值行向量。矩阵TAB的第一列必须是单调的。 >>tab = [(1:4)' hilb(4)] >>y = table1(tab,[1 2.3 3.6 4]) 查表结果为： tab = 1.0000 1.0000 0.5000 0.3333 0.2500 2.0000 0.5000 0.3333 0.2500 0.2000 3.0000 0.3333 0.2500 0.2000 0.1667 4.0000 0.2500 0.2000 0.1667 0.1429 Warning: TABLE1 is obsolete and will be removed in future versions. Use INTERP1 or INTERP1Q instead. > In D:\MATLABR12\toolbox\matlab\polyfun\table1.m at line 31 y = 1.0000 0.5000 0.3333 0.2500 0.4500 0.3083 0.2350 0.1900 0.2833 0.2200 0.1800 0.1524 0.2500 0.2000 0.1667 0.1429 由上面结果可知，table1是一将要废弃的命令。 table2(TAB,X0,Y0) %返回用表格矩阵TAB中的行与列交叉线性线性插值元素，对X0（TAB的第一列查找X0）进行线性插值，对Y0（TAB的第一行查找Y0）进行线性插值，对上述两个数值进行交叉线性插值，得到的结果为Z。矩阵TAB是第一列与第一行列都包含关键值，而其他的元素包含数据的矩阵。TAB(1,1)的关键值将被忽略。[X0,Y0]中的每点将相应地返回一线性插值。矩阵TAB的第一行与第一列必须是单调的。 >>tab = [NaN 1:4; (1:4)' magic(4)] >>y = table2(tab,[2 3 3.7],[1.3 2.3 4]) 查表的结果为： tab = NaN 1 2 3 4 1 16 2 3 13 2 5 11 10 8 3 9 7 6 12 4 4 14 15 1 Warning: TABLE2 is obsolete and will be removed in future versions. Use INTERP2 instead. > In D:\MATLABR12\toolbox\matlab\polyfun\table2.m at line 24 Warning: TABLE1 is obsolete and will be removed in future versions. Use INTERP1 or INTERP1Q instead. > In D:\MATLABR12\toolbox\matlab\polyfun\table1.m at line 31 In D:\MATLABR12\toolbox\matlab\polyfun\table2.m at line 29 Warning: TABLE1 is obsolete and will be removed in future versions. Use INTERP1 or INTERP1Q instead. > In D:\MATLABR12\toolbox\matlab\polyfun\table1.m at line 31 In D:\MATLABR12\toolbox\matlab\polyfun\table2.m at line 31 y = 6.8000 10.7000 8.0000 8.4000 6.7000 12.0000 7.4200 12.0200 4.3000 由上面的结果可知，table2是将要废弃的命令。 数值积分 一元函数的数值积分 quad、quadl q = quad(fun,a,b,tol,trace,p1,p2,…) %将可选参数p1,p2,…等传递给函数fun(x,p1,p2,…)，再作数值积分。若tol=[]或trace=[]，则用缺省值进行计算。 [q,n] = quad(fun,a,b,…) %同时返回函数计算的次数n … = quadl(fun,a,b,…) %用高精度进行计算，效率可能比quad更好。 >>fun = inline(‘3*x.^2./(x.^3-2*x.^2+3)’); >>Q1 = quad(fun,0,2) >>Q2 = quadl(fun,0,2) 计算结果为： Q1 = 3.7224 Q2 = 3.7224 梯形法数值积分 T = trapz(Y) %用等距梯形法近似计算Y的积分。若Y是一向量，则trapz(Y)为Y的积分；若Y是一矩阵，则trapz(Y)为Y的每一列的积分；若Y是一多维阵列，则trapz(Y)沿着Y的第一个非单元集的方向进行计算。 T = trapz(X,Y) %用梯形法计算Y在X点上的积分。若X为一列向量，Y为矩阵，且size(Y,1) = length(X)，则trapz(X,Y)通过Y的第一个非单元集方向进行计算。 T = trapz(…,dim) %沿着dim指定的方向对Y进行积分。若参量中