第四章 MATLAB的数值计算功能(2)

更新时间：2024-06-09 00:46:01 阅读量：综合文库文档下载

说明：文章内容仅供预览，部分内容可能不全。下载后的文档，内容与下面显示的完全一致。下载之前请确认下面内容是否您想要的，是否完整无缺。

第四章净身出户离婚的小说推荐度：
相关推荐

第四章 MATLAB 的数值计算功能

Chapter 4: Numerical computation of MATLAB

数值计算是MATLAB最基本、最重要的功能，是MATLAB最具代表性的特点。MATLAB在数值计算过程中以数组和矩阵为基础。数组是MATLAB运算中的重要数据组织形式。前面章节对数组、矩阵的特征及其创建与基本运算规则等相关知识已作了较详尽的介绍，本章重点介绍常用的数值计算方法。

一、多项式(Polynomial)`

多项式在众多学科的计算中具有重要的作用，许多方程和定理都是多项式的形式。MATLAB提供了标准多项式运算的函数，如多项式的求根、求值和微分，还提供了一些用于更高级运算的函数，如曲线拟合和多项式展开等。

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函数直接输入多项式的系数矢量，就可方便的建立符号形式的多项式。

例1：创建给定的多项式x3-4x2+3x+2

poly2sym([1 -4 3 2]) ans =

x^3-4*x^2+3*x+2

也可以用poly2str.求一个方阵对应的符号形式的多项式。例 2：

a=[3 4 6;5 7 1;8 3 9]

p=poly(a); %求方阵的特征多项式系数矢量 poly2str(p,'x') %建立符号形式的多项式 ans =

x^3 - 19 x^2 + 40 x + 214

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.

POLY2SYM Polynomial coefficient vector to symbolic polynomial.

POLY2SYM(C) is a symbolic polynomial in 'x' with coefficients from the vector C.POLY2SYM(C,'V') and POLY2SYM(C,SYM('V') both use the symbolic variable specified by the second argument. POLY2STR Return polynomial as string.

S = POLY2STR(P,'s') or S=POLY2STR(P,'z') POLY2SYM Polynomial coefficient vector to symbolic polynomial. POLY2SYM(C,'V') and

POLY2SYM(C,SYM('V') both use the symbolic variable specified by the second argument.

b) 由根矢量创建多项式

多项式行向量可以由命令poly创建，如A为矩阵，则poly（A）将创建A矩阵的特征多项式，如果A为向量[ ar1 ar2 … arn-1 arn]，则创建（x-ar1）(x- ar2) …(x-arn-1) (x- arn)生成的多项式的系数矢量。

若已知根矢量ar, 通过调用函数 p=poly(ar) 可以产生多项式的系数矢量, 再利用poly2sym函数就可方便的建立符号形式的多项式。例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）根矢量元素为n ，则多项式系数矢量元素为n+1；

（2）函数poly2sym(pa) 把多项式系数矢量表达成符号形式的多项

式，缺省情况下自变量符号为x，可以指定自变量。（3）使用简单绘图函数ezplot可以直接绘制符号形式多项式的曲线。

例 2：由给定复数根矢量求多项式系数矢量。 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

说明：含复数根的根矢量所创建的多项式要注意：

（1）要形成实系数多项式，根矢量中的复数根必须共轭成对；

（2）含复数根的根矢量所创建的多项式系数矢量中，可能带有很小的

虚部，此时可采用取实部的命令real把虚部滤掉。

c) 多项式的求根

如果需要进行系数表示形式的多项式的求根运算，有两种方法可以实现，一是直接调用求根函数roots，poly和 roots 互为逆函数。另一种是先把多项式转化为伴随矩阵，然后再求其特征值，该特征值即是多项式的根。

d) 特征多项式输入法

用poly函数可实现由矩阵的特征多项式系数创建多项式。条件：特征多项式系数矢量的第一个元素必须为1。

例 1：求三阶方阵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阶的。

例 2：将多项式的系数表示形式转换为根表现形式（求根）。求 x3-6x2-72x-27的根 a=[1 -6 -72 -27] r=roots(a) r =

12.1229 -5.7345 -0.3884

MATLAB约定，多项式系数矢量用行矢量表示，根矢量用列矢量表示。

2. 多项式的运算(Computation of polynomial)

a, 多项式的加减法运算(addition and subtraction of polynomial)

MATLAB没有提供多项式直接的加减法函数，若两个多项式向量大小相同，标准的数组加法有效，当两个多项式阶次不同时，低阶多项式必须用首零填补使其与高阶多项式有相同的阶次。例1：a=[3 2 5 8];

b=[3 5 1 9];

c=a+b %同价求和

d=[2 4 1 -4 5];

e=[0 c]+d %不同价求和

c =

6 7 6 17 e =

2 10 8 2 22

减法的操作与加法类似，相当于将后项变量取反。

b,多项式的乘除运算(Multiplication and division of polynomial)

多项式的乘除法对应着向量的卷积与解卷积，多项式乘法（卷积）用函数conv(a,b)实现，除法（解卷积）用函数deconv(a,b)实现。

长度为m的向量a和长度为n的向量b的卷积定义为： C(k)=错误！未找到引用源。 C向量的长度为：m+n-1

解卷积是卷积的逆运算，向量a对向量c进行解卷积将得到商向量 q和余量r，并且满足：错误！未找到引用源。

例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’) %建立指定变量为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

例3：求多项式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

3. 其他常用的多项式运算命令(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分别是分子（numerator）分

母（denominator）多项式系数矢量，r,p,k分别是留数、极点和直项矢量

p=polyfit(x,y,n) 用n阶多项式拟合x，y矢量给定的数据。 polyder(p) 多项式微分。注：对于多项式b(s)与不重根的n阶多项式a(s)之比，其部分分式展开为：

b(s)rrr?1?2?L?n?k(s) a(s)s?p1s?p2s?pn

式中：p1,p2,…,pn称为极点(poles)，r1,r2,…,rn 称为留数(residues)，k(s)称为直项(direct terms)，假如a(s)含有m重根pj，则相应部分应写成：

rjs?pj?rj?1(s?pj)2?L?rj?m?1(s?pj)m

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.

例1：对 (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将是空矩阵，表示无此项)

例2：对一组实验数据进行多项式最小二乘拟合(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 的解，右除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

MATLAB对不同形式的参数矩阵，采用不同的运算法则来处理，它会自动检测参数矩阵，以区别下面几种形式： *三角矩阵(Triangular Matrix)

*对称正定矩阵(symmetrical positive determined matrix) *非奇异方阵(Nonsingular matrix) *超定系统(Overdetermind system) *欠定系统(Underdetermind system)

1. 方阵系统：(Square array) m=n

最常见的方程式为ax=b形式，系数矩阵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) m﹥n 实验数据较多，寻求他们的曲线拟合。如在t内测得一组数据y： t y 0.0 0.82 0.3 0.72 0.8 0.63 1.1 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) m﹤n

欠定系统为线性相关系统，其解都不唯一，MATLAB会计算一组构成通解的基解，而方程的特解则用QR分解法决定。两种解法：最少元素解a\\b，最小范数解pinv(a)*b. 例1：用两种方法求解欠定系统。

对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.

例1：计算方阵的行列式和逆矩阵。 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函数来计算广义逆矩阵。例1：计算广义逆矩阵。 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分解，针对对称正定矩阵；

（２）LU分解，高斯消元法，针对一般矩阵；

（３）QR分解，正交化，针对一般长方形矩阵（行数≠列数）这三种分解运算分别由chol, lu和 qr三个函数来分解. 1． Cholesky分解(Cholesky Decomposition)

仅适用于对称和上三角矩阵，如果A为对称正定矩阵，则Cholesky分解可将矩阵A分解为上三角矩阵和其转置的乘积，即：A=R’*R, R为上三角阵。

方程A*X=b变成R’*R*X=b, 所以有 X=R\\(R’\\b) 例1：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：对下列矩阵进行cholesky分解，并求解方程组的解。 16x1+4x2+8x3=28 4x1+5x2-4x3=5 8x1-4x2+22x3=26

A=[16 4 8;4 5 -4;8 -4 22]; b=[28 5 26]';

R=chol(A) %Cholesky分解

X=R\\(R'\\b) %根据R矩阵求解方程组的解 R =

4 1 2 0 2 -3 0 0 3 X = 1 1 1

2． LU分解(LU factorization).

用lu函数完成LU分解，将矩阵分解为上、下两个三角阵，其调用格式为：

[l,u]=lu(a) l代表下三角阵，u代表上三角阵。例1： 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

\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阶的。例1：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 \

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的特征值。

例1：计算特征值和特征矢量。

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’。

[x,fual,exitflag,output]=fminsearch(fun,x0,options): 单纯形法求多元函数极值点指令的最完整格式。

[x,fual,exitflag,output，grad，hessian]=fminunc(fun,x0,options): 拟牛顿法求多元函数极值点指令的最完整格式。

其中x1,x2分别表示被研究区间的左、右边界。 fminsearch和fminunc中的x0是一个向量，表示极值点位置的初步推测。输入变量options的默认值可以用options=optimset('FunFun_Name')获得。例1，求解函数的极值 %求最小值

fn='2*exp(-x)*sin(x)'; %定义fn函数

xmin=fminbnd(fn,2,5) %在区间（2，5）寻找最小值 emin=5*pi/4-xmin %求最小值的误差 x=xmin; %令x为最小值点

ymin=eval(fn) %计算最小值点的函数值 %求最大值

fx='-(2*exp(-x)+sin(x))'; % 为求最大值，函数值取负 xmax=fminbnd(fx,0,3) %在区间（0，3）寻找最大值 emax=pi/4-xmax %求最大值的误差

ymax=eval(fx) %求最大值点的函数值

x=fminsearch(f,[1;2;3]) %用单纯法求解，从点[1，2，3]开始搜索

4、数值积分：数值积分函数

quad('fname',a,b,tol,trace,p1,p2,…) 低阶法--自适应递推辛普生法 quadl('fname',a,b,tol,trace,p1,p2,…) 高阶法--自适应递推牛顿-柯西法 trapz（x,y）计算变量为x，函数y的积分 fname是被积函数表达式字符串或函数文件名，a,b分别为积分的上、下限，tol用来控制积分精度，默认为tol=0.001,trace 取1，表示用图形表示积分结果，取0则无图形，p1,p2…是向函数传递的参数，可取默认值。例1：用trapz在区间-1＜x＜3上计算y=humps(x)下面的面积 x=-1:0.12:3;

y=humps(x);% Y = HUMPS(X) is a function with strong maxima near

% x = .3 and x = .9.

plot(x,y)

area=trapz(x,y) %以更好的精度计算 x=-1:0.02:3; y=humps(x); area=trapz(x,y)

%用quad和quadl在区间-1＜x＜3上计算y=humps(x)下面的面积 format long

area=quad('humps',-1,3) area=quadl('humps',-1,3)

5. 数值微分与梯度 (Difference and approximate derivative ,gradient). 数值微分函数diff(x) 例1：按列求微分。

x=[1,10,20;2,12,23;3,14,26;3,16,29] d=diff(x) figure(1) plot(x) figure(2)

plot(d) %求一阶微分 x =

1 10 20 2 12 23 3 14 26 3 16 29 d =

1 2 3 1 2 3

0 2 3

数值微分比较困难，函数的微小变化都容易产生相邻斜率的大的变化。因而应尽可能避免数值微分。特别是对实验数据进行微分的时候，最好用最小二次曲线拟合数据，然后对所得的多项式进行微分。

x=0:0.1:1;

y=[-0.443 1.783 3.231 6.21 7.08 7.73 7.98 8.99 9.23 9.93 10.99 ]; n=2;

p=polyfit(x,y,n) %多项式拟合 xi=linspace(0,1,100); z=polyval(p,xi); figure(1) hold on

plot(x,y,'o',x,y,xi,z,':') xlabel('x') ylabel('y')

title('二次曲线拟合')

pd=polyder(p) %多项式微分 w=polyval(pd,xi); plot(xi,w,'*') xlabel('x') ylabel('y')

title('曲线微分')

%函数diff是对一些描述某函数的数据进行粗略计算的微 %分函数，用于计算数组中元素间的差分。f=f(x) 的微分近 %似为 dy/dx=f(x+h)-f(x)/(x+h)-x dy=diff(y)./diff(x); xd=x(1:max(size(x))-1); figure(2) plot(xd,dy)

title('用diff进行近似差分') ylabel('dy/dx轴')

xlabel('x轴')

结果可见用有限差分近似微分导致很差的结果。例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

[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.

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’ 例1：二维插值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.)

例1：三维插值实例。

[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.

*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,Z] = MESHGRID(x,y,z) produces 3-D arrays that can be used to evaluate functions of three variables and 3-D volumetric plots. For example, to evaluate the function x*exp(-x^2-y^2) over the range -2 < x < 2, -2 < y < 2,

[X,Y] = meshgrid(-2:.2:2, -2:.2:2); Z = X .* exp(-X.^2 - Y.^2); mesh(Z)

七．数字信号处理初步(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操作。

注：快速傅立叶变换的结果为复数。

例1：产生一个正弦衰减曲线，进行快速傅立叶变换，并画出幅值(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

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

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 的整数次幂，则运算速度最快；