《通信系统中MATLAB基础与仿真应用》课件-第6章

第6章MATLAB在工程教学中的应用6.1解线性方程组6.2多项式运算6.3曲线拟合6.4插值和样条6.5数值积分和微分6.6数据分析函数和傅立叶变换小结

习题

本章作为本书的基础知识核心部分，从工程教学的角度，系统地介绍了MATLAB在高等数学、线性代数和数据处理等方面的应用，同时给出了大量的实例，并结合图形图像处理和数据可视化等MATLAB功能进行讲述。通过本章的学习，读者可熟练掌握MATLAB关于数学计算方面的内容。需要说明的一点是，本章是本着解决问题的目的进行编写的，因此不再区分此解决方法是属于数值计算还是属于符号计算，无论采取何种计算方法，都是解决该类问题的一种选择。在确实需要区分时，会予以说明。6.1解线性方程组在工程教学和实践中，对线性方程组的求解是一个重要的内容。在实际工程应用、实验数据分析以及进行理论研究等情况下，很多问题都可归结为线性方程组的求解。因而，线性方程组求解的应用非常广泛。求解线性方程组不可避免地要用到矩阵分解的概念，矩阵分解函数在第3章已有简单介绍，本节首先介绍用于解线性方程组的几个矩阵分解函数，然后具体讲解各类线性方程组的解法。6.1.1矩阵的分解

矩阵分解是一个非常重要的概念，如在求方程组的解、矩阵的特征值、特征向量和矩阵的秩等重要参数时，都要用到矩阵的分解。在实际工程中的某些场合也要对矩阵进行特定形式的分解。总之，无论在理论证明或是科学计算上，对矩阵进行分解操作都十分重要。MATLAB提供了许多矩阵分解函数，利用这些函数可以比较容易地进行矩阵分解。下面着重介绍其中几种矩阵分解的方法，其相关函数如表6-1所示。表6-1矩阵分解函数在MATLAB中，线性方程组的求解主要用到三种基本的矩阵分解，即对称正定矩阵的cholesky分解、一般方程的gaussian消去法和矩阵的正交分解，这三种分解由函数chol、lu和qr完成。分解时都使用三角矩阵，其中所有位于对角线以上或以下的元素都为0。矩阵经分解成为三角矩阵后，线性方程组不论是用左除还是右除，都可以简单、快速地求解。

1.对称正定矩阵的cholesky分解并不是所有的对称矩阵都可以进行cholesky分解，能进行此种分解的矩阵必须是正定的，即矩阵X的所有对角线元素都是正的，而非对角线元素不会太大。此种分解可用于复矩阵，此时复矩阵必须满足X′=X是hermitian正定的。设X是一个n×n的对称正定矩阵，则存在对角线为正的上三角矩阵R，使

X＝R′R。在MATLAB中，cholesky分解的函数调用格式有两种:

·R＝chol(X):X是对称正定矩阵，R是上三角矩阵，使R′R＝X。如果矩阵X是非正定的，则会给出错误信息。

·［R，p］=chol(X):返回两个参数，不会给出出错信息。如果X是正定的，则P等于0，R同上；如果X是非正定的，则P等于正整数，R是一个阶数为q=p-1的上三角阵，使得R＇*R=X(1:q，1:q)。下面通过例6-1来说明其用法。

【例6-1】cholesky分解。

X＝［22-2；25-4；-2-45］%输入对称正定矩阵X

R＝chol(X)%进行cholesky分解

X＝

2

2

-2

2

5

-4

-2

-4

5

R＝

1.4142

1.4142

-1.4142

0

1.7321

-1.1547

001.2910

R＇*R

ans＝

2.00002.0000-2.0000

2.00005.0000-4.0000

-2.0000-4.00005.0000

2.lu分解矩阵的lu分解在工程教学和实践中非常重要，许多运算都以lu分解为基础，如方阵的求逆操作inv()、行列式求值det()以及解线性方程组等。lu分解是除法运算的基础。gaussian消去法或lu分解是将任何方阵表示为一个下三角矩阵L和一个上三角矩阵U的乘积。其调用格式有两种:

·［L，U］=lu(X):此调用格式将得到一个上三角矩阵且存储在U中和一个准下三角矩阵且存储在L中，使得X＝LU。准下三角矩阵L实际上是下三角矩阵的转置矩阵。

·［L，U，P］=lu(X):此调用格式将得到一个主对角元素为1的下三角矩阵L、上三角矩阵U和一个由0和1组成的置换矩阵P，使得PX=LU。注意:进行lu分解时，矩阵X必须是方阵。在MATLAB中，lu分解允许线性方程组Ax=b进行如下快速运算:

X=U＼(L＼b)矩阵行列式的值和矩阵求逆也可以利用lu分解进行计算:

det(X)=det(L)*det(U)

inv(X)=inv(U)*inv(L)下面通过例6-2说明如何进行lu分解。

【例6-2】lu分解。[HJ2.3mm]

X＝［3-12；12-1；-214］%输入矩阵X

X＝

3

-1

2

1

2

1

-2

1

4［L，U，P］＝lu(X)%进行lu分解

L＝

1.0000

0

0

0.3333

1.0000

0

-0.6667

0.1429

1.0000

U＝

3.0000

-1.0000

2.0000

0

2.3333

-1.6667

0

0

5.5714

P＝

1

0

0

0

1

0

0

0

1

P*X

ans＝

3

-1

2

1

2

-1

-2

1

4

L*U

ans＝

3.0000

-1.0000

2.0000

1.0000

2.0000

-1.0000

-2.0000

1.0000

4.0000

3.qr分解

qr分解即矩阵的正交分解，它可将矩阵分解为一个正交矩阵和一个上三角矩阵的乘积。此种分解适用于任意矩阵，是非常重要的分解形式。其调用格式为:

·［Q，R］=qr(X):此调用格式生成一个与X同阶的上三角矩阵R和一个正交矩阵Q，使得X＝QR。

·［Q，R，E］=qr(X):此调用格式将得到一个置换矩阵E、上三角矩阵R和正交矩阵Q，使得XE=QR，选择置换矩阵E使得abs(diag(R))递减。

·［Q，R］=qr(X，0):此调用格式将生成一种“经济”的分解，如果矩阵X是m×n阶，并且m＞n，则仅计算出Q的前n列。·［Q，R，E］=qr(X，0)：此调用格式将生成一种“经济”的分解，其中E是一个置换向量，使得QR=A(:，E)，选择列置换向量使得abs(diag(R))递减。

【例6-3】qr分解。

X=［3-125；12-17；-21-24］％输入矩阵X

X=

3

-1

2

5

1

2

-1

7

-2

1

-2

4［Q，R］＝qr(X)%进行qr分解

Q＝

-0.8018

0.1543

0.5754

-0.2673

-0.9567

-0.1155

0.5345

-0.2469

0.8083

R＝

-3.7417

0.8018

-2.4054

-3.7417

0

-2.3146

1.7591

-6.9128

0

0

-0.3464

5.3116［Q，R，E］=qr(X，0)

Q＝

-0.5270

0.6463

0.5518

-0.7379

0.0259

-0.6754

-0.4216

-0.7626

0.4905

R＝

-9.4868

-1.4754

-1.3703

0.5270

0

3.4383

-1.4607

2.8437

0

0

-1.4102

0.7971

E＝

4123

4.奇异值分解奇异值分解在矩阵分析中占有极为重要的位置。奇异值分解的定义为:[HJ2.3mm]对于任意的矩阵A∈Cm×n，存在酉矩阵(UnitaryMatrix)，U=［u1，u2，…，um］，V=［v1，v2，…，vn］使得:

UTAV=diag(σ1，σ2，…，σp)σ1≥σ2≥…≥σp≥0，p=min{m，n}其中，σi、ui、vi分别称为矩阵A的第i个奇异值、左奇异向量和右奇异向量，而它们的组合就称为奇异值分解三对组。这里的上标“T”表示共轭转置。

MATLAB中提供的奇异值分解函数svd的调用格式为:

·［U，S，V］=svd(X):产生一个与矩阵X具有相同维数的矩阵S，其对角线元素为递减的非负值，同时得到酉矩阵U和V，使得X=U*S*V′。

·S=svd(X):得到矩阵X的奇异值组成的向量。

·［U，S，V］=svd(X，0):得到一个“经济大小”的分解结果，如果X是m×n矩阵且m＞n，则只计算U矩阵的前n行，且S矩阵是n×n阶的。

【例6-4】奇异值分解。

X=［314；224；1-3-2；123］%输入矩阵X

X=

3

1

4

2

2

4

1

-3

-2

1

2

3［U，S，V］=svd(X)%进行奇异值分解

U=

0.5873

-0.5075

-0.5954

0.2073

0.5951

-0.0691

0.4057

-0.6903

-0.3132

-0.8424

0.4331

0.06881

0.4503

0.1674

0.5417

0.68981

S=

8.2230

0

0

0

3.2221

0

0

0

0

0

0

0

V=

0.3757

-0.7249

0.5774

0.4400

0.6878

0.5774

0.8157

-0.0371

-0.5774［U，S，V］=svd(X，0)%进行奇异值分解

U=

0.5873

-0.5075

-0.5954

0.5951

-0.0691

0.4057

-0.3132

-0.8424

0.4331

0.4503

0.1674

0.5417

S=

8.2230

0

0

0

3.2221

0

0

0

0

V=

0.3757

-0.7249

0.5774

0.4400

0.6878

0.5774

0.8157

-0.0371

-0.5774奇异值分解也是矩阵求秩运算的基础，对矩阵A进行奇异值分解S＝svd(A)，得到向量S的非零元素的个数就是矩阵A的秩。如:

X=［314；224；1-3-2；123］；

S=svd(X)

S=

8.2230

3.2221

0可见矩阵X的秩为2，用求秩运算rank(X)可以验证这一结果。

5.Schur分解

schur分解的定义为:

A==USUT

其中，A是方阵，U是一个正交矩阵，S是一个块上三角矩阵，由对角线上的1×1和2×2块组成。特征值由S的对角线和块给出，而U的列比一系列特征向量给出了更多的数值特性。schur分解可以对缺陷矩阵进行。矩阵的复schur形式矩阵是一个特征值在对角线上的上三角阵；而实schur形式是实特征值在对角线上，复特征值以2×2的块矩阵排列在对角线上。用函数T=schur(X)求矩阵X的schur形式矩阵时，如果X是实矩阵，则函数返回实数schur形式矩阵；如果X是复矩阵，则函数返回复数schur形式矩阵。函数rsf2csf()可以将实数形式转化为复数形式。矩阵X必须是方阵。函数schur的调用格式为:

·T=schur(X):仅仅返回schur形式矩阵T。

·［U，T］=schur(X):得到schur形式矩阵T和酉矩阵U，使得X＝U*T*U′和U′*U＝eye(size(X))。

【例6-5】schur分解。

A=［1

2

-1；3

4

0；1

2

3］％输入方阵A

A=

1

2

-1

3

4

0

1

2

3［U，T］=schur(A)％进行schur分解

U=

0.8262

-0.2294

-0.5145

-0.5571

-0.4680

-0.6860

0.0835

-0.8534

0.5145

T=

-0.4495

0.8575

0.3760

-0.0000

4.4495

2.8501

0

0

4.0000［UU，TT］=rsf2csf(U，T)％对U和T进行转换

UU=

0.8262

-0.2294

-0.5145

-0.5571

-0.4680

-0.6860

0.0835

-0.8534

0.5145

TT=

-0.44495

0.8575

0.3760

0

4.4495

2.8501

0

0

4.0000

6.矩阵的hessenberg分解

矩阵的hessenberg分解由函数hess完成，该函数命令的调用格式为:

·H=hess(A):可求矩阵的hessenberg形式矩阵H，H的第一子对角线以下的元素为0。如果矩阵A是对称的或是hermitian矩阵，则H是对角三角阵。矩阵A必须是方阵。

·［P，H］=hess(A):得到一个酉矩阵P和一个hessenberg形式矩阵H，使得A＝P*H*P′和P′*P＝eye(slze(A))。

【例6-6】hessenberg分解。

A=［1

2

-1；3

4

0；1

2

3］％输入方阵A

A=

1

2

-1

3

4

0

1

2

3［P，H］=hess(A)％进行hessenberg分解

P=

1.0000

0

0

0

-0.9487

-0.3162

0

-0.3162

0.9487

H=

1.0000

-1.5811

-1.5811

-3.1623

4.5000

0.5000

0

-1.5000

2.50006.1.2线性方程组的求解在工程教学和计算中，线性方程组的求解是一个很重要的问题。在矩阵表示方法中，线性方程组可以表示为给定两个矩阵A和B，求X的唯一解，使得:

AX=B或XA=B在MATLAB中，求解线性方程组时，主要采用前面章节介绍的除法运算符“/”和“＼”求解。如:

·X＝A＼B表示求矩阵方程AX＝B的解；

·X＝B/A表示求矩阵方程XA=B的解。对方程X＝A＼B，要求矩阵A和B有相同的行数，X和B有相同的列数，它的行数等于矩阵A的列数。方程X＝B/A同理。如果矩阵A不是方阵，其维数是m×n，则有:

·当m=n时，为恰定方程，寻求精确解；

·当m>n时，为超定方程，寻求最小二乘解；

·当m<n时，为不定方程，寻求基本解，其中至多有m个非零元素。针对不同的情况，MATLAB将采用不同的算法来求解。6.1.3恰定方程组

恰定方程组由n个未知数的n个方程构成，方程有唯一的一组解。其一般形式可用矩阵、向量写成如下形式:

Ax=b其中，A是方阵，b是一个列向量。在线性代数教科书中，最常用的方程解法有:

·利用cramer公式求解法；

·利用矩阵求逆解法，即x=A-1b；

·利用gaussian消去法；

·利用lu法求解。一般来说，对于维数不高、条件数不大的矩阵，上面四种解法所得的结果差别不大。前三种解法的真正意义是在其理论上，而不是实际的数值计算。在MATLAB中，出于对算法稳定性的考虑，行列式及其逆的计算大都在lu分解的基础上进行。在MATLAB中，求解这类方程组的命令十分简单，直接采用表达式:x=A＼b。

MATLAB的指令解释器在确认变量A非奇异后，就对它进行lu分解，并最终给出解x；若矩阵A的条件数很大，则MATLAB会提醒用户注意所得解的可靠性。如果矩阵A是奇异的，则Ax=b的解不存在，或者存在但不唯一；如果矩阵A接近奇异，则MATLAB将给出警告信息；如果发现A是奇异的，则计算结果为inf，并且给出警告信息；如果矩阵A是病态矩阵，则也会给出警告信息。注意:在求解方程时，尽量不要用inv(A)*b命令，而应采用A＼b的解法。因为后者的计算速度比前者快、精度高，尤其当矩阵A的维数比较大时。另外，除法命令的适用性较强，对于非方阵的A，也能给出最小二乘解。

【例6-7】“左除”法与“求逆”法解恰定方程组时的性能比较。

rand(′state′，12)；％选定随机种子，产生随机矩阵A

A=rand(100，100)＋1.e8；％生成(100×100)均匀分布的随机阵％加大数值是要使矩阵A的条件数增大

X=ones(100，1) ％生成100行1列的向量

b＝A*x；

flops(0) ％浮点运算计数器置0

tic ％启动计时器StopwatchTime，开始计时

y=inv(A)*b； ％求逆法解方程运算次数

tic=toc ％关闭计时器，并显示解方程所用的时间

ci=flops ％求逆法解方程运算次数

erri=norm(y-x) ％结果与精确解的范-2误差

errif=norm(A*y-b) ％方程的范-2误差％％％

flops(0)

tic

y=A＼b； ％用除法解方程

td=toc ％关闭计时器，并显示解方程所用的时间

cd=flops ％除法解方程运算次数

errd=norm(y-x) ％结果与精确解的范-2误差

errdf=norm(A*y-b) ％方程的范-2误差运算结果:

ti=

0.0600 ％求逆法解方程所用的时间

ci=

2070322 ％求逆法解方程所用的运算次数

erri=

3.0708e-004

errif=

6.6280e+004

td= ％除法解方程所用的时间

0

cd= ％除法解方程所用的运算次数

741872

errd=

3.2243e-004

errdf=

2.0095e-005

从本例执行的结果可知，求逆法解方程所需运算次数是除法解方程的2.5倍，时间上除法几乎是“机器零”。6.1.4超定方程组

对于方程组Ax=b，A为n×m矩阵，如果A列满秩，且n>m，则方程没有精确解，此时称方程组为超定方程组。线性超定方程经常遇到的问题是数据的曲线拟合。对于超定方程，在MATLAB中，利用左除命令(x=A＼b)来寻求它的最小二乘解；还可以用广义逆来求解，即x=pinv(A)，所得的解不一定满足Ax=b，x只是最小二乘意义上的解。左除的方法建立在奇异值分解基础之上，由此获得的解最可靠；广义逆法建立在对原超定方程直接进行householder变换的基础上，其算法可靠性稍逊于奇异值分解，但速度较快。【例6-8】求超定方程组的解

A=［2

-1

3；3

1

-5；4

-1

1；1

3-1

3］％输入矩阵A

A=

2

-1

3

3

1

-5

4

-1

1

1

3-

13

b=［3

0

3

-6］′；

rank(A)

ans=

3

x1=A＼b％左除解方程

x2=pinv(A)b ％广义逆求解

x1=

1.0000

2.0000

1.0000

x2=

1.0000

2.0000

1.0000

A*x1-b

ans=

1.0e-014*

-0.0888

-0.0888

-0.1332

0可见x1并不是方程Ax=b的精确解，用x2=pinv(A)*b所得的解与x1相同。6.1.5欠定方程组

欠定方程组未知量个数多于方程个数，但理论上有无穷个解。MATLAB将寻求一个基本解，其中最多只能有m个非零元素。特解由列主元qr分解求得。【例6-9】解欠定方程

A=［1-211；1-21-1；1-215］ ％输入矩阵A

A=

1

-2

1

1

1

-2

1

-1

1

-2

1

5

b=［1-15］＇；

x1=A＼b％左除法解方程

Warning:Rankdeficient，rank=2tol=4.6151e-015

x1=

0

-0.0000

0

1.0000

x2=pinv(A)*b ％用广义逆解方程

x2=

0

-0.0000

0.0000

1.00006.1.6方程组的非负最小二乘解在某种情况下，所求的线性方程组的解出现负数是没有意义的。虽然方程组可以得到精确解，但却不能取负值解。在这种情况下，其非负最小二乘解比方程的精确解更有意义。在MATLAB中，求非负最小二乘解通常用函数nnls，其调用格式为:

·X=nnls(A，b):返回方程Ax=b的最小二乘解，方程的求解过程被限制在x≥0的条件下。

·X=nnls(A，b，TOL):指定误差TOL来求解，TOL的默认值为TOL=max(size(A))*norm(A，1)*eps，矩阵的-1范数越大，求解的误差越大。

·［X，W］=nnls(A，b):当x(i)=0时，w(i)<0；当x(i)>0时，w(i)≈0，同时返回一个双向量w。

【例6-10】求方程组的非负最小二乘解。

A=［3.4336 -0.5238

0.6710

-0.5238 3.2833

-0.7302

0.6710 -0.7302

4.0261］；％输入矩阵A

b=［-1.0000 1.5000

2.5000］；［X，W］=nnls(A，b)％求方程的非负最小二乘解

X=

0

0.6563

0.6998

W=

-3.6820

-0.0000

-0.0000

x1=A＼b ％用除法解方程

x1=

-0.3569

0.5744

0.7846

A*X-b ％验证非负最小二乘解

ans=

1.1258

0.1437

-0.616

A*x1-b ％验证除法解

ans=

1.0e-015*

-0.2220

0.4441

06.1.7方程解的精度人们在用计算机计算一个问题时，最关心的是该问题的数值解。而所得的解是否可靠，除非计算所用的是一些特殊的整数或有理数，否则计算中的圆整误差、截断误差等都不可避免。在MATLAB中，数与数之间的最小分辨率用eps表示。表达式中任何数的相对误差都不可能小于eps。此外，在程序执行过程中，都不可避免地使这个最小误差比放大。范数和条件数对求解过程中误差放大现象有重要的定量描述，在这里仅重点说明方程的精度问题。在数值分析中，方程精度的估计可利用解的相对误差来进行，其公式为其中，K是矩阵的条件数，在MATLAB中的命令为cond()。由于eps是机器精度，所以可以用K*eps的大小粗略判断所得的方程解是否可靠。

【例6-11】对方程Ax=b的近似解和精确解进行比较，其中A是Hilbert矩阵。

N=［6，14］；％计算的矩阵阶数

fori=1:length(N)

n=N(i) ％所要计算的矩阵阶数放入n中

A=hilb(n)； ％产生n阶Hilbert矩阵

Ai=invhilb(n)； ％产生完全准确的n阶逆Hilbert矩阵

b=ones(n，1)； ％产生n阶全是1的向量

x_app=A＼b； ％利用左除来求近似解

x_exa=Ai*b； ％利用准确求逆来求方程的准确解

fdb=norm(A*x_app-b)；

fb=norm(b)；

fdx=norm(x_app-x_exa)；

fx=norm(x_app)；

err_actual(i)=fdx/fx； ％实际的相对误差

K=cond(A)；

err_app(i)=K*eps； ％最大可能的近似相对误差

err_max(i)=K*fdb/fb； ％最大可能的相对误差

end

disp(′Hilbert矩阵阶数′)，disp(N)

formatshorte

disp(′实际的相对误差′)，disp(err_actual)，disp(′′)

disp(′最大可能的近似相对误差′)，disp(err_app)，disp(′′)

disp(′最大可能的相对误差′)，disp(err_max)，disp(′′)运行结果如下:

Hilbert矩阵阶数

6 14实际的相对误差

5.0339e-011 3.9641e+000最大可能的近似相对误差

3.3198e-009 9.4718e+001最大可能的相对误差

6.0095e-007 6.6239e+009

从该例可以看出，14阶的Hilbert矩阵是严重错误的。6.1.8用函数零点求方程的解

在MATLAB中，用函数零点法求方程或方程组解的函数命令有两个，即fzero和fsolve。在此，首先需对方程和方程组进行转化，比如将方程f(x)=g(x)转化为F(x)=f(x)-g(x)=0，方程组也是如此；然后将函数F(x)写成MATLAB的m函数，以便在fzero和fsolve命令中调用。本小节重点讲述其使用方法。

1.一元方程转化的函数零点求法在一元方程中，任意函数方程F(x)=0可能有零点，也可能没有零点；可能有一个零点，也可能有多个零点，很难有一个通用的解法。一般来说，其求解的过程为:先猜测一个初始零点或者该零点大概所在的区间，然后通过计算使猜测值不断精确化，或使猜测区间不断收缩，直至达到预先指定的精度为止。在猜测一个初始的零点时，一般用MATLAB的作图命令来获取初始近似解。其具体步骤为:先确定一个零点可能存在的自变量区间，然后利用fplot命令绘出f(x)在该区间的图形，用眼观察F(x)与横轴的交点坐标，或者更细些用zoom命令对交点局部放大来读数，或借助ginput命令来获得更精确的交点坐标。

1)利用fzero函数求一元函数零点命令fzero的调用格式为:

(1)x=fzero(fun，x0):求一元函数零点命令的最简形式；

(2)［x，fval，exitflag］=fzero(fun，x0，options，P1，P2，…):从MATLAB5.3版本起，求一元函数零点命令最完整的格式。其中参数x0是初始猜测的零点，同时也是预定搜索零点的大致位置。它可以是标量或是二元向量，当x0是标量时，该命令将在它的两侧寻找一个与之最靠近的零点:当x0是一个二元向量［a，b］时，该命令将在区间［a，b］内寻找零点。options是优化迭代所采用的参数选项。它是MATLAB设计优化程序fzero、fsolve、fminbnd、fminsearch和fminunc时都需要的一个“模块”，采用结构数组存放优化参数。如果没有优化参数要设置，则可以在optons的位置上用“［］”作为占位符。在此命令中，options的缺省设置可以用命令options=optimset(′fzero′)获得。它有两个域:Display和TolX。其中:

options.Display:显示设置，有三个选取值［off|iter|{final}］；

options.TolX:自变量计算的终止误差，可选取［opsitivescalar］，缺省值为2.2204e-016。若要修改缺省值，可用options.TolX=0.001或用options=optim-set(′TolX′，0.001)来改动。

P1和P2是向函数fun传递的附加参数。它的具体取名和函数fun中一致。

x是输出参数，为所求的零点自变量值。

fval是输出参数，为函数fun在x处的值。

exitflag用于描述函数fun的退出情况，表明程序的终止条件。若exitflag>0，则表示找到函数零点后退出；若exitflag<0，则表示没有找到零点或在搜索过程中遇到了无穷大的函数值。

【例6-12】试用fzero命令求解函数f(x)=x4-4x-5的零点。在MATLAB内，建立M文件C6L12。

(1)建立函数f(x)M文件。

funtiony=fun1(x)%例6-12用M函数文件

y=x.^4-4*x-5；

(2)建立水平横轴的M文件。

functiony=fun2(x) %例6-12用M函数文件(绘出水平轴的函数)

y=0；

(3)用作图法估计函数的零点位置。

fplot(′fun1′，［-5，5］)

holdon

fplot(′fun2′，［-5，5］，′r′)程序的运行结果如图6-1所示。图6-1函数零点估计图

2)用zoom和ginput命令获得零点的初始近似值在MATLAB的绘图窗口中，可实现zoom和ginput命令进行图形和鼠标的交互操作，或直接在程序中输入下列命令，便可得到如图6-2所示的局部放大图及鼠标操作线。

zoomon％局部放大命令［tt，yy］=ginput(2) ％用鼠标获取2个零点猜测值

zoomoff ％恢复原来图形大小显示所得零点初始猜测值，结果为:

tt=

-1.0138

1.9124

yy=

0.5848

0.5848图6-2局部放大图和利用鼠标取值图

3)用函数fzero求函数的精确零点［x，fval，exitflag］=fzero(′fun1′，tt(1)，［］)％靠近tt(1)点处的精确零点［x，fval，exitflag］=fzero(′fun1′，tt(2)，［］) ％靠近tt(2)处的精确零点

结果为:

Zerofoundintheinterval:［-0.98515，-1.0138］

x=

-1

fval=

0

fxitflag=

1

Zerofoundintheinterval:［1.8584，1.9124］

x=

1.8812

fval=

-6.2172e-015

exitflag=

1

2.多元函数的零点解非线性方程组

对于采用多元函数的零点解非线性方程组，MATLAB提供有fsolve命令来完成这类工作。一般来说，多元函数的零点问题更难解决。而一旦知道零点的大致位置，就可用数值方法搜索到精确零点。同一元函数用函数曲线与横轴的交点来获得零点的大致位置一样，二元函数的零等位线有助于问题的解决。但更高维的初始零点则很难通过图形获得。有了初始零点后，即可借助函数fsolve命令来求零点的精确解。其调用格式为:

·x=fsolve(fun，x0):解非线性方程组最简单的调用格式。该式中除两个输入参数外，其余输入输出参数都可以缺省。

·［x，fval，exitflag，output，jacob］=fsolve(fun，x0，options，P1，P2…):解非线性方程组最完整的调用格式。各参数的含义如下:

·x0：零点数是猜测值的向量。

·options:fsolve的优化迭代所采用参数的结构数组。它的缺省值可用options=optimset(′fsolve′)获得。它含有5个域:Display、MaxFunEvals、Maxlter、TolFun和TolX。其关键参数为options.Display。

options.Display:显示设置，有3个选取值［off|iter|{final}］；

options.MaxfunEvals:允许函数计算的最多次数，可选取［positiveinteger］，缺省值是自变量的100倍；

options.Maxlter:允许的最多迭代次数，可选取［positiveinteger］，缺省值为400；

options.TolFun:函数计算的终止误差，可选取［positivescalar］，缺省值为1.0000e-006；

options.TolX:自变量计算的终止误差，可选取［positivescalar］，缺省值为1.0000e-006。

·P1和P2:向函数fun传递的参数。

·x和fval:输出参数，分别是所求零点的自变量值和函数值。·exitflag:若exitflag>0，则表示找到零点后退出；若exitflag<0，则表示没有找到零点或在搜索过程中遇到了无穷大的函数值。

·output:输出本命令所用的计算方法、迭代次数等信息，它是结构数组。

·jacob:函数在x处的jacobian。注意:上两例中输入参数fun有三种编写方式:字符串方式、内联函数方式和M函数文件方式。

【例6-13】试求方程组的解。

(1)观察量函数0等位线的交点情况。

x=-1:0.5:1；

y=x；［X，Y］=meshgrid(x，y)；％产生x-y平面上的网点坐标

fun1=sin(X)+Y； ％函数fun1

fun2=X+6*Y ； ％函数fun2

v=［-0.2，0，0.2］ ； ％指定三个等位线，是为了更可靠地判断0等位线的位置

contour(X，Y，fun1，v) ％绘出fun1的三条等位线

holdon

contour(X，Y，fun2，v) ％绘出fun2的三条等位线

holdoff

(2)从图6-3中用鼠标获取零点的初始近似值。［x0，y0］=ginput(1)； ％用鼠标在图形上取一个点的坐标(三线组中间那条线)

disp(［x0，y0］) ％显示初始零点猜测值

0.0000 -0.0058图6-3两个二元函数0等位线的交点图有了初始零点后，便可利用fsolve函数命令求精确解。

fun=′［sin(x(1))+x(2)，x(1)+6*x(2)］′%用字符串表达式形式命令。注意自变量必须写成x(1)和x(2)

XXYY=fsolve(fun，［x0(1)，y0(1)］)%解此非线性方程组

disp(XXYY)

XXYY=

1.0e-016*

0.169 -0.0347

1.0e-016*

0.1691-0.0347检验:

FF1=sin(XXYY(1))+XXYY(2) ％检验方程1

FF2=XXYY(1)+6*XXYY(2) ％检验方程2

FF1=

1.3444e-017

FF2=

-3.9031e-018上面的FF1和FF2表明，所求的零点相当准确。说明:在使用函数fsolve命令时，fun也可是内联函数的形式，此时

fun=inline(′［sin(x(1))+x(2)，x(1)+6*x(2)］′，′x′)；%′x′项必须有

XXYY=fsolve(fun，［x0(1)，y0(1)］)

fun函数也可用M函数文件的形式，此时应先用fun.m表示被解函数(并在搜索路径上):

functionyy=fun(x)

yy(1)=sin(x(1))+x(2)′

yy(2)=x(1)+6*x(2)；

XXYY=fsolve(′fun′，［x0(1)，y0(1)］)6.1.9符号方程及方程组的求解

在工程教学或工程实践中，会遇到一些带有符号的方程，或某些问题需求它的符号解，这类问题就是符号方程和方程组的求解问题。一般在实际中遇到的代数方程组包括线性方程组、非线性方程组和超越方程组。在MATLAB中，用于求解代数方程组的命令有linsolve命令和solve命令。

1.线性方程组的符号解

矩阵计算是求解线性方程组最简单有效的方法。从MATLAB5.x版本起，不管数据对象是数值还是符号，实现矩阵运算的指令几乎相同。因而，求线性方程组符号解时，可以套用数值解的命令的编写方法进行。而在MATLAB中，函数命令linsolve专门用来求解线性方程组，其使用等同于数值计算方法。对方程A*X=B，函数命令linsolve的调用格式为:

X=linsolve(A，B)等同于

X=sym(A)＼sym(B)矩阵A必须至少是行满秩的。当A的列数大于行数时，将给出解不唯一的警告提示。使用该格式可得到方程组的特解X。方程组的通解XX为

XX==X+′k′*null(A)其中，k是任意常数；null命令将求出A的“零空间”的基。

【例6-14】求给定线性方程组的解。

A=sym(′［12/31/2；351；121］′) %输入矩阵A

A=［1，2/3，1/2］［3，5，1］［1，2，1］

b=sym(′［12；1/33；11/7］′)；

X=linsolve(A，b)%计算线性方程的解

X=［23/39，283/91］［-8/13，-102/91］［64/39，-66/91］

【例6-15】求欠定方程组。

A=sym(′［12/31/21；3511；1211］′)％输入矩阵A

A=［1，2/3，1/2，1］［3，5，1，1］［1，2，1，1］

b=sym(′［1；1/3；1］′)；

X=linsolve(A，b) %特解

X=［-1/3］［0］［0］［4/3］

XX=X+′k′null(A)%通解

XX=［-1/3-3/2*k］［k］［-8/3*k］［4/3+13/6*k］说明:欠定方程组的解不唯一，而MATLAB总是给出一个特解，并且解向量中含非零元素最少。

2.一般代数方程组的符号解

slove命令可以解一般代数方程，包括线性方程、非线性方程和超越方程。当方程不存在符号解，且又无其他自由参数时，函数solve将给出数值解。其命令调用格式为:

·solve(′eqn1′，′eqn2′，…，′eqnN′):对N个方程的默认变量求解。

·solve(′eqn1′，′eqn2′，…，′eqnN′，′var1，var2，…，varN′):对N个方程的var1，var2，…，varN变量求解，但该式要注意变量的英文字母顺序，并且在变量前不可有空格。

·S=solve(′eqn1′，′eqn2′，…，′eqnN′，′var1′，′var2′，…，′varN′):对N个方程的′var1′，′var2′，…，′varN′变量求解。

·［x1，x2，…，xn］=solve(′eqn1′，′eqn2′，…，′eqnN′，′var1′，′var2′，…，′varN′):对变量

var1，var2，…，varN求解，并且求解的结果分别赋给x1，x2，…，xn。此式中，MATLAB是按照变量var1，var2，…，varN在英文字母中的顺序赋值给x1，x2，…，xn的。提示:′eqn1′，′eqn2′，…，′eqnN′是字符串表达的方程，或是字符串表达式；′eqnN′，′var1′，′var2′，…，′varN′是字符串表达的求解变量名。第三个命令中，S是一个结构数组，如果要显示结果，则必须使用S.var1，S.var2，…，S.varN的引用形式。在得不到“封闭性解析解”时，将给出数值解。

【例6-16】求非线性方程组的解。［x，y，z］=solve(′x^2+sqrt(2)*x+2=0′，′x+3*z=4′，′y*z=-1′，′x′，′y′，′z′)

x=［(-1/2+1/2*i*3^(1/2))*2^(1/2)］［(-1/2-1/2*i*3^(1/2))*2^(1/2)］

y=［-51/73+3/73*i*3^(1/2)-27/146*(-1/2+1/2*i*3^(1/2))*2^(1/2)-3/146*2^(1/2)］［-51/73-3/73*i*3^(1/2)-27/146*(-1/2-1/2*i*3^(1/2))*2^(1/2)-3/146*2^(1/2)］

z=［-1/3*(-1/2+1/2*i*3^(1/2))*2^(1/2)+4/3］［-1/3*(-1/2-1/2*i*3^(1/2))*2^(1/2)+4/3］

S=solve(′x^2+sqrt(2)*x+2=0′，′x+3*z=4′，′y*z=-1′，′x′，′y′，′z′)

S=

x:［2x1sym］

y:［2x1sym］

z:［2x1sym］

disp(′S.x′)，disp(S.x)％显示结构数组S中x的内容

S.x［(-1/2+1/2*i*3^(1/2))*2^(1/2)］［(-1/2-1/2*i*3^(1/2))*2^(1/2)］

disp(′S.y′)，disp(S.y)％显示结构数组S中y的内容

S.y［-51/73+3/73*i*3^(1/2)-27/146*(-1/2+1/2*i*3^(1/2))*2^(1/2)-3/146*2^(1/2)］［-51/73-3/73*i*3^(1/2)-27/146*(-1/2-1/2*i*3^(1/2))*2^(1/2)-3/146*2^(1/2)］

disp(′S.z′)，disp(S.z)％显示结构数组S中z的内容

S.z［-1/3*(-1/2+1/2*i*3^(1/2))*2^(1/2)+4/3］［-1/3*(-1/2-1/2*i*3^(1/2))*2^(1/2)+4/3］［x，y］=solve(′x^x=2′，′x/y=3′)

x=

log(2)/lambertw(log(2))

y=

1/3*log(2)/lambertw(log(2))上述符号解中的lambertw(ω)代表ω函数ω(x)eω(x)=x的解。有关详细内容，请用mhelp命令访问MAPLE库的lambertw条目。

【例6-17】解超越方程组的解。6.1.10矩阵的特征值和特征向量

特征值和特征向量是线性代数中非常重要的概念，在实际的工程应用和求解数学问题中占有非常重要的地位。比如工程中的振动问题、稳定问题等，从数学关系上常常归结为矩阵的特征值和特征向量的求解问题。在求解微分方程组以及简化矩阵时都要用到特征理论。而在与振动有关的各学科中，特征值和特征向量的问题都有着广泛的应用。矩阵A与向量x相乘，即表示矩阵对向量的变换(Transformation)。对于任何一个矩阵来说，总存在一些特殊的向量，在变换的作用下，向量的方向不变，仅是长短发生变化。这种向量就是所谓的特征向量(Eigenvector)。它满足方程:Ax=λx其中，A是n×n的方阵，λ称为特征值(Eigenvalue)，是个标量。x是对应λ的一个长度为n的列向量。该方程就称为特征方程(EigenvalueEquation)。上式的广义特征方程是:Ax=λBx。此式中，A和B都是n×n的矩阵，λ是标量。当B＝1时，就退化为Ax=λx。而广义特征值问题就是方程Ax=λBx的非平凡解问题。一些有缺陷的矩阵不能对角化，不能进行特征值分解。由A的特征值构成的对角矩阵以及由对应的特征向量构成的矩阵V的各列，必须满足:AX＝VD。如果V是非奇异的，则这就是矩阵A的特征值分解。

1.矩阵的数值特征值和数值特征向量在数学教科书上，最常见的求解特征方程Ax=λx的方法是:先根据|A-λxi|=0求特征值λxi(i＝1，2，…，n)，然后由Ax=λx求对应的特征向量xi。而在MATLAB中，其计算特征值和特征向量的算法取自EISPACK程序库，相应的计算命令较为简单，调用格式如下:

·D＝eig(A):仅计算A的特征值组成的列向量。

·［V，D］=eig(A):生成由矩阵A的特征值、矩阵D和特征向量构成的矩阵V，使A*V=V*D。

·［V，D］=eig(A，′nobalance′):计算时不采用预先平衡。通常，预先平衡增加了特征值和特征向量的计算精度，然而，当一个矩阵包含有与截断误差数量级相差不远的元素时，平衡过程有可能将它们放大，从而导致错误的特征值。该指令可使精度增加。

·D＝eig(A，B):如果A和B是方阵，则生成广义特征值D。·［V，D］=eig(A，B):计算广义特征值矩阵D和广义特征向量矩阵V，使得A*V=B*V*D

【例6-18】计算矩阵A的特征值和特征向量。

A=［-110；-430；102］；

B=［123；139；402］；

D＝eig(A)％计算矩阵A的特征值

D＝

2

1

1［V，D］＝eig(A)

V＝

0

0.4082

-0.4082

0

0.8165

-0.8165

1.0000

-0.4082

0.4082

D＝

2

0

0

0

1

1

0

0

1［V，D］＝eig(A，B)

V＝

0.5940-0.0854i

-0.0356+0.5465i

-0.5983-0.0466i

-0.2068+0.3408

-0.0538+0.8255i

0.2285+0.3267i

-0.4748+0.5054i

0.0081-0.1251i

0.4410-0.5352i

D＝

0.2524-0.5295i

0

0

0

0.1531-0.0000i

0

0

0

0.2524+0.5292i

A*V

ans＝

-0.8008+0.4262i

-0.0182+0.2790i

0.8296+0.3733i

-2.9965+1.3640i

-0.0189+0.2905i

3.0789+1.1663i

-0.3556+1.0962i

-0.0193+0.2964i

0.2836-1.1170i

B*V*D

ans＝

-0.8008+0.4262i

-0.0182+0.2790i

0.8269+0.3733i

-2.9965+1.3640i

-0.0189+0.2905i

3.0789+1.1663i

-0.3556+1.0962i

-0.0193+0.2964i

0.2836-1.1170i

2.符号特征值和特征向量

在MATLAB中，也是使用函数eig计算方阵A的符号特征值和特征向量。其调用格式为:

Lambda＝eig(A)［V，D］=eig(A)计算任意精度的矩阵特征值和特征向量的调用格式为:

Lambds＝eig(vpa(A))［V，D］=eig(vpa(A))上式中，各参数的具体含义与数值特征值和特征向量相同，在此不再赘述。

【例6-19】计算符号特征值和特征向量。

A＝［8/91/21/3；1/21/31/4；1/31/41/5］；

A＝sym(A)%转化为符号矩阵

A＝［8/9，1/2，1/3］［1/2，1/3，1/4］［1/3，1/4，1/5］［V，D］＝eig(A)

V＝［28/153+2/153*12589^(1/2)，28/153-2/153*12589^(1/2)，1］［1，1，-4］［292/255-1/255*12589^(1/2)，299/255+1/255*12589^(1/2)，10/3］

D＝［32/45+1/180*12589^(1/2)，0，0］［0，32/45-1/180*12589^(1/2)，0］［0，0，0］

eig(vpa(A))

ans＝［-.31796273547608369503007e-31］［.87773816609932252988056841486003e-1］［1.3344484056122899692341653807363］［V，D］=eig(vpa(A))

V＝［.18860838403857944292978827467058，.80318501188318836927227766194002，-.56508469644519509785929796275283］［.75443353615431777171915309868189，.48687247435830155621160586316434，.44021044199100581223229247499660］［.62869461346193147643262758223468，.34329146566500535667224373721531，.69777793932276550403654035882202］

D=［-.33889316547608368542146e-31，0，0］［0，1.3344484056122899692341653907363，0］［0，0，.087773816609932252988056841485995e-1］6.1.11矩阵的对角化和其他矩阵函数矩阵运算经过矩阵对角化后变得十分简单，特别对于维数较多的矩阵更是如此。矩阵对角化在实际的工程教学和实践中有着非常广泛的用途。几乎所有的对角化都基于特征值和特征向量的求解，特征值和特征向量的求解目的也是为了对角化。

1.矩阵的PAP对角化由线性代数可知，对于任意可对角化的矩阵A，都存在一个可逆矩阵P，使得P-1AP为对角阵，并且对角阵的对角线元素为矩阵A的特征值。在MATLAB中，求出矩阵特征值和特征向量D和V后，即可满足上述条件(P＝V)。

【例6-20】矩阵的PAP对角化。

A＝［32-1；-2-22；361］；［P，D］＝eig(A)

P＝

0.8890

0.2673

-0.0531

-0.2540

-0.5345

0.4677

0.3810

0.8081

0.8823

D＝

2.0000

0

0

0

-4.0000

0

0

02.000

0

Inv(P)*A*P

ans＝

2.0000

0.0000

-0.0000

-0.0000

-4.0000

0.0000

-0.0000

-0.0000

2.0000

2.实对称矩阵的QRQ对角化

实对称矩阵A都可对角化，并且都存在正交矩阵Q，使得Q-1AQ即(Q′AQ)为对角阵。对角阵的对角线元素均为矩阵A的特征值。由线性代数的知识可知，对于实对称矩阵A，特征值分解函数eig(A)返回的特征向量矩阵就是正交矩阵。

【例6-21】实矩阵的QRQ对角化。

A=［2

2

-2；2

5

-；-2

-4

5］；［Q，D］＝eig(A)

Q＝

0.8944

0.3333

-0.2981

-0.4472

0.6667

-0.5963

0-0.6667-0.7454

D＝

1.0000

0

0

0

10.0000

0

0

0

1.0000

Q*A*Q

ans＝

1.0000

0

0

0

10.0000

0

-0.0000

-0.0000

1.0000

inv(Q)*A*Q

ans＝

1.0000

-0.00000.0000

-0.0000

10.0000-0.0000

0.0000

-0.0000

1.0000

3.约旦(Jordan)标准型矩阵当用相似变换对角化矩阵时，就会产生约旦标准型矩阵。即给定矩阵A，寻找非奇异矩阵V使inv(V)*A*V(更简洁地说，是使J=V＼A*V)尽可能地接近对角阵。约旦标准型是特征值矩阵，它的变换矩阵的每一列就是特征向量。对于有重特征值的非对称矩阵，不能进行对角化。约旦标准型对角线上的元素是特征值，但一些对角线以上的元素是1，而不是0。在MATLB中，计算约旦标准型的函数是jordan，其调用格式为:

J=jordan［V，J］=Jordan(A)式中:J是约旦标准型，V是相似变换矩阵，使得V＼A*V=J。V的列是A的一般化特征向量。

【例6-22】求矩阵的约旦标准型。

A=［111；021；131］；［V，J］=Jordan(A)

V＝

0.7500

0.0947

0.1553

-0.2500

0.0947

0.1553

0.2500

0.1479

-0.3979

J＝

1.0000

0

0

0

3.5616

0

0

0

-0.5616

d＝eig(A)