第四章 MATLAB 的数值计算功能
Chapter 4,Numerical computation of MATLAB
一、多项式(Polynomial)`
1．多项式的表达与创建(Expression and Creating of polynomial)
(1) 多项式的表达（expression of polynomial）_
Matlab用行矢量表达多项式系数(Coefficient)，各元素按变量的降幂顺序排列，如多项式为：
P(x)=a0xn+a1xn-1+a2xn-2…an-1x+an
则其系数矢量(Vector of coefficient)为：P=[a0 a1 … an-1 an]
如将根矢量(Vector of root)表示为：
ar=[ ar1 ar2 … arn]
则根矢量与系数矢量之间关系为：
（x-ar1）(x- ar2) … (x- arn)= a0xn+a1xn-1+a2xn-2…an-1x+an
(2)多项式的创建（polynomial creating）
a)系数矢量的直接输入法利用poly2sym函数直接输入多项式的系数矢量，就可方便的建立符号形式的多项式。
例：创建多项式x3-4x2+3x+2
poly2sym([1 -4 3 2])
ans =
x^3-4*x^2+3*x+2
POLY Convert roots to polynomial.
POLY(A),when A is an N by N matrix,is a row vector with
N+1 elements which are the coefficients of the characteristic polynomial,DET(lambda*EYE(SIZE(A)) - A),
POLY(V),when V is a vector,is a vector whose elements are
the coefficients of the polynomial whose roots are the
elements of V,For vectors,ROOTS and POLY are inverse
functions of each other,up to ordering,scaling,and
roundoff error.
b) 由根矢量创建多项式通过调用函数 p=poly(ar)产生多项式的系数矢量,再利用poly2sym函数就可方便的建立符号形式的多项式。
注：（1）根矢量元素为n,则多项式系数矢量元素为n+1；
（2）函数poly2sym(pa) 把多项式系数矢量表达成符号形式的多项式，缺省情况下自变量符号为x，可以指定自变量。
（3）使用简单绘图函数ezplot可以直接绘制符号形式多项式的曲线。
例 1：由根矢量创建多项式。将多项式(x-6)(x-3)(x-8)表示为系数形式
a=[6 3 8] %根矢量
pa=poly(a) %求系数矢量
ppa=poly2sym(pa) %以符号形式表示原多项式
ezplot(ppa,[-50,50])
pa =
1 -17 90 -144
ppa =
x^3-17*x^2+90*x-144
注：含复数根的根矢量所创建的多项式要注意：
（1）要形成实系数多项式，根矢量中的复数根必须共轭成对；
（2）含复数根的根矢量所创建的多项式系数矢量中，可能带有很小的虚部，此时可采用取实部的命令(real)把虚部滤掉。
进行多项式的求根运算时，有两种方法，一是直接调用求根函数roots，poly和 roots 互为逆函数。另一种是先把多项式转化为伴随矩阵，然后再求其特征值，该特征值即是多项式的根。
例 3,由给定复数根矢量求多项式系数矢量。
r=[-0.5 -0.3+0.4i -0.3-0.4i];
p=poly(r)
pr=real(p)
ppr=poly2sym(pr)
p =
1.0000 1.1000 0.5500 0.1250
pr =
1.0000 1.1000 0.5500 0.1250
ppr =
x^3+11/10*x^2+11/20*x+1/8
c) 特征多项式输入法用poly函数可实现由矩阵的特征多项式系数创建多项式。
条件：特征多项式系数矢量的第一个元素必须为一。
例 2,求三阶方阵A的特征多项式系数，并转换为多项式形式。
a=[6 3 8;7 5 6; 1 3 5]
Pa=poly(a) %求矩阵的特征多项式系数矢量
Ppa=poly2sym(pa)
Pa =
1.0000 -16.0000 38.0000 -83.0000
Ppa =
x^3-17*x^2+90*x-144
注：n 阶方阵的特征多项式系数矢量一定是n +1阶的。
注：（1）要形成实系数多项式，根矢量中的复数根必须共轭成对；
（2）含复数根的根矢量所创建的多项式系数矢量中，可能带有很小的虚部，此时可采用取实部的命令(real)把虚部滤掉。
进行多项式的求根运算时，有两种方法，一是直接调用求根函数roots，poly和 roots 互为逆函数。另一种是先把多项式转化为伴随矩阵，然后再求其特征值，该特征值即是多项式的根。
例 4,将多项式的系数表示形式转换为根表现形式。
求 x3-6x2-72x-27的根
a=[1 -6 -72 -27]
r=roots(a)
r =
12.1229
-5.7345
-0.3884
MATLAB约定，多项式系数矢量用行矢量表示，根矢量用列矢量表示。
>>
1,多项式的乘除运算(Multiplication and division of polynomial)
多项式乘法用函数conv(a,b)实现,除法用函数deconv(a,b)实现。
例1：a(s)=s2+2s+3,b(s)=4s2+5s+6,计算 a(s)与 b(s)的乘积。
a=[1 2 3]; b=[4 5 6];
c=conv(a,b)
cs=poly2sym(c,’s’)
c =
4 13 28 27 18
cs =
4*s^4+13*s^3+28*s^2+27*s+18
例2,展开(s2+2s+2)(s+4)(s+1) （多个多项式相乘）
c=conv([1,2,2],conv([1,4],[1,1]))
cs=poly2sym(c,’s’) %（指定变量为s）
c =
1 7 16 18 8
cs =
s^4+7*s^3+16*s^2+18*s+8
例2：求多项式s^4+7*s^3+16*s^2+18*s+8分别被(s+4),(s+3)除后的结果。
c=[1 7 16 18 8];
[q1,r1]=deconv(c,[1,4]) %q—商矢量,r—余数矢量
[q2,r2]=deconv(c,[1,3])
cc=conv(q2,[1,3]) %对除(s+3)结果检验
test=((c-r2)==cc)
q1 =
1 3 4 2
r1 =
0 0 0 0 0
q2 =
1 4 4 6
r2 =
0 0 0 0 -10
cc =
1 7 16 18 18
test =
1 1 1 1 1
1,其他常用的多项式运算命令(Other computation command of polynomial)
pa=polyval(p,s) 按数组运算规则计算给定s时多项式p的值。
pm=polyvalm(p,s) 按矩阵运算规则计算给定s时多项式p的值。
[r,p,k]=residue(b,a) 部分分式展开，b,a分别是分子分母多项式系数矢量，r,p,k分别是留数、极点和直项矢量
p=polyfit(x,y,n) 用n阶多项式拟合x，y矢量给定的数据。
polyder(p) 多项式微分。
注,对于多项式b(s)与不重根的n阶多项式a(s)之比，其部分分式展开为：
式中：p1,p2,…,pn称为极点(poles)，r1,r2,…,rn 称为留数(residues)，k(s)称为直项(direct terms)，假如a(s)含有m重根pj，则相应部分应写成：
RESIDUE Partial-fraction expansion (residues).
[R,P,K] = RESIDUE(B,A) finds the residues,poles and direct term of a partial fraction expansion of the ratio of two polynomials B(s)/A(s),If there are no multiple roots,
B(s) R(1) R(2) R(n)
---- = -------- + -------- +,.,+ -------- + K(s)
A(s) s - P(1) s - P(2) s - P(n)
Vectors B and A specify the coefficients of the numerator and denominator polynomials in descending powers of s,The residues
are returned in the column vector R,the pole locations in column vector P,and the direct terms in row vector K,The number of poles is n = length(A)-1 = length(R) = length(P),The direct term coefficient vector is empty if length(B) < length(A),otherwise
length(K) = length(B)-length(A)+1.
If P(j) =,.,= P(j+m-1) is a pole of multplicity m,then the expansion includes terms of the form
R(j) R(j+1) R(j+m-1)
-------- + ------------ +,.,+ ------------
s - P(j) (s - P(j))^2 (s - P(j))^m
[B,A] = RESIDUE(R,P,K),with 3 input arguments and 2 output arguments,converts the partial fraction expansion back to the polynomials with coefficients in B and A.
例3：对 (3x4+2x3+5x2+4x+6)/(x5+3x4+4x3+2x2+7x+2) 做部分分式展开
a=[1 3 4 2 7 2];
b=[3 2 5 4 6];
[r,s,k]=residue(b,a)
r =
1.1274 + 1.1513i
1.1274 - 1.1513i
-0.0232 - 0.0722i
-0.0232 + 0.0722i
0.7916
s =
-1.7680 + 1.2673i
-1.7680 - 1.2673i
0.4176 + 1.1130i
0.4176 - 1.1130i
-0.2991
k =
[] (分母阶数高于分子阶数时，k将是空矩阵，表示无此项)
例5：对一组实验数据进行多项式最小二乘拟合(least square fit)
x=[1 2 3 4 5]; % 实验数据
y=[5.5 43.1 128 290.7 498.4];
p=polyfit(x,y,3) %做三阶多项式拟合
x2=1:.1:5;
y2=polyval(p,x2); % 根据给定值计算多项式结果
plot(x,y,’o’,x2,y2)
线性代数(Linear Algebra)
解线性方程(Linear equation)就是找出是否存在一个唯一的矩阵x，使得a,b满足关系：
ax=b 或 xa=b
MALAB中x=a\b 是方程 ax=b 的解,x=b/a是方程式xa=b的解。
通常线性方程多写成ax=b，“\”较多用，两者的关系为：
(b/a)’=(a’\b’)
系数矩阵a可能是m行n列的，有三种情况：
*方阵系统,(Square matrix) m=n 可求出精确解（a必须是非奇异(nonsingular)，即满秩(full rank)）
*超定系统：(Overdetermind system)m>n 可求出最小二乘解
*欠定系统：(Underdetermind system) m<n 可尝试找出含有最少m个基解或最小范数解
MATLAB对不同形式的参数矩阵，采用不同的运算法则来处理，它会自动检测参数矩阵，以区别下面几种形式：
*三角矩阵(Triangular Matrix)
*对称正定矩阵(symmetrical positive determined matrix)
*非奇异方阵(Nonsingular matrix)
*超定系统(Overdetermind system)
*欠定系统(Underdetermind system)
方阵系统：(Square array)
最常见的是系数矩阵为方阵a，常数项b为列矢量,其解x可写成x=a\b,x和b大小相同。
例1,求方阵系统的根。
a=[11 6 7; 5 13 9; 17 1 8]
b=[16 13 4]’
x=a\b
a =
11 6 7
5 13 9
17 1 8
b =
16
13
4
x =
3.9763
5.4455
-8.6303
例2：假如a,b为两个大小相同的矩阵，求方阵系统的根。
a=[4 5 9; 18 19 5; 1 4 13]
b=[1 5 12; 3 15 19; 7 6 10]
x=a\b
C=a*x
a =
4 5 9
18 19 5
1 4 13
b =
1 5 12
3 15 19
7 6 10
x =
-3.6750 -0.7333 2.9708
3.7250 1.4667 -2.1292
-0.3250 0.0667 1.1958
C =
1.0000 5.0000 12.0000
3.0000 15.0000 19.0000
7.0000 6.0000 10.0000
若方阵a的各个行矢量线性相关(linear correlation)，则称方阵a为奇异矩阵。这时线性方程将有无穷多组解。若方阵是奇异矩阵，则反斜线运算因子将发出警告信息。
2．超定系统(Overdetermind system)
实验数据较多，寻求他们的曲线拟合。
如在t内测得一组数据y：
t y
0.0 0.82
0.3 0.72
0.8 0.63
0.60
1.6 0.55
2.2 0.50
这些数据显然有衰减指数趋势,y(t)~c1+c2e-t
此方程意为y矢量可以由两个矢量逐步逼近而得，一个是单行的常数矢量，一个是由指数e-t项构成，两个参数c1和c2可用最小二乘法求得，它们表示实验数据与方程y(t)~c1+c2e-t之间距离的最小平方和。
例1,求上述数据的最小二乘解。将数据带入方程式y(t)~c1+c2e-t中，可得到含有两个未知数的6个等式，可写成6行2 列的矩阵e.
t=[0 0.3 0.8 1.1 1.6 2.2]’;
y=[0.82 0.72 0.63 0.60 0.55 0.50]’;
e=[ones(size(t)) exp(-t)] %求6个y(t)方程的系数矩阵
c=e\y % 求方程的解
e =
1.0000 1.0000
1.0000 0.7408
1.0000 0.4493
1.0000 0.3329
1.0000 0.2019
1.0000 0.1108
c =
0.4744
0.3434
带入方程得：y(t)~0.4744+0.3434e-t
用此方程可绘制曲线：
t=[0 0.3 0.8 1.1 1.6 2.2]’;
y=[0.82 0.72 0.63 0.60 0.55 0.50]’;
t1=[0:0.1:2.5]’; y1=[ones(size(t1)),exp(-t1)]*c
plot(t1,y1,’b’,t,y,’ro’)
如果一个矩阵的行矢量是线性相关的，则它的最小二乘解并不唯一，因此，a\b运算将给出警告，并产生含有最少元素的基解。
3,欠定系统,(Underdetermind system)
欠定系统为线性相关系统，其解都不唯一，MATLAB会计算一组构成通解的基解，而方程的特解则用QR分解法决定。
两种解法：最少元素解a\b，最小范数解pinv(a)*b.
例,用两种方法求解欠定系统。
对a和矢量b分别用a\b和pinv(a)*b求解:
a=[1 1 1; 1 1 -1]
b=[10 6]’
p=a\b
q=pinv(a)*b
a =
1 1 1
1 1 -1
b =
10
6
p =
8.0000
0
2.0000
q =
4.0000
4.0000
2.0000
三,逆矩阵及行列式(Revers and determinant of matrix)
1,方阵的逆和行列式(Revers and determinant of square matrix)
若a是方阵，且为非奇异阵，则方程ax=I和 xa=I有相同的解X。X称为a的逆矩阵，记做a-1，在MATLAB中 用inv 函数来计算矩阵的逆。计算方阵的行列式则用det函数。
DET Determinant.
DET(X) is the determinant of the square matrix X.
Use COND instead of DET to test for matrix singularity.
INV Matrix inverse.
INV(X) is the inverse of the square matrix X,A warning message is printed if X is badly scaled or nearly singular.
例：计算方阵的行列式和逆矩阵。
a=[3 -3 1;-3 5 -2;1 -2 1];
b=[14 13 5; 5 1 12;6 14 5];
d1=det(a)
x1=inv(a)
d2=det(b)
x2=inv(b)
d1 =
1
x1 =
1.0000 1.0000 1.0000
1.0000 2.0000 3.0000
1.0000 3.0000 6.0000
d2 =
-1351
x2 =
0.1207 -0.0037 -0.1118
-0.0348 -0.0296 0.1058
-0.0474 0.0873 0.0377
2,广义逆矩阵（伪逆）(Generalized inverse matrix)
一般非方阵无逆矩阵和行列式，方程ax=I 和xa=I至少有一个无解，这种矩阵可以求得特殊的逆矩阵，成为广义逆矩阵(generalized inverse matrix)（或伪逆 pseudoinverse）。矩阵amn存在广义逆矩阵xnm，使得 ax=Imn,MATLAB用pinv函数来计算广义逆矩阵。
例：计算广义逆矩阵。
a=[8 14; 1 3; 9 6]
x=pinv(a)
b=x*a
c=a*x
d=c*a %d=a*x*a=a
e=x*c %e=x*a*x=x
a =
8 14
1 3
9 6
x =
-0.0661 -0.0402 0.1743
0.1045 0.0406 -0.0974
b =
1.0000 -0.0000
-0.0000 1.0000
c =
0.9334 0.2472 0.0317
0.2472 0.0817 -0.1177
0.0317 -0.1177 0.9849
d =
8.0000 14.0000
1.0000 3.0000
9.0000 6.0000
e =
-0.0661 -0.0402 0.1743
0.1045 0.0406 -0.0974
PINV Pseudoinverse.
X = PINV(A) produces a matrix X of the same dimensions as A' so that A*X*A = A,X*A*X = X and A*X and X*A are Hermitian,The computation is based on SVD(A) and any singular values less than a tolerance are treated as zero.
The default tolerance is MAX(SIZE(A)) * NORM(A) * EPS.
PINV(A,TOL) uses the tolerance TOL instead of the default.
四,矩阵分解(Matrix decomposition)
MATLAB求解线性方程的过程基于三种分解法则：
（１）Cholesky分解，针对对称正定矩阵；
（２）高斯消元法,针对一般矩阵；
（３）正交化,针对一般矩阵（行数≠列数）
这三种分解运算分别由chol,lu和 qr三个函数来分解.
Cholesky分解(Cholesky Decomposition)
仅适用于对称和上三角矩阵例：cholesky分解。
a=pascal(6)
b=chol(a)
a =
1 1 1 1 1 1
1 2 3 4 5 6
1 3 6 10 15 21
1 4 10 20 35 56
1 5 15 35 70 126
1 6 21 56 126 252
b =
1 1 1 1 1 1
0 1 2 3 4 5
0 0 1 3 6 10
0 0 0 1 4 10
0 0 0 0 1 5
0 0 0 0 0 1
CHOL Cholesky factorization.
CHOL(X) uses only the diagonal and upper triangle of X,The lower triangular is assumed to be the (complex conjugate) transpose of the upper,If X is positive definite,then R = CHOL(X) produces an upper triangular R so that R'*R = X,If X is not positive definite,an error message is printed.
[R,p] = CHOL(X),with two output arguments,never produces an
error message,If X is positive definite,then p is 0 and R is the same as above,But if X is not positive definite,then p is a positive integer.
When X is full,R is an upper triangular matrix of order q = p-1
so that R'*R = X(1:q,1:q),When X is sparse,R is an upper triangular matrix of size q-by-n so that the L-shaped region of the first q rows and first q columns of R'*R agree with those of X.
2,LU分解(LU factorization).
用lu函数完成LU分解，将矩阵分解为上、下两个三角阵，其调用格式为：
[l,u]=lu(a) l代表下三角阵，u代表上三角阵。
例：
LU分解。
a=[47 24 22; 11 44 0;30 38 41]
[l,u]=lu(a)
a =
47 24 22
11 44 0
30 38 41
l =
1.0000 0 0
0.2340 1.0000 0
0.6383 0.5909 1.0000
u =
47.0000 24.0000 22.0000
0 38.3830 -5.1489
0 0 30.0000
LU LU factorization.
[L,U] = LU(X) stores an upper triangular matrix in U and a "psychologically lower triangular matrix" (i.e,a product of lower triangular and permutation matrices) in L,so that X = L*U,X can be rectangular,
[L,U,P] = LU(X) returns unit lower triangular matrix L,upper triangular matrix U,and permutation matrix P so that P*X = L*U.
3,QR分解(Orthogonal-triangular decomposition).
函数调用格式：[q,r]=qr(a),q代表正规正交矩阵，r代表三角形矩阵。原始阵a不必一定是方阵。如果矩阵a是m×n阶的，则矩阵q是m×m阶的，矩阵r是m×n阶的。
例：QR分解.
A=[22 46 20 20; 30 36 46 44;39 8 45 2];
[q,r]=qr(A)
q =
-0.4082 -0.7209 -0.5601
-0.5566 -0.2898 0.7786
-0.7236 0.6296 -0.2829
r =
-53.8981 -44.6027 -66.3289 -34.1014
0 -38.5564 0.5823 -25.9097
0 0 11.8800 22.4896
QR Orthogonal-triangular decomposition.
[Q,R] = QR(A) produces an upper triangular matrix R of the same
dimension as A and a unitary matrix Q so that A = Q*R.
[Q,R,E] = QR(A) produces a permutation matrix E,an upper
triangular R and a unitary Q so that A*E = Q*R,The column
permutation E is chosen so that abs(diag(R)) is decreasing,
[Q,R] = QR(A,0) produces the "economy size" decomposition,If A is m-by-n with m > n,then only the first n columns of Q are computed.
4,特征值与特征矢量(Eigenvalues and eigenvectors).
MATLAB中使用函数eig计算特征值和 特征矢量，有两种调用方法：
*e=eig(a),其中e是包含特征值的矢量；
*[v,d]=eig(a),其中v是一个与a相同的n×n阶矩阵，它的每一列是矩阵a的一个特征值所对应的特征矢量，d为对角阵，其对角元素即为矩阵a的特征值。
例：计算特征值和特征矢量。
a=[34 25 15; 18 35 9; 41 21 9]
e=eig(a)
[v,d]=eig(a)
a =
34 25 15
18 35 9
41 21 9
e =
68.5066
15.5122
-6.0187
v =
-0.6227 -0.4409 -0.3105
-0.4969 0.6786 -0.0717
-0.6044 -0.5875 0.9479
d =
68.5066 0 0
0 15.5122 0
0 0 -6.0187
EIG Eigenvalues and eigenvectors.
E = EIG(X) is a vector containing the eigenvalues of a square matrix X.
[V,D] = EIG(X) produces a diagonal matrix D of eigenvalues and a full matrix V whose columns are the corresponding eigenvectors so that X*V = V*D.
[V,D] = EIG(X,'nobalance') performs the computation with balancing
disabled,which sometimes gives more accurate results for certain
problems with unusual scaling,If X is symmetric,EIG(X,'nobalance')
is ignored since X is already balanced.
5,奇异值分解.( Singular value decomposition).
如存在两个矢量u,v及一常数c,使得矩阵A满足：Av=cu,A’u=cv
称c为奇异值，称u,v为奇异矢量。
将奇异值写成对角方阵∑，而相对应的奇异矢量作为列矢量则可写成两个正交矩阵U，V,使得,AV=U∑,A‘U=V∑ 因为U，V正交，所以可得奇异值表达式：
A=U∑V’。
一个m行n列的矩阵A经奇异值分解，可求得m行m列的U,m行n列的矩阵∑和n行n列的矩阵V.。
奇异值分解用svd函数实现，调用格式为；
[u,s,v]=svd(a)
SVD Singular value decomposition.
[U,S,V] = SVD(X) produces a diagonal matrix S,of the same dimension as X and with nonnegative diagonal elements in decreasing order,and unitary matrices U and V so that X = U*S*V'.
S = SVD(X) returns a vector containing the singular values.
[U,S,V] = SVD(X,0) produces the "economy size" decomposition,If X is m-by-n with m > n,then only the first n columns of U are computed and S is n-by-n.
例,奇异值分解。
a=[8 5; 7 3;4 6];
[u,s,v]=svd(a) % s为奇异值对角方阵
u =
-0.6841 -0.1826 -0.7061
-0.5407 -0.5228 0.6591
-0.4895 0.8327 0.2589
s =
13.7649 0
0 3.0865
0 0
v =
-0.8148 -0.5797
-0.5797 0.8148
五,数据分析(Data Analyaia)
MATLAB对数据分析有两条约定：
（1） 若输入量X是矢量，则不论是行矢量还是列矢量，运算是对整个矢量进行的；
（2）若输入量X是数组，（或称矩阵），则命令运算是按列进行的。即默认每个列是有一个变量的不同“观察“所得的数据组成。
1,基本统计命令 （表4-1）
例,做各种基本统计运算。
A=[5 -10 -6 0;2 6 3 -3;-9 5 -10 11;-22 17 0 -19;-1 6 -4 4]
Amax=max(A) %找A各列的最大元素
Amin=min(A) %找A各列的最小元素
Amed=median(A) %找A各列的中位元素
Amean=mean(A) %找A各列的平均值
Astd=std(A) %求A各列的标准差
Aprod=prod(A) %求A各列元素的积
Asum=sum(A) %求A各列元素的和
S=cumsum(A) %求A各列元素的累积和
P=cumprod(A) %求A各列元素的累积j积
I=sort(A) %使A的各列元素按递增排列
A =
5 -10 -6 0
2 6 3 -3
-9 5 -10 11
-22 17 0 -19
-1 6 -4 4
Amax =
5 17 3 11
Amin =
-22 -10 -10 -19
Amed =
-1 6 -4 0
Amean =
-5.0000 4.8000 -3.4000 -1.4000
Astd =
10.8397 9.6281 5.0794 11.1490
Aprod =
-1980 -30600 0 0
Asum =
-25 24 -17 -7
S =
5 -10 -6 0
7 -4 -3 -3
-2 1 -13 8
-24 18 -13 -11
-25 24 -17 -7
P =
5 -10 -6 0
10 -60 -18 0
-90 -300 180 0
1980 -5100 0 0
-1980 -30600 0 0
I =
-22 -10 -10 -19
-9 5 -6 -3
-1 6 -4 0
2 6 0 4
5 17 3 11
>>
求矩阵元素的最大值、最小值可用：
Amax=max(maxA)) 或 Amax=max(A(:))，
Amin=min(min(A)) 或 Amin=min(A(:))
2．协方差阵和相关阵(Covariance matrix and Correlation coefficients).
（表 4—2）
例,计算协方差和相关阵。
x=rand(10,3);
y=rand(10,3);
cx=cov(x) %求协方差阵
cy=cov(y)
cxy=cov(x,y) %求两随机变量的协方差
px=corrcoef(x) %求相关阵
pxy=corrcoef(x,y) %求两随机变量的（2×2）相关系数
cx =
0.0483 -0.0066 0.0146
-0.0066 0.0283 0.0154
0.0146 0.0154 0.0978
cy =
0.1177 0.0073 -0.0127
0.0073 0.0239 -0.0230
-0.0127 -0.0230 0.0772
cxy =
0.0550 0.0023
0.0023 0.0697
px =
1.0000 -0.1783 0.2118
-0.1783 1.0000 0.2934
0.2118 0.2934 1.0000
pxy =
1.0000 0.0372
0.0372 1.0000
COV Covariance matrix.
COV(X),if X is a vector,returns the variance,For matrices,where each row is an observation,and each column a variable,COV(X) is the covariance matrix,DIAG(COV(X)) is a vector of variances for each column,and SQRT(DIAG(COV(X))) is a vector of standard deviations,
COV(X,Y),where X and Y are vectors of equal length,is equivalent to COV([X(:) Y(:)]),
COV(X) or COV(X,Y) normalizes by (N-1) where N is the number of observations,This makes COV(X) the best unbiased estimate of the
covariance matrix if the observations are from a normal distribution.
CORRCOEF Correlation coefficients.
R=CORRCOEF(X) calculates a matrix R of correlation coefficients for an array X,in which each row is an observation and each column is a variable.
R=CORRCOEF(X,Y),where X and Y are column vectors,is the same as
R=CORRCOEF([X Y]),If C is the covariance matrix,C = COV(X),then CORRCOEF(X) is the matrix whose (i,j)'th element is
C(i,j)/SQRT(C(i,i)*C(j,j)).
2,微分与梯度 (Difference and approximate derivative,gradient).
(表4—3)
例1：按列求微分。
x=[1,10,20;2,12,23;3,14,26;3,16,29]
d=diff(x) %求一阶微分
x =
1 10 20
2 12 23
3 14 26
3 16 29
d =
1 2 3
1 2 3
0 2 3
例2,对于（u=x2+y2和Δ2=4）求5点差分。
[x,y]=meshgrid(-4:4,-3:3);
u=x.^2+y.^2
v4=4*del2(u) %求m×n阶矩阵U的五点差分矩阵
u =
25 18 13 10 9 10 13 18 25
20 13 8 5 4 5 8 13 20
17 10 5 2 1 2 5 10 17
16 9 4 1 0 1 4 9 16
17 10 5 2 1 2 5 10 17
20 13 8 5 4 5 8 13 20
25 18 13 10 9 10 13 18 25
v4 =
4 4 4 4 4 4 4 4 4
4 4 4 4 4 4 4 4 4
4 4 4 4 4 4 4 4 4
4 4 4 4 4 4 4 4 4
4 4 4 4 4 4 4 4 4
4 4 4 4 4 4 4 4 4
4 4 4 4 4 4 4 4 4
MESHGRID X and Y arrays for 3-D plots.
[X,Y] = MESHGRID(x,y) transforms the domain specified by vectors
x and y into arrays X and Y that can be used for the evaluation
of functions of two variables and 3-D surface plots.
The rows of the output array X are copies of the vector x and
the columns of the output array Y are copies of the vector y.
[X,Y] = MESHGRID(x) is an abbreviation for [X,Y] = MESHGRID(x,x).
[X,Y,Z] = MESHGRID(x,y,z) produces 3-D arrays that can be used to
evaluate functions of three variables and 3-D volumetric plots.
DEL2 Discrete Laplacian.
L = DEL2(U) when U is a matrix,is an discrete approximation of 0.25*del^2 u = (d^2u/dx^2 + d^2/dy^2)/4,The matrix L is the same size as U with each element equal to the difference between an element of U and the average of its four neighbors.
L = DEL2(U) when U is an N-D array,returns an approximation of
(del^2 u)/2/n where n is ndims(u).
L = DEL2(U,H),where H is a scalar,uses H as the spacing between
points in each direction (H=1 by default).
L = DEL2(U,HX,HY) when U is 2-D,uses the spacing specified by HX
and HY,If HX is a scalar,it gives the spacing between points in
the x-direction,If HX is a vector,it must be of length SIZE(U,2)
and specifies the x-coordinates of the points,Similarly,if HY
is a scalar,it gives the spacing between points in the
y-direction,If HY is a vector,it must be of length SIZE(U,1) and
specifies the y-coordinates of the points.
L = DEL2(U,HX,HY,HZ,...) when U is N-D,uses the spacing given by
HX,HY,HZ,etc.
例3：产生一个二元函数偏导数和梯度。
x=-2:0.2:2;
y=-2:0.2:2;
[xx,yy]=meshgrid(x,y);
z=xx.*exp(-xx.^2-yy.^2);
[Gx,Gy]=gradient(z,0.2,0.2); %Gx,Gy分别是二元函数的偏导
contour(x,y,z,'k'),hold on,
quiver(xx,yy,Gx,Gy,'r'),hold off
DIFF Difference and approximate derivative.
DIFF(X),for a vector X,is [X(2)-X(1) X(3)-X(2),.,X(n)-X(n-1)].
DIFF(X),for a matrix X,is the matrix of row differences,
[X(2:n,:) - X(1:n-1,:)].
DIFF(X),for an N-D array X,is the difference along the first
non-singleton dimension of X.
DIFF(X,N) is the N-th order difference along the first non-singleton dimension (denote it by DIM),If N >= size(X,DIM),DIFF takes successive differences along the next non-singleton dimension.
（[FX,FY] = GRADIENT(F,HX,HY),when F is 2-D,uses the spacing
specified by HX and HY,HX and HY can either be scalars to specify
the spacing between coordinates or vectors to specify the coordinates of the points,If HX and HY are vectors,their length must match the corresponding dimension of F.）
（QUIVER(X,Y,U,V) plots velocity vectors as arrows with components (u,v)
at the points (x,y),The matrices X,Y,U,V must all be the same size
and contain corresponding position and velocity components (X and Y
can also be vectors to specify a uniform grid),QUIVER automatically
scales the arrows to fit within the grid.）
GRADIENT Approximate gradient.
[FX,FY] = GRADIENT(F) returns the numerical gradient of the
matrix F,FX corresponds to dF/dx,the differences in the
x (column) direction,FY corresponds to dF/dy,the differences
in the y (row) direction,The spacing between points in each
direction is assumed to be one,When F is a vector,DF = GRADIENT(F) is the 1-D gradient.
[FX,FY] = GRADIENT(F,H),where H is a scalar,uses H as the
spacing between points in each direction.
[FX,FY] = GRADIENT(F,HX,HY),when F is 2-D,uses the spacing
specified by HX and HY,HX and HY can either be scalars to specify
the spacing between coordinates or vectors to specify the
coordinates of the points,If HX and HY are vectors,their length
must match the corresponding dimension of F.
六,插值,(Interpolation)
在已知数据之间计算估计值的过程。
1,一维插值(1D Interpolation)
由interp1实现，用多项式技术计算插值点。
Yi=interp1(x,y,xi,method) y—函数值矢量,x—自变量取值范围，xi—插值点的自变量矢量，
Method—插值方法选项。
MATLAB6.1的4种方法：
*临近点插值：method= ‘nearest’
*线性插值,method= ‘linear’
*三次样条插值：method= ‘spline’
*立方插值,method= ‘pchip’ or ‘cubic’
选择插值方法时主要考虑因素,运算时间、占用计算机内存和插值的光滑程度。比较：
运算时间,占用计算机内存 光滑程度。
*临近点插值,快 少 差
*线性插值,稍长 较多 稍好
*三次样条插值,最长 较多 最好
*立方插值,较长 多 较好
例1：一维插值函数插值方法的对比。
x=0:10;
y=sin(x);
xi=0:0.25:10;
strmod={'nearest','linear','spline','cubic'} %将插值方法定义为单元数组
str1b={'(a) method=nearest','(b) method=linear',...
'(c) method=spline','(d) method=cubic'} %将图标定义为单元数组
for i=1:4
yi=interp1(x,y,xi,strmod{i});
subplot(2,2,i)
plot(x,y,'ro',xi,yi,'b'),xlabel(str1b(i))
end
strmod =
'nearest' 'linear' 'spline' 'cubic'
例2,三次样条插值
x0=0:10;
y0=sin(x0);
x=0:.25:10;
y=spline(x0,y0,x);
plot(x0,y0,'or',x,y,'k')
与interp1结果一样
2,二维插值(2D Interpolation)
用于图形图象处理和三维曲线拟合等领域，由interp2实现，一般格式为：
ZI=interp2(X,Y,Z,XI,YI,method) X,Y—自变量组成的数组，尺寸相同。
XI,YI—插值点的自变量数组
Method—插值方法选项，4种
*临近点插值：method= ‘nearest’
*线性插值,method= ‘linear’ 该方法是interp2函数的缺省方法
*三次样条插值：method= ‘spline’
*立方插值,method= ‘pchip’ or ‘cubic’
例,二维插值4种方法的对比。
[x,y,z]=peaks(7); figure(1),mesh(x,y,z)
[xi,yi]=meshgrid(-3:0.2:3,-3:0.2:3);
z1=interp2(x,y,z,xi,yi,'nearest');
z2=interp2(x,y,z,xi,yi,'linear');
z3=interp2(x,y,z,xi,yi,'spline');
z4=interp2(x,y,z,xi,yi,'cubic');
figure(2),subplot(2,2,1)
mesh(xi,yi,z1)
title('nearest')
subplot(2,2,2)
mesh(xi,yi,z2)
title('linear')
subplot(2,2,3)
mesh(xi,yi,z3)
title('spiine')
subplot(2,2,4)
mesh(xi,yi,z4)
title('cubic')
3,多维插值,(3D Interpolation)
包括三维插值函数interp3和多维插值函数interpn,函数调用格式与一、二维插值基本相同。
VI=interp3(X,Y,Z,V,XI,YI,ZI,method)
其中,X,Y,Z—自变量组成的数组；
V—三维函数数组
XI,YI,ZI—插值点的自变量数组
Method --插值方法选项。
（FLOW A simple function of 3 variables.
FLOW,a function of three variables,is the speed profile of a
submerged jet within a infinite tank,FLOW is useful for
demonstrating SLICE and INTERP3.）
(SLICE(X,Y,Z,V,XI,YI,ZI) draws slices through the volume V along the
surface defined by the arrays XI,YI,ZI.)
(SHADING FLAT sets the shading of the current graph to flat.
SHADING INTERP sets the shading to interpolated.
SHADING FACETED sets the shading to faceted,which is the default.)
例：三维插值实例。
[x,y,z,v]=flow(10); %三变量无限大容器淹没射流场的速度剖面图
figure(1),
slice(x,y,z,v,[6 9.5],2,[-2,2]) %在X=6,9.5,Y=2,z=-2,0.2五处取切面
[xi,yi,zi]=meshgrid(.1:.25:10,-3:.25:3,-3:.25:3);
vi=interp3(x,y,z,v,xi,yi,zi); %全场沿坐标面插值
figure(2),
slice(xi,yi,zi,vi,[6 9.5],2,[-2,2]) %切面插值
shading flat
FLOW A simple function of 3 variables.
FLOW,a function of three variables,is the speed profile of a submerged jet within a infinite tank,FLOW is useful for demonstrating SLICE and INTERP3.
There are several variants of the calling sequence:
V = FLOW; produces a 50-by-25-by-25 array.
V = FLOW(N) produces a 2N-by-N-by-N array.
V = FLOW(X,Y,Z) evaluates the speed profile at the points (X,Y,Z).
[X,Y,Z,V] = FLOW(...) returns the coordinates as well.
七,数字信号处理初步(Introduction to Signal Processing )
MATLAB主要包括信号处理（Signal Processing Toolbox）和滤波器设计(Filter Design Toolbox)两部分.
基本平台中提供了一些常用的信号处理函数（表 4-4），可实现数字信号处理的基本功能。
1,快速傅立叶变换(FFT),MATLAB6.1 采用了新的快速傅立叶变换计算方法，速度高，可以作到实时处理。
函数fft的调用格式,
*Y=fft(X) 返回应用快速傅立叶方法计算得到的矢量X的离散傅立叶变换（DFT）,如果 X为矩阵，fft返回矩阵每一列的傅立叶变换，如果X为多维数组，fft运算从第一个非独立维开始执行。
*Y=fft(X,n) 返回n点的离散傅立叶变换，如果X的长度小于n，X中补0使其与n的长度相同，；如果X的长度大于n，则X的多出部分将被删除；如X为矩阵，用同样方法处理矩阵列的长度。
*Y=fft(X,[],dim) 和Y=fft(X,n,dim)沿dim维进行FFT操作。
注：快速傅立叶变换的结果为复数。
例,产生一个正弦衰减曲线，进行快速傅立叶变换，并画出幅值(amplitude)图，相位(phase)图、实部(real)图和虚部(image)图。
tp=0:2048; %时域数据点数N
yt=sin(0.08*pi*tp).*exp(-tp/80); %生成正弦衰减函数
figure(1),
plot(tp,yt),axis([0,400,-1,1]),%绘正弦衰减曲线
t=0:800/2048:800; %频域点数Nf
f=0:1.25:1000;
yf=fft(yt); %快速傅立叶变换
ya=abs(yf(1:801)); %幅值
yp=angle(yf(1:801))*180/pi; %相位
yr=real(yf(1:801)); %实部
yi=imag(yf(1:801)); %虚部
figure(2),
subplot(2,2,1)
plot(f,ya),axis([0,200,0,60]) %绘幅值曲线
subplot(2,2,2)
plot(f,yp),axis([0,200,-200,10]) %绘相位曲线
subplot(2,2,3)
plot(f,yr),axis([0,200,-40,40]) %绘实部曲线
subplot(2,2,4)
plot(f,yi),axis([0,200,-60,10]) %绘虚部曲线
FFT Discrete Fourier transform.
FFT(X) is the discrete Fourier transform (DFT) of vector X,For
matrices,the FFT operation is applied to each column,For N-D
arrays,the FFT operation operates on the first non-singleton
dimension.
FFT(X,N) is the N-point FFT,padded with zeros if X has less
than N points and truncated if it has more.
FFT(X,[],DIM) or FFT(X,N,DIM) applies the FFT operation across the dimension DIM.
For length N input vector x,the DFT is a length N vector X,with elements
N
X(k) = sum x(n)*exp(-j*2*pi*(k-1)*(n-1)/N),1 <= k <= N.
n=1
The inverse DFT (computed by IFFT) is given by
N
x(n) = (1/N) sum X(k)*exp( j*2*pi*(k-1)*(n-1)/N),1 <= n <= N.
k=1
2,快速傅立叶变换的长度与运算速度(Relation of the length of FFT and the velocity of calculation speed )
使用fft函数时，可输入第二个参数n以指定变换点的数量：
y=fft(x,n)
fft的运算速度取决于变换的长度。
*如n是2 的整数次幂，则运算速度最快；
*如n是合数，fft采用质因数分解的算法，速度取决于质因数的大小。
*如n是质数(prime number)，计算速度最慢。
例,创建70000×1阶的随机矢量x，取快速傅立叶变换的长度分别为质数65539,2的16次方,两个合数(composite number)66000和5535，分别计算使用这些长度fft所占用的cpu时间。
x=rand(70000,1);
isprime(65539) %质数
ans =
1
2^16
ans =
65536
（FACTOR(N) returns a vector containing the prime factors of N.）
factor(66000) %合数
ans =
2 2 2 2 3 5 5 5 11
factor(65535) %合数
ans =
3 5 17 257
t=cputime; %开始运行时间
y=fft(x,65539);
e=cputime-t %运行结束与开始运行时的时间差
e=
9900
t=cputime;
y=fft(x,65536);
e=cputime-t
e=
1600
t=cputime;
y=fft(x,65535);
e=cputime-t
e =
3800
由上例可见，变换长度的差别对运算速度的影响很大，这正是大多数的数字信号处理系统中FFT的长度取为2 的n次方的原因。
Chapter 4,Numerical computation of MATLAB
一、多项式(Polynomial)`
1．多项式的表达与创建(Expression and Creating of polynomial)
(1) 多项式的表达（expression of polynomial）_
Matlab用行矢量表达多项式系数(Coefficient)，各元素按变量的降幂顺序排列，如多项式为：
P(x)=a0xn+a1xn-1+a2xn-2…an-1x+an
则其系数矢量(Vector of coefficient)为：P=[a0 a1 … an-1 an]
如将根矢量(Vector of root)表示为：
ar=[ ar1 ar2 … arn]
则根矢量与系数矢量之间关系为：
（x-ar1）(x- ar2) … (x- arn)= a0xn+a1xn-1+a2xn-2…an-1x+an
(2)多项式的创建（polynomial creating）
a)系数矢量的直接输入法利用poly2sym函数直接输入多项式的系数矢量，就可方便的建立符号形式的多项式。
例：创建多项式x3-4x2+3x+2
poly2sym([1 -4 3 2])
ans =
x^3-4*x^2+3*x+2
POLY Convert roots to polynomial.
POLY(A),when A is an N by N matrix,is a row vector with
N+1 elements which are the coefficients of the characteristic polynomial,DET(lambda*EYE(SIZE(A)) - A),
POLY(V),when V is a vector,is a vector whose elements are
the coefficients of the polynomial whose roots are the
elements of V,For vectors,ROOTS and POLY are inverse
functions of each other,up to ordering,scaling,and
roundoff error.
b) 由根矢量创建多项式通过调用函数 p=poly(ar)产生多项式的系数矢量,再利用poly2sym函数就可方便的建立符号形式的多项式。
注：（1）根矢量元素为n,则多项式系数矢量元素为n+1；
（2）函数poly2sym(pa) 把多项式系数矢量表达成符号形式的多项式，缺省情况下自变量符号为x，可以指定自变量。
（3）使用简单绘图函数ezplot可以直接绘制符号形式多项式的曲线。
例 1：由根矢量创建多项式。将多项式(x-6)(x-3)(x-8)表示为系数形式
a=[6 3 8] %根矢量
pa=poly(a) %求系数矢量
ppa=poly2sym(pa) %以符号形式表示原多项式
ezplot(ppa,[-50,50])
pa =
1 -17 90 -144
ppa =
x^3-17*x^2+90*x-144
注：含复数根的根矢量所创建的多项式要注意：
（1）要形成实系数多项式，根矢量中的复数根必须共轭成对；
（2）含复数根的根矢量所创建的多项式系数矢量中，可能带有很小的虚部，此时可采用取实部的命令(real)把虚部滤掉。
进行多项式的求根运算时，有两种方法，一是直接调用求根函数roots，poly和 roots 互为逆函数。另一种是先把多项式转化为伴随矩阵，然后再求其特征值，该特征值即是多项式的根。
例 3,由给定复数根矢量求多项式系数矢量。
r=[-0.5 -0.3+0.4i -0.3-0.4i];
p=poly(r)
pr=real(p)
ppr=poly2sym(pr)
p =
1.0000 1.1000 0.5500 0.1250
pr =
1.0000 1.1000 0.5500 0.1250
ppr =
x^3+11/10*x^2+11/20*x+1/8
c) 特征多项式输入法用poly函数可实现由矩阵的特征多项式系数创建多项式。
条件：特征多项式系数矢量的第一个元素必须为一。
例 2,求三阶方阵A的特征多项式系数，并转换为多项式形式。
a=[6 3 8;7 5 6; 1 3 5]
Pa=poly(a) %求矩阵的特征多项式系数矢量
Ppa=poly2sym(pa)
Pa =
1.0000 -16.0000 38.0000 -83.0000
Ppa =
x^3-17*x^2+90*x-144
注：n 阶方阵的特征多项式系数矢量一定是n +1阶的。
注：（1）要形成实系数多项式，根矢量中的复数根必须共轭成对；
（2）含复数根的根矢量所创建的多项式系数矢量中，可能带有很小的虚部，此时可采用取实部的命令(real)把虚部滤掉。
进行多项式的求根运算时，有两种方法，一是直接调用求根函数roots，poly和 roots 互为逆函数。另一种是先把多项式转化为伴随矩阵，然后再求其特征值，该特征值即是多项式的根。
例 4,将多项式的系数表示形式转换为根表现形式。
求 x3-6x2-72x-27的根
a=[1 -6 -72 -27]
r=roots(a)
r =
12.1229
-5.7345
-0.3884
MATLAB约定，多项式系数矢量用行矢量表示，根矢量用列矢量表示。
>>
1,多项式的乘除运算(Multiplication and division of polynomial)
多项式乘法用函数conv(a,b)实现,除法用函数deconv(a,b)实现。
例1：a(s)=s2+2s+3,b(s)=4s2+5s+6,计算 a(s)与 b(s)的乘积。
a=[1 2 3]; b=[4 5 6];
c=conv(a,b)
cs=poly2sym(c,’s’)
c =
4 13 28 27 18
cs =
4*s^4+13*s^3+28*s^2+27*s+18
例2,展开(s2+2s+2)(s+4)(s+1) （多个多项式相乘）
c=conv([1,2,2],conv([1,4],[1,1]))
cs=poly2sym(c,’s’) %（指定变量为s）
c =
1 7 16 18 8
cs =
s^4+7*s^3+16*s^2+18*s+8
例2：求多项式s^4+7*s^3+16*s^2+18*s+8分别被(s+4),(s+3)除后的结果。
c=[1 7 16 18 8];
[q1,r1]=deconv(c,[1,4]) %q—商矢量,r—余数矢量
[q2,r2]=deconv(c,[1,3])
cc=conv(q2,[1,3]) %对除(s+3)结果检验
test=((c-r2)==cc)
q1 =
1 3 4 2
r1 =
0 0 0 0 0
q2 =
1 4 4 6
r2 =
0 0 0 0 -10
cc =
1 7 16 18 18
test =
1 1 1 1 1
1,其他常用的多项式运算命令(Other computation command of polynomial)
pa=polyval(p,s) 按数组运算规则计算给定s时多项式p的值。
pm=polyvalm(p,s) 按矩阵运算规则计算给定s时多项式p的值。
[r,p,k]=residue(b,a) 部分分式展开，b,a分别是分子分母多项式系数矢量，r,p,k分别是留数、极点和直项矢量
p=polyfit(x,y,n) 用n阶多项式拟合x，y矢量给定的数据。
polyder(p) 多项式微分。
注,对于多项式b(s)与不重根的n阶多项式a(s)之比，其部分分式展开为：
式中：p1,p2,…,pn称为极点(poles)，r1,r2,…,rn 称为留数(residues)，k(s)称为直项(direct terms)，假如a(s)含有m重根pj，则相应部分应写成：
RESIDUE Partial-fraction expansion (residues).
[R,P,K] = RESIDUE(B,A) finds the residues,poles and direct term of a partial fraction expansion of the ratio of two polynomials B(s)/A(s),If there are no multiple roots,
B(s) R(1) R(2) R(n)
---- = -------- + -------- +,.,+ -------- + K(s)
A(s) s - P(1) s - P(2) s - P(n)
Vectors B and A specify the coefficients of the numerator and denominator polynomials in descending powers of s,The residues
are returned in the column vector R,the pole locations in column vector P,and the direct terms in row vector K,The number of poles is n = length(A)-1 = length(R) = length(P),The direct term coefficient vector is empty if length(B) < length(A),otherwise
length(K) = length(B)-length(A)+1.
If P(j) =,.,= P(j+m-1) is a pole of multplicity m,then the expansion includes terms of the form
R(j) R(j+1) R(j+m-1)
-------- + ------------ +,.,+ ------------
s - P(j) (s - P(j))^2 (s - P(j))^m
[B,A] = RESIDUE(R,P,K),with 3 input arguments and 2 output arguments,converts the partial fraction expansion back to the polynomials with coefficients in B and A.
例3：对 (3x4+2x3+5x2+4x+6)/(x5+3x4+4x3+2x2+7x+2) 做部分分式展开
a=[1 3 4 2 7 2];
b=[3 2 5 4 6];
[r,s,k]=residue(b,a)
r =
1.1274 + 1.1513i
1.1274 - 1.1513i
-0.0232 - 0.0722i
-0.0232 + 0.0722i
0.7916
s =
-1.7680 + 1.2673i
-1.7680 - 1.2673i
0.4176 + 1.1130i
0.4176 - 1.1130i
-0.2991
k =
[] (分母阶数高于分子阶数时，k将是空矩阵，表示无此项)
例5：对一组实验数据进行多项式最小二乘拟合(least square fit)
x=[1 2 3 4 5]; % 实验数据
y=[5.5 43.1 128 290.7 498.4];
p=polyfit(x,y,3) %做三阶多项式拟合
x2=1:.1:5;
y2=polyval(p,x2); % 根据给定值计算多项式结果
plot(x,y,’o’,x2,y2)
线性代数(Linear Algebra)
解线性方程(Linear equation)就是找出是否存在一个唯一的矩阵x，使得a,b满足关系：
ax=b 或 xa=b
MALAB中x=a\b 是方程 ax=b 的解,x=b/a是方程式xa=b的解。
通常线性方程多写成ax=b，“\”较多用，两者的关系为：
(b/a)’=(a’\b’)
系数矩阵a可能是m行n列的，有三种情况：
*方阵系统,(Square matrix) m=n 可求出精确解（a必须是非奇异(nonsingular)，即满秩(full rank)）
*超定系统：(Overdetermind system)m>n 可求出最小二乘解
*欠定系统：(Underdetermind system) m<n 可尝试找出含有最少m个基解或最小范数解
MATLAB对不同形式的参数矩阵，采用不同的运算法则来处理，它会自动检测参数矩阵，以区别下面几种形式：
*三角矩阵(Triangular Matrix)
*对称正定矩阵(symmetrical positive determined matrix)
*非奇异方阵(Nonsingular matrix)
*超定系统(Overdetermind system)
*欠定系统(Underdetermind system)
方阵系统：(Square array)
最常见的是系数矩阵为方阵a，常数项b为列矢量,其解x可写成x=a\b,x和b大小相同。
例1,求方阵系统的根。
a=[11 6 7; 5 13 9; 17 1 8]
b=[16 13 4]’
x=a\b
a =
11 6 7
5 13 9
17 1 8
b =
16
13
4
x =
3.9763
5.4455
-8.6303
例2：假如a,b为两个大小相同的矩阵，求方阵系统的根。
a=[4 5 9; 18 19 5; 1 4 13]
b=[1 5 12; 3 15 19; 7 6 10]
x=a\b
C=a*x
a =
4 5 9
18 19 5
1 4 13
b =
1 5 12
3 15 19
7 6 10
x =
-3.6750 -0.7333 2.9708
3.7250 1.4667 -2.1292
-0.3250 0.0667 1.1958
C =
1.0000 5.0000 12.0000
3.0000 15.0000 19.0000
7.0000 6.0000 10.0000
若方阵a的各个行矢量线性相关(linear correlation)，则称方阵a为奇异矩阵。这时线性方程将有无穷多组解。若方阵是奇异矩阵，则反斜线运算因子将发出警告信息。
2．超定系统(Overdetermind system)
实验数据较多，寻求他们的曲线拟合。
如在t内测得一组数据y：
t y
0.0 0.82
0.3 0.72
0.8 0.63
0.60
1.6 0.55
2.2 0.50
这些数据显然有衰减指数趋势,y(t)~c1+c2e-t
此方程意为y矢量可以由两个矢量逐步逼近而得，一个是单行的常数矢量，一个是由指数e-t项构成，两个参数c1和c2可用最小二乘法求得，它们表示实验数据与方程y(t)~c1+c2e-t之间距离的最小平方和。
例1,求上述数据的最小二乘解。将数据带入方程式y(t)~c1+c2e-t中，可得到含有两个未知数的6个等式，可写成6行2 列的矩阵e.
t=[0 0.3 0.8 1.1 1.6 2.2]’;
y=[0.82 0.72 0.63 0.60 0.55 0.50]’;
e=[ones(size(t)) exp(-t)] %求6个y(t)方程的系数矩阵
c=e\y % 求方程的解
e =
1.0000 1.0000
1.0000 0.7408
1.0000 0.4493
1.0000 0.3329
1.0000 0.2019
1.0000 0.1108
c =
0.4744
0.3434
带入方程得：y(t)~0.4744+0.3434e-t
用此方程可绘制曲线：
t=[0 0.3 0.8 1.1 1.6 2.2]’;
y=[0.82 0.72 0.63 0.60 0.55 0.50]’;
t1=[0:0.1:2.5]’; y1=[ones(size(t1)),exp(-t1)]*c
plot(t1,y1,’b’,t,y,’ro’)
如果一个矩阵的行矢量是线性相关的，则它的最小二乘解并不唯一，因此，a\b运算将给出警告，并产生含有最少元素的基解。
3,欠定系统,(Underdetermind system)
欠定系统为线性相关系统，其解都不唯一，MATLAB会计算一组构成通解的基解，而方程的特解则用QR分解法决定。
两种解法：最少元素解a\b，最小范数解pinv(a)*b.
例,用两种方法求解欠定系统。
对a和矢量b分别用a\b和pinv(a)*b求解:
a=[1 1 1; 1 1 -1]
b=[10 6]’
p=a\b
q=pinv(a)*b
a =
1 1 1
1 1 -1
b =
10
6
p =
8.0000
0
2.0000
q =
4.0000
4.0000
2.0000
三,逆矩阵及行列式(Revers and determinant of matrix)
1,方阵的逆和行列式(Revers and determinant of square matrix)
若a是方阵，且为非奇异阵，则方程ax=I和 xa=I有相同的解X。X称为a的逆矩阵，记做a-1，在MATLAB中 用inv 函数来计算矩阵的逆。计算方阵的行列式则用det函数。
DET Determinant.
DET(X) is the determinant of the square matrix X.
Use COND instead of DET to test for matrix singularity.
INV Matrix inverse.
INV(X) is the inverse of the square matrix X,A warning message is printed if X is badly scaled or nearly singular.
例：计算方阵的行列式和逆矩阵。
a=[3 -3 1;-3 5 -2;1 -2 1];
b=[14 13 5; 5 1 12;6 14 5];
d1=det(a)
x1=inv(a)
d2=det(b)
x2=inv(b)
d1 =
1
x1 =
1.0000 1.0000 1.0000
1.0000 2.0000 3.0000
1.0000 3.0000 6.0000
d2 =
-1351
x2 =
0.1207 -0.0037 -0.1118
-0.0348 -0.0296 0.1058
-0.0474 0.0873 0.0377
2,广义逆矩阵（伪逆）(Generalized inverse matrix)
一般非方阵无逆矩阵和行列式，方程ax=I 和xa=I至少有一个无解，这种矩阵可以求得特殊的逆矩阵，成为广义逆矩阵(generalized inverse matrix)（或伪逆 pseudoinverse）。矩阵amn存在广义逆矩阵xnm，使得 ax=Imn,MATLAB用pinv函数来计算广义逆矩阵。
例：计算广义逆矩阵。
a=[8 14; 1 3; 9 6]
x=pinv(a)
b=x*a
c=a*x
d=c*a %d=a*x*a=a
e=x*c %e=x*a*x=x
a =
8 14
1 3
9 6
x =
-0.0661 -0.0402 0.1743
0.1045 0.0406 -0.0974
b =
1.0000 -0.0000
-0.0000 1.0000
c =
0.9334 0.2472 0.0317
0.2472 0.0817 -0.1177
0.0317 -0.1177 0.9849
d =
8.0000 14.0000
1.0000 3.0000
9.0000 6.0000
e =
-0.0661 -0.0402 0.1743
0.1045 0.0406 -0.0974
PINV Pseudoinverse.
X = PINV(A) produces a matrix X of the same dimensions as A' so that A*X*A = A,X*A*X = X and A*X and X*A are Hermitian,The computation is based on SVD(A) and any singular values less than a tolerance are treated as zero.
The default tolerance is MAX(SIZE(A)) * NORM(A) * EPS.
PINV(A,TOL) uses the tolerance TOL instead of the default.
四,矩阵分解(Matrix decomposition)
MATLAB求解线性方程的过程基于三种分解法则：
（１）Cholesky分解，针对对称正定矩阵；
（２）高斯消元法,针对一般矩阵；
（３）正交化,针对一般矩阵（行数≠列数）
这三种分解运算分别由chol,lu和 qr三个函数来分解.
Cholesky分解(Cholesky Decomposition)
仅适用于对称和上三角矩阵例：cholesky分解。
a=pascal(6)
b=chol(a)
a =
1 1 1 1 1 1
1 2 3 4 5 6
1 3 6 10 15 21
1 4 10 20 35 56
1 5 15 35 70 126
1 6 21 56 126 252
b =
1 1 1 1 1 1
0 1 2 3 4 5
0 0 1 3 6 10
0 0 0 1 4 10
0 0 0 0 1 5
0 0 0 0 0 1
CHOL Cholesky factorization.
CHOL(X) uses only the diagonal and upper triangle of X,The lower triangular is assumed to be the (complex conjugate) transpose of the upper,If X is positive definite,then R = CHOL(X) produces an upper triangular R so that R'*R = X,If X is not positive definite,an error message is printed.
[R,p] = CHOL(X),with two output arguments,never produces an
error message,If X is positive definite,then p is 0 and R is the same as above,But if X is not positive definite,then p is a positive integer.
When X is full,R is an upper triangular matrix of order q = p-1
so that R'*R = X(1:q,1:q),When X is sparse,R is an upper triangular matrix of size q-by-n so that the L-shaped region of the first q rows and first q columns of R'*R agree with those of X.
2,LU分解(LU factorization).
用lu函数完成LU分解，将矩阵分解为上、下两个三角阵，其调用格式为：
[l,u]=lu(a) l代表下三角阵，u代表上三角阵。
例：
LU分解。
a=[47 24 22; 11 44 0;30 38 41]
[l,u]=lu(a)
a =
47 24 22
11 44 0
30 38 41
l =
1.0000 0 0
0.2340 1.0000 0
0.6383 0.5909 1.0000
u =
47.0000 24.0000 22.0000
0 38.3830 -5.1489
0 0 30.0000
LU LU factorization.
[L,U] = LU(X) stores an upper triangular matrix in U and a "psychologically lower triangular matrix" (i.e,a product of lower triangular and permutation matrices) in L,so that X = L*U,X can be rectangular,
[L,U,P] = LU(X) returns unit lower triangular matrix L,upper triangular matrix U,and permutation matrix P so that P*X = L*U.
3,QR分解(Orthogonal-triangular decomposition).
函数调用格式：[q,r]=qr(a),q代表正规正交矩阵，r代表三角形矩阵。原始阵a不必一定是方阵。如果矩阵a是m×n阶的，则矩阵q是m×m阶的，矩阵r是m×n阶的。
例：QR分解.
A=[22 46 20 20; 30 36 46 44;39 8 45 2];
[q,r]=qr(A)
q =
-0.4082 -0.7209 -0.5601
-0.5566 -0.2898 0.7786
-0.7236 0.6296 -0.2829
r =
-53.8981 -44.6027 -66.3289 -34.1014
0 -38.5564 0.5823 -25.9097
0 0 11.8800 22.4896
QR Orthogonal-triangular decomposition.
[Q,R] = QR(A) produces an upper triangular matrix R of the same
dimension as A and a unitary matrix Q so that A = Q*R.
[Q,R,E] = QR(A) produces a permutation matrix E,an upper
triangular R and a unitary Q so that A*E = Q*R,The column
permutation E is chosen so that abs(diag(R)) is decreasing,
[Q,R] = QR(A,0) produces the "economy size" decomposition,If A is m-by-n with m > n,then only the first n columns of Q are computed.
4,特征值与特征矢量(Eigenvalues and eigenvectors).
MATLAB中使用函数eig计算特征值和 特征矢量，有两种调用方法：
*e=eig(a),其中e是包含特征值的矢量；
*[v,d]=eig(a),其中v是一个与a相同的n×n阶矩阵，它的每一列是矩阵a的一个特征值所对应的特征矢量，d为对角阵，其对角元素即为矩阵a的特征值。
例：计算特征值和特征矢量。
a=[34 25 15; 18 35 9; 41 21 9]
e=eig(a)
[v,d]=eig(a)
a =
34 25 15
18 35 9
41 21 9
e =
68.5066
15.5122
-6.0187
v =
-0.6227 -0.4409 -0.3105
-0.4969 0.6786 -0.0717
-0.6044 -0.5875 0.9479
d =
68.5066 0 0
0 15.5122 0
0 0 -6.0187
EIG Eigenvalues and eigenvectors.
E = EIG(X) is a vector containing the eigenvalues of a square matrix X.
[V,D] = EIG(X) produces a diagonal matrix D of eigenvalues and a full matrix V whose columns are the corresponding eigenvectors so that X*V = V*D.
[V,D] = EIG(X,'nobalance') performs the computation with balancing
disabled,which sometimes gives more accurate results for certain
problems with unusual scaling,If X is symmetric,EIG(X,'nobalance')
is ignored since X is already balanced.
5,奇异值分解.( Singular value decomposition).
如存在两个矢量u,v及一常数c,使得矩阵A满足：Av=cu,A’u=cv
称c为奇异值，称u,v为奇异矢量。
将奇异值写成对角方阵∑，而相对应的奇异矢量作为列矢量则可写成两个正交矩阵U，V,使得,AV=U∑,A‘U=V∑ 因为U，V正交，所以可得奇异值表达式：
A=U∑V’。
一个m行n列的矩阵A经奇异值分解，可求得m行m列的U,m行n列的矩阵∑和n行n列的矩阵V.。
奇异值分解用svd函数实现，调用格式为；
[u,s,v]=svd(a)
SVD Singular value decomposition.
[U,S,V] = SVD(X) produces a diagonal matrix S,of the same dimension as X and with nonnegative diagonal elements in decreasing order,and unitary matrices U and V so that X = U*S*V'.
S = SVD(X) returns a vector containing the singular values.
[U,S,V] = SVD(X,0) produces the "economy size" decomposition,If X is m-by-n with m > n,then only the first n columns of U are computed and S is n-by-n.
例,奇异值分解。
a=[8 5; 7 3;4 6];
[u,s,v]=svd(a) % s为奇异值对角方阵
u =
-0.6841 -0.1826 -0.7061
-0.5407 -0.5228 0.6591
-0.4895 0.8327 0.2589
s =
13.7649 0
0 3.0865
0 0
v =
-0.8148 -0.5797
-0.5797 0.8148
五,数据分析(Data Analyaia)
MATLAB对数据分析有两条约定：
（1） 若输入量X是矢量，则不论是行矢量还是列矢量，运算是对整个矢量进行的；
（2）若输入量X是数组，（或称矩阵），则命令运算是按列进行的。即默认每个列是有一个变量的不同“观察“所得的数据组成。
1,基本统计命令 （表4-1）
例,做各种基本统计运算。
A=[5 -10 -6 0;2 6 3 -3;-9 5 -10 11;-22 17 0 -19;-1 6 -4 4]
Amax=max(A) %找A各列的最大元素
Amin=min(A) %找A各列的最小元素
Amed=median(A) %找A各列的中位元素
Amean=mean(A) %找A各列的平均值
Astd=std(A) %求A各列的标准差
Aprod=prod(A) %求A各列元素的积
Asum=sum(A) %求A各列元素的和
S=cumsum(A) %求A各列元素的累积和
P=cumprod(A) %求A各列元素的累积j积
I=sort(A) %使A的各列元素按递增排列
A =
5 -10 -6 0
2 6 3 -3
-9 5 -10 11
-22 17 0 -19
-1 6 -4 4
Amax =
5 17 3 11
Amin =
-22 -10 -10 -19
Amed =
-1 6 -4 0
Amean =
-5.0000 4.8000 -3.4000 -1.4000
Astd =
10.8397 9.6281 5.0794 11.1490
Aprod =
-1980 -30600 0 0
Asum =
-25 24 -17 -7
S =
5 -10 -6 0
7 -4 -3 -3
-2 1 -13 8
-24 18 -13 -11
-25 24 -17 -7
P =
5 -10 -6 0
10 -60 -18 0
-90 -300 180 0
1980 -5100 0 0
-1980 -30600 0 0
I =
-22 -10 -10 -19
-9 5 -6 -3
-1 6 -4 0
2 6 0 4
5 17 3 11
>>
求矩阵元素的最大值、最小值可用：
Amax=max(maxA)) 或 Amax=max(A(:))，
Amin=min(min(A)) 或 Amin=min(A(:))
2．协方差阵和相关阵(Covariance matrix and Correlation coefficients).
（表 4—2）
例,计算协方差和相关阵。
x=rand(10,3);
y=rand(10,3);
cx=cov(x) %求协方差阵
cy=cov(y)
cxy=cov(x,y) %求两随机变量的协方差
px=corrcoef(x) %求相关阵
pxy=corrcoef(x,y) %求两随机变量的（2×2）相关系数
cx =
0.0483 -0.0066 0.0146
-0.0066 0.0283 0.0154
0.0146 0.0154 0.0978
cy =
0.1177 0.0073 -0.0127
0.0073 0.0239 -0.0230
-0.0127 -0.0230 0.0772
cxy =
0.0550 0.0023
0.0023 0.0697
px =
1.0000 -0.1783 0.2118
-0.1783 1.0000 0.2934
0.2118 0.2934 1.0000
pxy =
1.0000 0.0372
0.0372 1.0000
COV Covariance matrix.
COV(X),if X is a vector,returns the variance,For matrices,where each row is an observation,and each column a variable,COV(X) is the covariance matrix,DIAG(COV(X)) is a vector of variances for each column,and SQRT(DIAG(COV(X))) is a vector of standard deviations,
COV(X,Y),where X and Y are vectors of equal length,is equivalent to COV([X(:) Y(:)]),
COV(X) or COV(X,Y) normalizes by (N-1) where N is the number of observations,This makes COV(X) the best unbiased estimate of the
covariance matrix if the observations are from a normal distribution.
CORRCOEF Correlation coefficients.
R=CORRCOEF(X) calculates a matrix R of correlation coefficients for an array X,in which each row is an observation and each column is a variable.
R=CORRCOEF(X,Y),where X and Y are column vectors,is the same as
R=CORRCOEF([X Y]),If C is the covariance matrix,C = COV(X),then CORRCOEF(X) is the matrix whose (i,j)'th element is
C(i,j)/SQRT(C(i,i)*C(j,j)).
2,微分与梯度 (Difference and approximate derivative,gradient).
(表4—3)
例1：按列求微分。
x=[1,10,20;2,12,23;3,14,26;3,16,29]
d=diff(x) %求一阶微分
x =
1 10 20
2 12 23
3 14 26
3 16 29
d =
1 2 3
1 2 3
0 2 3
例2,对于（u=x2+y2和Δ2=4）求5点差分。
[x,y]=meshgrid(-4:4,-3:3);
u=x.^2+y.^2
v4=4*del2(u) %求m×n阶矩阵U的五点差分矩阵
u =
25 18 13 10 9 10 13 18 25
20 13 8 5 4 5 8 13 20
17 10 5 2 1 2 5 10 17
16 9 4 1 0 1 4 9 16
17 10 5 2 1 2 5 10 17
20 13 8 5 4 5 8 13 20
25 18 13 10 9 10 13 18 25
v4 =
4 4 4 4 4 4 4 4 4
4 4 4 4 4 4 4 4 4
4 4 4 4 4 4 4 4 4
4 4 4 4 4 4 4 4 4
4 4 4 4 4 4 4 4 4
4 4 4 4 4 4 4 4 4
4 4 4 4 4 4 4 4 4
MESHGRID X and Y arrays for 3-D plots.
[X,Y] = MESHGRID(x,y) transforms the domain specified by vectors
x and y into arrays X and Y that can be used for the evaluation
of functions of two variables and 3-D surface plots.
The rows of the output array X are copies of the vector x and
the columns of the output array Y are copies of the vector y.
[X,Y] = MESHGRID(x) is an abbreviation for [X,Y] = MESHGRID(x,x).
[X,Y,Z] = MESHGRID(x,y,z) produces 3-D arrays that can be used to
evaluate functions of three variables and 3-D volumetric plots.
DEL2 Discrete Laplacian.
L = DEL2(U) when U is a matrix,is an discrete approximation of 0.25*del^2 u = (d^2u/dx^2 + d^2/dy^2)/4,The matrix L is the same size as U with each element equal to the difference between an element of U and the average of its four neighbors.
L = DEL2(U) when U is an N-D array,returns an approximation of
(del^2 u)/2/n where n is ndims(u).
L = DEL2(U,H),where H is a scalar,uses H as the spacing between
points in each direction (H=1 by default).
L = DEL2(U,HX,HY) when U is 2-D,uses the spacing specified by HX
and HY,If HX is a scalar,it gives the spacing between points in
the x-direction,If HX is a vector,it must be of length SIZE(U,2)
and specifies the x-coordinates of the points,Similarly,if HY
is a scalar,it gives the spacing between points in the
y-direction,If HY is a vector,it must be of length SIZE(U,1) and
specifies the y-coordinates of the points.
L = DEL2(U,HX,HY,HZ,...) when U is N-D,uses the spacing given by
HX,HY,HZ,etc.
例3：产生一个二元函数偏导数和梯度。
x=-2:0.2:2;
y=-2:0.2:2;
[xx,yy]=meshgrid(x,y);
z=xx.*exp(-xx.^2-yy.^2);
[Gx,Gy]=gradient(z,0.2,0.2); %Gx,Gy分别是二元函数的偏导
contour(x,y,z,'k'),hold on,
quiver(xx,yy,Gx,Gy,'r'),hold off
DIFF Difference and approximate derivative.
DIFF(X),for a vector X,is [X(2)-X(1) X(3)-X(2),.,X(n)-X(n-1)].
DIFF(X),for a matrix X,is the matrix of row differences,
[X(2:n,:) - X(1:n-1,:)].
DIFF(X),for an N-D array X,is the difference along the first
non-singleton dimension of X.
DIFF(X,N) is the N-th order difference along the first non-singleton dimension (denote it by DIM),If N >= size(X,DIM),DIFF takes successive differences along the next non-singleton dimension.
（[FX,FY] = GRADIENT(F,HX,HY),when F is 2-D,uses the spacing
specified by HX and HY,HX and HY can either be scalars to specify
the spacing between coordinates or vectors to specify the coordinates of the points,If HX and HY are vectors,their length must match the corresponding dimension of F.）
（QUIVER(X,Y,U,V) plots velocity vectors as arrows with components (u,v)
at the points (x,y),The matrices X,Y,U,V must all be the same size
and contain corresponding position and velocity components (X and Y
can also be vectors to specify a uniform grid),QUIVER automatically
scales the arrows to fit within the grid.）
GRADIENT Approximate gradient.
[FX,FY] = GRADIENT(F) returns the numerical gradient of the
matrix F,FX corresponds to dF/dx,the differences in the
x (column) direction,FY corresponds to dF/dy,the differences
in the y (row) direction,The spacing between points in each
direction is assumed to be one,When F is a vector,DF = GRADIENT(F) is the 1-D gradient.
[FX,FY] = GRADIENT(F,H),where H is a scalar,uses H as the
spacing between points in each direction.
[FX,FY] = GRADIENT(F,HX,HY),when F is 2-D,uses the spacing
specified by HX and HY,HX and HY can either be scalars to specify
the spacing between coordinates or vectors to specify the
coordinates of the points,If HX and HY are vectors,their length
must match the corresponding dimension of F.
六,插值,(Interpolation)
在已知数据之间计算估计值的过程。
1,一维插值(1D Interpolation)
由interp1实现，用多项式技术计算插值点。
Yi=interp1(x,y,xi,method) y—函数值矢量,x—自变量取值范围，xi—插值点的自变量矢量，
Method—插值方法选项。
MATLAB6.1的4种方法：
*临近点插值：method= ‘nearest’
*线性插值,method= ‘linear’
*三次样条插值：method= ‘spline’
*立方插值,method= ‘pchip’ or ‘cubic’
选择插值方法时主要考虑因素,运算时间、占用计算机内存和插值的光滑程度。比较：
运算时间,占用计算机内存 光滑程度。
*临近点插值,快 少 差
*线性插值,稍长 较多 稍好
*三次样条插值,最长 较多 最好
*立方插值,较长 多 较好
例1：一维插值函数插值方法的对比。
x=0:10;
y=sin(x);
xi=0:0.25:10;
strmod={'nearest','linear','spline','cubic'} %将插值方法定义为单元数组
str1b={'(a) method=nearest','(b) method=linear',...
'(c) method=spline','(d) method=cubic'} %将图标定义为单元数组
for i=1:4
yi=interp1(x,y,xi,strmod{i});
subplot(2,2,i)
plot(x,y,'ro',xi,yi,'b'),xlabel(str1b(i))
end
strmod =
'nearest' 'linear' 'spline' 'cubic'
例2,三次样条插值
x0=0:10;
y0=sin(x0);
x=0:.25:10;
y=spline(x0,y0,x);
plot(x0,y0,'or',x,y,'k')
与interp1结果一样
2,二维插值(2D Interpolation)
用于图形图象处理和三维曲线拟合等领域，由interp2实现，一般格式为：
ZI=interp2(X,Y,Z,XI,YI,method) X,Y—自变量组成的数组，尺寸相同。
XI,YI—插值点的自变量数组
Method—插值方法选项，4种
*临近点插值：method= ‘nearest’
*线性插值,method= ‘linear’ 该方法是interp2函数的缺省方法
*三次样条插值：method= ‘spline’
*立方插值,method= ‘pchip’ or ‘cubic’
例,二维插值4种方法的对比。
[x,y,z]=peaks(7); figure(1),mesh(x,y,z)
[xi,yi]=meshgrid(-3:0.2:3,-3:0.2:3);
z1=interp2(x,y,z,xi,yi,'nearest');
z2=interp2(x,y,z,xi,yi,'linear');
z3=interp2(x,y,z,xi,yi,'spline');
z4=interp2(x,y,z,xi,yi,'cubic');
figure(2),subplot(2,2,1)
mesh(xi,yi,z1)
title('nearest')
subplot(2,2,2)
mesh(xi,yi,z2)
title('linear')
subplot(2,2,3)
mesh(xi,yi,z3)
title('spiine')
subplot(2,2,4)
mesh(xi,yi,z4)
title('cubic')
3,多维插值,(3D Interpolation)
包括三维插值函数interp3和多维插值函数interpn,函数调用格式与一、二维插值基本相同。
VI=interp3(X,Y,Z,V,XI,YI,ZI,method)
其中,X,Y,Z—自变量组成的数组；
V—三维函数数组
XI,YI,ZI—插值点的自变量数组
Method --插值方法选项。
（FLOW A simple function of 3 variables.
FLOW,a function of three variables,is the speed profile of a
submerged jet within a infinite tank,FLOW is useful for
demonstrating SLICE and INTERP3.）
(SLICE(X,Y,Z,V,XI,YI,ZI) draws slices through the volume V along the
surface defined by the arrays XI,YI,ZI.)
(SHADING FLAT sets the shading of the current graph to flat.
SHADING INTERP sets the shading to interpolated.
SHADING FACETED sets the shading to faceted,which is the default.)
例：三维插值实例。
[x,y,z,v]=flow(10); %三变量无限大容器淹没射流场的速度剖面图
figure(1),
slice(x,y,z,v,[6 9.5],2,[-2,2]) %在X=6,9.5,Y=2,z=-2,0.2五处取切面
[xi,yi,zi]=meshgrid(.1:.25:10,-3:.25:3,-3:.25:3);
vi=interp3(x,y,z,v,xi,yi,zi); %全场沿坐标面插值
figure(2),
slice(xi,yi,zi,vi,[6 9.5],2,[-2,2]) %切面插值
shading flat
FLOW A simple function of 3 variables.
FLOW,a function of three variables,is the speed profile of a submerged jet within a infinite tank,FLOW is useful for demonstrating SLICE and INTERP3.
There are several variants of the calling sequence:
V = FLOW; produces a 50-by-25-by-25 array.
V = FLOW(N) produces a 2N-by-N-by-N array.
V = FLOW(X,Y,Z) evaluates the speed profile at the points (X,Y,Z).
[X,Y,Z,V] = FLOW(...) returns the coordinates as well.
七,数字信号处理初步(Introduction to Signal Processing )
MATLAB主要包括信号处理（Signal Processing Toolbox）和滤波器设计(Filter Design Toolbox)两部分.
基本平台中提供了一些常用的信号处理函数（表 4-4），可实现数字信号处理的基本功能。
1,快速傅立叶变换(FFT),MATLAB6.1 采用了新的快速傅立叶变换计算方法，速度高，可以作到实时处理。
函数fft的调用格式,
*Y=fft(X) 返回应用快速傅立叶方法计算得到的矢量X的离散傅立叶变换（DFT）,如果 X为矩阵，fft返回矩阵每一列的傅立叶变换，如果X为多维数组，fft运算从第一个非独立维开始执行。
*Y=fft(X,n) 返回n点的离散傅立叶变换，如果X的长度小于n，X中补0使其与n的长度相同，；如果X的长度大于n，则X的多出部分将被删除；如X为矩阵，用同样方法处理矩阵列的长度。
*Y=fft(X,[],dim) 和Y=fft(X,n,dim)沿dim维进行FFT操作。
注：快速傅立叶变换的结果为复数。
例,产生一个正弦衰减曲线，进行快速傅立叶变换，并画出幅值(amplitude)图，相位(phase)图、实部(real)图和虚部(image)图。
tp=0:2048; %时域数据点数N
yt=sin(0.08*pi*tp).*exp(-tp/80); %生成正弦衰减函数
figure(1),
plot(tp,yt),axis([0,400,-1,1]),%绘正弦衰减曲线
t=0:800/2048:800; %频域点数Nf
f=0:1.25:1000;
yf=fft(yt); %快速傅立叶变换
ya=abs(yf(1:801)); %幅值
yp=angle(yf(1:801))*180/pi; %相位
yr=real(yf(1:801)); %实部
yi=imag(yf(1:801)); %虚部
figure(2),
subplot(2,2,1)
plot(f,ya),axis([0,200,0,60]) %绘幅值曲线
subplot(2,2,2)
plot(f,yp),axis([0,200,-200,10]) %绘相位曲线
subplot(2,2,3)
plot(f,yr),axis([0,200,-40,40]) %绘实部曲线
subplot(2,2,4)
plot(f,yi),axis([0,200,-60,10]) %绘虚部曲线
FFT Discrete Fourier transform.
FFT(X) is the discrete Fourier transform (DFT) of vector X,For
matrices,the FFT operation is applied to each column,For N-D
arrays,the FFT operation operates on the first non-singleton
dimension.
FFT(X,N) is the N-point FFT,padded with zeros if X has less
than N points and truncated if it has more.
FFT(X,[],DIM) or FFT(X,N,DIM) applies the FFT operation across the dimension DIM.
For length N input vector x,the DFT is a length N vector X,with elements
N
X(k) = sum x(n)*exp(-j*2*pi*(k-1)*(n-1)/N),1 <= k <= N.
n=1
The inverse DFT (computed by IFFT) is given by
N
x(n) = (1/N) sum X(k)*exp( j*2*pi*(k-1)*(n-1)/N),1 <= n <= N.
k=1
2,快速傅立叶变换的长度与运算速度(Relation of the length of FFT and the velocity of calculation speed )
使用fft函数时，可输入第二个参数n以指定变换点的数量：
y=fft(x,n)
fft的运算速度取决于变换的长度。
*如n是2 的整数次幂，则运算速度最快；
*如n是合数，fft采用质因数分解的算法，速度取决于质因数的大小。
*如n是质数(prime number)，计算速度最慢。
例,创建70000×1阶的随机矢量x，取快速傅立叶变换的长度分别为质数65539,2的16次方,两个合数(composite number)66000和5535，分别计算使用这些长度fft所占用的cpu时间。
x=rand(70000,1);
isprime(65539) %质数
ans =
1
2^16
ans =
65536
（FACTOR(N) returns a vector containing the prime factors of N.）
factor(66000) %合数
ans =
2 2 2 2 3 5 5 5 11
factor(65535) %合数
ans =
3 5 17 257
t=cputime; %开始运行时间
y=fft(x,65539);
e=cputime-t %运行结束与开始运行时的时间差
e=
9900
t=cputime;
y=fft(x,65536);
e=cputime-t
e=
1600
t=cputime;
y=fft(x,65535);
e=cputime-t
e =
3800
由上例可见，变换长度的差别对运算速度的影响很大，这正是大多数的数字信号处理系统中FFT的长度取为2 的n次方的原因。
