MATLAB实验指导书

《MATLAB实验指导书》由会员分享，可在线阅读，更多相关《MATLAB实验指导书（27页珍藏版）》请在装配图网上搜索。

1、MATLAB实验指导书线性代数部分北京石油化工学院、基础知识1.1 常见数学函数函数名数学计算功能函数名数学计算功能abs (x)实数的绝对值或复数的幅值floor (x)对x朝-g方向取整acos (x)反余弦arcsin xgcd (m, n)求正整数m和n的最大公约数acosh (x)反双曲余弦arccosh ximag (x)求复数x的虚部angle (x)在四象限内求复数x的相角lcm (m, n)求正整数m和n的最小公倍数asin (x)反正弦arcsin xlog (x)自然对数（以e为底数）asinh (x)反双曲正弦arcsinh xlog10 (x)常用对数（以10为底数）

2、atan (x)反正切arctan xreal (x)求复数x的实部atan2 (x,y)在四象限内求反正切rem (m, n)求正整数m和n的m/n之余数atanh (x)反双曲正切arctanh xround (x)对x四舍五入到最接近的整数ceil (x)对x朝+e方向取整sign (x)符号函数：求出x的符号conj (x)求复数x的共轭复数sin (x)正弦sin xcos (x)余弦cos xsinh (x)反双曲正弦sinh xcosh (x)双曲余弦cosh xsqrt (x)求实数x的平方根：Jxexp (x)指数函数extan (x)正切tan xfix (x)对x朝原点方

3、向取整tanh (x)双曲正切tanh x女口：输入 x=-4.85 -2.3 -0.21.3 4.566.75,则:ceil(x)= -4-20257fix(x) = -4-20146floor(x) = -5-3-1146round(x) =-5-201571.2 系统的在线帮助1 help 命令：1. 当不知系统有何帮助内容时，可直接输入help以寻求帮助： help （回车）2. 当想了解某一主题的内容时，如输入： help syntax （了解 Matlab 的语法规定）3. 当想了解某一具体的函数或命令的帮助信息时，如输入： help sqrt （了解函数 sqrt 的相关信息）2

4、 lookfor 命令现需要完成某一具体操作，不知有何命令或函数可以完成，如输入 lookfor line （查找与直线、线性问题有关的函数）1.3 常量与变量系统的变量命名规则：变量名区分字母大小写；变量名必须以字母打头，其后 可以是任意字母，数字，或下划线的组合。此外，系统内部预先定义了几个有特殊意 义和用途的变量，见下表：特殊的变量、常量取值ans用于结果的缺省变量名pi圆周率n的近似值（3.1416）eps数学中无穷小（epsilon）的近似值（2.2204e - 016）inf无穷大，如 1/0 = inf ( infinity)NaN非数，如 0/0 = NaN (Not a Nu

5、mber)， inf / inf = NaNi，j虚数单位：i = j八11 数值型向量（矩阵）的输入1任何矩阵（向量），可以直接按行方式输入每个元素：同一行中的元素用逗号（，） 或者用空格符来分隔；行与行之间用分号（；）分隔。所有元素处于一方括号（ ） 内；例 1 ： Time = 11 12 1 2 3 4 5 6 7 8 9 10 X_Data = 2.32 3.43； 4.37 5.982系统中提供了多个命令用于输入特殊的矩阵：函数功能函数功能compan伴随阵toeplitzToeplitz 矩阵diag对角阵vanderVandermonde 矩阵hadamardHadamard

6、矩阵zeros元素全为0的矩阵hankelHankel矩阵ones元素全为1的矩阵invhilbHilbert矩阵的逆阵rand元素服从均匀分布的随机矩阵kronKronercker 张量积randn元素服从正态分布的随机矩阵magic魔方矩阵eye对角线上元素为1的矩阵pascalPascal矩阵meshgrid由两个向量生成的矩阵上面函数的具体用法，可以用帮助命令help得到。如:meshgrid(x,y)输入 x=1 2 3 4;y=1 0 5; X,Y=meshgrid(x, y)，则X =Y =12341111123400001 2 3 4 5 5 5 5 目的是将原始数据X, y转

7、化为矩阵数据X, Y。2 符号向量（矩阵）的输入1用函数 sym 定义符号矩阵：函数 sym 实际是在定义一个符号表达式,这时的符号矩阵中的元素可以是任何的符号 或者是表达式,而且长度没有限制。只需将方括号置于单引号中。例 2 ： sym_matrix = sym（a b c；Jack Help_Me NO_WAY） sym_matriX = a, b, cJack, Help_Me, NO_WAY2用函数 syms 定义符号矩阵 先定义矩阵中的每一个元素为一个符号变量,而后像普通矩阵一样输入符号矩阵。例 3 ： syms a b c ； M1 = sym（Classical）； M2 = s

8、ym（ Jazz）； M3 = sym（Blues）； A = a b c； M1, M2, M3；sym（2 3 5）A =a,b,cClassical,Jazz,Blues2,3,51.4 数组（矩阵）的点运算运算符：+ （加）、-（减）、./ （右除）、. （左除）、人（乘方），例 4 ： g = 1 2 3 4；h = 4 3 2 1； si = g + h, s2 = g.*h, s3 = g.Ah,s4 = g.A2, s5 = 2.Ah1.5 矩阵的运算运算符：+ （加）、-（减）* （乘）/ （右除） （左除）人（乘方）、（转置、等； 常用函数：det （行列式）、inv （逆

9、矩阵）、rank （秩）、eig （特征值、特征向量）、rref （化矩阵为行最简形）例 5 ： A=2 0-1;1 3 2; B=l 7-1;4 2 3;2 0 1; M = A*B%矩阵A与B按矩阵运算相乘 det_B = det(B)%矩阵 A 的行列式 rank_A = rank(A )%矩阵A的秩 inv_B = inv(B )%矩阵 B 的逆矩阵 V,D = eig(B)%矩阵B的特征值矩阵V与特征向量构成的矩阵D X = A/B%A/B = A*B-1,即 XB=A，求 X Y = BA%BA = B-1*A，即 BY=A，求 Y上机练习(一)：1练习数据和符号的输入方式，将前面

10、的命令在命令窗口中执行通过；2. 输入 A=7 1 5； 2 5 6； 3 1 5, B=1 1 1; 2 2 2; 3 3 3，在命令 窗口中执行下列表达式，掌握其含义：A(2, 3)A(:,2)A(3,:)A(:,1:2:3)A(:,3).*B(:,2)A(:,3)*B(2,:) A*BA.*BAA2A.A2B/AB./A3. 输入 C=1:2:20,贝IC (i)表示什么？其中 i=l,2,3,.,10；4. 查找已创建变量的信息，删除无用的变量；5. 欲通过系统做一平面图,请查找相关的命令与函数,获取函数的帮助信息。、编程2.1 无条件循环当需要无条件重复执行某些命令时，可以使用for

11、循环: for循环变量t=表达式1 :达式2 :表达式3语句体end说明：表达式 1 为循环初值，表达式 2 为步长，表达式 3 为循环终值；当表达式 2 省 略时则默认步长为1; for语句允许嵌套。例6：生成3X4阶的Hiltber矩阵。for 1=1 : 3for j=1 : 4H G, j) =1/ (i+j-1); endend女口：矩阵输入程序m=input(矩阵行数：m=);n= input(矩阵列数：n=);for i=1:mfor j=1:ndisp(输入第,num2str(i),行，第，,num2str(j),列兀素)A(i, j) = input ()end end2.2

12、 条件循环1) if-else-then 语句if-else-then语句的常使用三种形式为:I(1) if逻辑表达式:(3) if逻辑表达式 1语句体11语句体 1end；elseif逻辑表达式 21.11语句体 2(2) if逻辑表达式1；elseif逻辑表达式 3语句体1111 else1:else语句体2111语句体 nend1end2) while循环语句while 循环的一般使用形式为：while 表达式语句体end例 7 ：用二分法计算多项式方程x3 - 2x -5 = 0在0, 3内的一个根。解：a = 0；fa = -inf； b = 3；fb = inf； while b-

13、a eps*bx =(a+b)/2；fx = xA3-2*x-5；if sign(fx)= sign(fa) a =x； fa = fx；elseb = x； fb = fx；endendx运行结果为：x = 2.09455151481542332.3 分支结构若需要对不同的情形执行不同的操作，可用switch分支语句： switch 表达式(标量或字符串)case 值 1语句体1case 值 2语句体2otherwise语句体 nend说明：当表达式不是“case”所列值时，执行otherwise语句体。2.4 建立 M 文件将多个可执行的系统命令，用文本编辑器编辑后并存放在后缀为 .m 的

14、文件中，若在 MATLAB命令窗口中输入该m-文件的文件名（不跟后缀.m!）,即可依次执行该文件中的多 个命令。这个后缀为.m的文件，也称为Matlab的脚本文件（Script File）。注意：文件存放路径必须在Mat lab能搜索的范围内。2.5 建立函数文件对于一些特殊用户函数，系统提供了一个用于创建用户函数的命令function，以备用户 随时调用。1格式：function 输出变量列表=fun_name（输入变量列表） 用户自定义的函数体2.函数文件名为：fun_name，注意：保存时文件名与函数名最好相同； 3存储路径：最好在系统的搜索路径上。4.调用方法：输出参量=fun_nam

15、e （输入变量）例 8 ：计算s = n!，在文本编辑器中输入：function s=pp（n）;s=1;for i=1:ns=s*i;ends;在 MATLAB 命令窗口中输入： s=pp（5）结果为 s = 120上机练习（二）：1编写程序，计算l+3+5+7+（2n+l）的值（用input语句输入n值）。x 0 x 12.编写分段函数f （x） = 2-x 1 x 2的函数文件，存放于文件ff.m中，计算出 、0 其它f（-3），f（迁），fe）的值。三、矩阵及其运算3.1 矩阵的创建1. 加、减运算 运算符：“”和“”分别为加、减运算符。 运算规则：对应元素相加、减，即按线性代数中矩阵

16、的“十”、“一”运算进行。例3-1 在Matlab编辑器中建立m文件：LX0701.mA=1, 1, 1; 1, 2, 3; 1, 3, 6B=8, 1, 6; 3, 5, 7; 4, 9, 2AB=A+BA-B=A-B在Matlab命令窗口建入LX0701，则结果显示： A+B=9274 7 105 128AB=-70-5-2-3-4-3-642.乘法运算符： * 运算规则：按线性代数中矩阵乘法运算进行，即放在前面的矩阵的各行元素，分别与 放在后面的矩阵的各列元素对应相乘并相加。（ 1 ）两个矩阵相乘例 3-2 在 Mtalab 编辑器中建立 M 文件： LX0702.mX= 2 3 4 5

17、1 2 2 1；Y=011110001100Z=X*Y存盘在命令行中建入LX0702，回车后显示：Z=856333(2) 矩阵的数乘：数乘矩阵上例中： a=2*X则显示： a =4 6 8 10 2442(3) 向量的点乘(内积)：维数相同的两个向量的点乘。命令：dot向量点乘函数例： X=-1 0 2；Y=-2 -1 1;Z=dot(X, Y)则显示： Z =4还可用另一种算法： sum(X.*Y) ans=4(4) 向量叉乘 在数学上，两向量的叉乘是一个过两相交向量的交点且垂直于两向量所在平面的向量在 Matlab 中，用函数 cross 实现。命令 cross 向量叉乘函数例 3-3 计

18、算垂直于向量(1, 2, 3)和(4, 5, 6)的向量。在Mtalab编辑器中建立M文件：LX0703.ma=123;b=4 5 6; c=cross(a,b)结果显示：c= -36-3可得垂直于向量(1, 2, 3)和(4, 5, 6)的向量为(-3, 6, -3)(5) 混合积混合积由以上两函数实现：例 3-4 计算向量 a=(l, 2,、3b=(4, 5,和)c=(-3, 6,-的混合积 a (b c) 在Matlab编辑器中建立M文件：LX0704.ma=1 2 3; b=45 6; c=-36 -3;x=dot(a, cross(b, c)结果显示：x =54注意：先叉乘后点乘，顺

19、序不可颠倒。3. 矩阵的除法Mat lab提供了两种除法运算：左除()和右除(/)。一般情况下，x=ab是方程a*x =b 的解，而x=b/a是方程x*a=b的解例：a=1 2 3; 4 2 6; 7 4 9b=4; 1; 2;x=ab则显示：x=-1.50002.00000.5000如果a为非奇异矩阵，则ab和 b/a可通过a的逆矩阵与b阵得到： ab = inv(a)*b b/a = b*inv(a)4. 矩阵乘方运算符：运算规则：(1)当A为方阵，p为大于0的整数时，AP表示A的P次方，即A自乘P次；p 为小于0的整数时，AP表示A-1的P次方。dp11V 1其中V为A的特征向 dpnn

20、量，11(2)当A为方阵，p为非整数时，则A P V为特征值矩阵dnn5. 矩阵的转置运算符： 运算规则：与线性代数中矩阵的转置相同。6.矩阵的逆矩阵3、1的逆矩阵3丿/1 2例 3-5 求 A = 2 213 4方法一：在Matlab编辑器中建立M文件：LX07051.mA=1 2 3; 2 2 1; 3 4 3;inv(A)或 AA(-1)则结果显示为ans =1.00003.0000-2.0000-1.5000-3.00002.50001.0000-1.000023100 、2 1 01 0进行初等行变换4300 1丿1.0000f 1方法二：由增广矩阵B= 213在Matlab编辑器中

21、建立M文件：LX07052.mB=1, 2, 3, 1, 0, 0; 2, 2, 1, 0, 1, 0; 3, 4, 3, 0, 0, 1；C=rref(B)%化行最简形X=C(:, 4:6)在Matlab命令窗口建入LX07052,则显示结果如下:C =1.0000001.00003.0000-2.000001.00000-1.5000-3.00002.5000001.00001.00001.0000-1.0000X =1.00003.0000-2.0000-1.5000-3.00002.50001.00001.0000-1.0000这就是 A 的逆矩阵。7.方阵的行列式命令： det 计算

22、行列式的值例3-6计算上例中A的行列式的值在Matlab编辑器中建立M文件：LX0706.mA=1 2 3; 2 2 1; 3 4 3;D=det(A)则结果显示为D =23.2 符号矩阵的运算1.符号矩阵的四则运算Matlab 5.x 抛弃了在 4.2 版中为符号矩阵设计的复杂函数形式，把符号矩阵的四则运算 简化为与数值矩阵完全相同的运算方式，其运算符为：加( + ),减(一)、乘(X)、除(/、 )等或：符号矩阵的和(symadd),差(symsub),乘(symmul)。例 3-7 A = sym(l/x, l/(x + l);l/(x + 2), 1/(x + 3);B = sym(x

23、, 1; x + 2, 0)；C=B-AD=ab则显示：C=x-l/x l-l/(x+l) x+2-l/(x+2)-l/(x+3)D=-6*x-2*xA3-7*xA21/2*xA3+x+3/2*xA26+2*xA3+10*xA2+14*x-2*xA2-3/2*x-1/2*xA32. 其他基本运算符号矩阵的其他一些基本运算包括转置()、行列式(det)、逆(inv)、秩(rank)、幕 (a)和指数(exp和expm)等都与数值矩阵相同3. 符号矩阵的简化 符号工具箱中提供了符号矩阵因式分解、展开、合并、简化及通分等符号操作函数。( 1)因式分解命令：factor符号表达式因式分解函数格式： f

24、actor(s)说明：s为符号矩阵或符号表达式。常用于多项式的因式分解例 3-8 将 x 9-1 分解因式在 Matlab 命令窗口建入syms xfactor(xA9-1)则显示： ans =(x-1)*(xA2+x+1)*(x6+xA3+1)例 3-9 问入取何值时，齐次方程组(1 九)x 2x + 4x = 0123 2x + (3 九)x + x = 0123x + x + (1 九)x = 0123有非0解解：在Matlab编辑器中建立M文件：LX0709.msyms kA=1-k -2 4;2 3-k 1;1 1 1-k;D=det(A) factor(D) 其结果显示如下：D =

25、-6*k+5*kA2-kA3ans = -k*(k-2)*(-3+k) 从而得到：当 k=0、k=2 或 k=3 时，原方程组有非 0 解。(2) 符号矩阵的展开 命令 expand 符号表达式展开函数 格式： expand(s)说明：s为符号矩阵或表达式。常用在多项式的因式分解中，也常用于三角函数，指数函数和对数函数的展开中例 3-10 将(x+l)3、sin(x+y)展开在Matlab编辑器中建立M文件：LX0710.msyms x y p=expand(x+1)A3) q=expand(sin(x+y)则结果显示为 p =xA3+3*xA2+3*x+1q =sin(x)*cos(y)+c

26、os(x)*sin(y)(3) 同类式合并命令：Collect 合并系数函数格式： Collect(s,v) 将 s 中的变量 v 的同幂项系数合并。Collect(s)s 一矩阵或表达式，此命令对由命令findsym函数返回的默认变量进行同类项合并。(4) 符号简化命令：simple或simplify寻找符号矩阵或符号表达式的最简型格式： Simple(s) s 一 矩阵或表达式说明：Simple(s)将表达式s的长度化到最短。若还想让表达式更加精美，可使用函数 Pretty。格式： Pretty(s) 使表达式 s 更加精美例 3-11 计算行列式1111abcda2 b2 c2 d 2a

27、4 b4 c4 d 4的值。在 Matlab 编辑器中建立 M 文件： LX0711.msyms a b c dA=l 1 1 1;a b c d;aA2 bA2 cA2 dA2;aA4 bA4 cA4 dA4;d1=det(A)d2=simple(d1)%化简表达式 d1pretty(d2)%让表达式 d2 符合人们的书写习惯则显示结果如下：d1 = b*cA2*dA4-b*dA2*cA4-bA2*c*dA4+bA2*d*cA4+bA4*c*dA2-bA4*d*cA2-a*cA2*dA4+a*dA2*cA4+a*bA2*dA4-a*bA2*cA4-a*bA4*dA2+a*bA4*cA2+aA

28、2*c*dA4-aA2*d*cA4-aA2*b*dA4+aA2*b*cA4+aA2*bA4*d-aA2*bA4*c-aA4*c*dA2+aA4*d*cA2+aA4*b*dA2-aA4*b*cA2-aA4*bA2*d+aA4*bA2*cd2 =(-d+c)*(b-d)*(b-c)*(-d+a)*(a-c)*(a-b)*(a+c+d+b)(-d+c)(b-d)(b-c)(-d+a)(a-c)(a-b)(a+c+d+b)例 3-12 设 A = 21，1、，3丿求矩阵X,使满足：AXB = C在 Matlab 编辑器中建立 M 文件： LX0712.mA=1 2 3;2 2 1;3 4 3;B=2,

29、1;5 3; C=1 3;2 0;3 1;X=AC/B 则结果显示如下：X =-2.00001.000010.0000-4.0000-10.00004.0000例 3-13计算(cos(t)(sin(t)-sin(t) )5cos(t) 丿在 Matlab 编辑器中建立 M 文件： LX0713.m syms tA =cos(t) -sin(t); sin(t), cos(t);B=sympow(A, 5)%计算 A 的 5 次幂C=simple(B)%化简pretty(C)则显示结果如下B =cos(t)*(cos(t)*(cos(t)*(cos(t)人2-sin(t)人2)-2*sin(t

30、)人2*cos(t)-sin(t)*(sin(t)*(cos(t)人2-sin(t) 人2)+2*cos(t)人2*sin(t)-sin(t)*(sin(t)*(cos(t)*(cos(t)人2-sin(t)人2)-2*sin(t)人2*cos(t)+cos(t)*(sin(t) *(cos(t)A2-sin(t)A2)+2*cos(t)A2*sin(t), cos(t)*(cos(t)*(-2*cos(t)A2*sin(t)-sin(t)*(cos(t)A2-sin(t)A2)-sin(t)*(cos(t)*(cos(t)A2-sin(t)A2)-2*s in(t)A2*cos(t)-sin

31、(t)*(sin(t)*(-2*cos(t)A2*sin(t)-sin(t)*(cos(t)A2-sin(t)A2)+cos(t)*(cos(t)*(cos(t )A2-sin(t)A2)-2*sin(t)A2*cos(t)sin(t)*(cos(t)*(cos(t)*(cos(t)A2-sin(t)A2)-2*sin(t)A2*cos(t)-sin(t)*(sin(t)*(cos(t)A2-sin(t) A2)+2*cos(t)A2*sin(t)+cos(t)*(sin(t)*(cos(t)*(cos(t)A2-sin(t)A2)-2*sin(t)A2*cos(t)+cos(t)*(sin(

32、t )*(cos(t)A2-sin(t)A2)+2*cos(t)A2*sin(t),sin(t)*(cos(t)*(-2*cos(t)A2*sin(t)-sin(t)*(cos(t)A2-sin(t)A2)-sin(t)*(cos(t)*(cos(t)A2-sin(t)A2 )-2*sin(t)A2*cos(t)+cos(t)*(sin(t)*(-2*cos(t)A2*sin(t)-sin(t)*(cos(t)A2-sin(t)A2)+cos(t)*(cos(t)* (cos(t)A2-sin(t)A2)-2*sin(t)A2*cos(t)C = cos(5*t), -sin(5*t) sin

33、(5*t), cos(5*t)cos(5 t)-sin(5 t) sin(5 t)cos(5 t) 四、秩与线性相关性4.1 矩阵和向量组的秩以及向量组的线性相关性。矩阵A的秩是矩阵A中最高阶非零子式的阶数；向量组的秩通常由该向量组构成的矩 阵来计算。命令：rank格式：rank(A) A为矩阵式向量组构成的矩阵例 4-1求向量组勺=(1- 223)，a2 = (-2 4 1 3)，叫=(T 2 0 3)，a4 = (0 62 3)，a5 = (2- 634)的秩，并判断其线性相关性。在Matlab编辑器中建立M文件：LX0714.mA=1 -2 2 3;-2 4 -1 3;-1 2 0 3;

34、0 6 2 3;2 -6 3 4;B=rank(A)运行后结果如下：B =3 由于秩为3 向量个数，因此向量组线性相关。4.2 向量组的最大无关组矩阵的初等行变换有三条：1交换两行ri分r.(第i、第j两行交换)2第i行的K倍 kri3.第i行的K倍加到第j行上去r. + kri通过这三条变换可以将矩阵化成行最简形，从而找出列向量组的一个最大无关组， Matlab 将矩阵化成行最简形的命令是：命令：rref格式：rref(A)A为矩阵例 4-2 求向量组 a1=(1,-2,2,3)，a2=(-2,4,-1,3)，a3=(-1,2,0,3)，a4=(0,6,2,3)，a5=(2,-6,3,4)

35、的一个最大无关组。在Matlab编辑器中建立M文件：LX0715.ma1=1-223;a2=-24-13;a3=-1203;a4=0623;a5=2-634;A=a1a2a3 a4 a5format rat%以有理格式输出B=rref(A) %求A的行最简形运行后的结果为A =1-2-102-2426-62-102333334B =101/3016/9012/30-1/90001-1/300000从B中可以得到：向量a1 a2 a4为其中一个最大无关组五、线性方程的组的求解我们将线性方程的求解分为两类：一类是方程组求唯一解或求特解，另一类是方程组 求无穷解即通解。可以通过系数矩阵的秩来判断：若

36、系数矩阵的秩r=n （n为方程组中未知变量的个数），则有唯一解若系数矩阵的秩rvn,则可能有无穷解。线性方程组的无穷解 = 对应齐次方程组的通解+非齐次方程组的一个特解；其特解的求法属于解的第一类问题,通解部分属第二类问题。5.1 求线性方程组的唯一解或特解（第一类问题）这类问题的求法分为两类：一类主要用于解低阶稠密矩阵 大型稀疏矩阵 迭代法。利用矩阵除法求线性方程组的特解直接法；另一类是解5.1.1或一个解）方程解法AX=bX=Ab例 5-1求方程组5x + 6x12x + 5x + 6x123x + 5x + 6x234x + 5x3+ 6x45 x + 5x45=1=0=0的解=0=1解

37、：在Matlab编辑器中建立M文件：LX0716.m A=5100065100065100065100065;%求秩%求解B=1 0 0 0 1; R_A=rank(A) X=AB 运行后结果如下 R_A =52.2662 -1.72181.0571 -0.59400.3188这就是方程组的解。例 5-2 求方程组X1 + X2 - 3X3 - X4 = 1 3x 一 x 一 3x + 4x = 41234x + 5x - 9x - 8x = 01234的一个特解解：在Matlab编辑器中建立M文件：LX0717.mA=1 1 -3 -1;3 -1 -3 4;1 5 -9 -8;B=1 4 0

38、;X=ABX =00-0.53330.60005.1.2利用矩阵的LU、QR和cholesky分解求方程组的解1.LU 分解：LU分解又称Gauss消去分解，可把任意方阵分解为下三角矩阵的基本变换形式(行交 换)和上三角矩阵的乘积。即A=LU, L为下三角阵，U为上三角阵。则： A*X=b 变成 L*U*X=b.X=U(Lb)这样可以大提高运算速度。命令L， U=lu (A)例 5-3 求方程组4x + 2x 一 x = 2111x12+ 3x23=8的一个特解解142-A =3一12b = 2,10, 8丫(1130 J In D:Matlabpujunlx0720.m at line 4

39、X =1.0e+016 * -0.40531.48621.3511说明：结果中的警告是由于系数行列式为零产生的。可以通过A*X验证其正确性。2.Cholesky 分解若A为对称正定矩阵，则Cholesky分解可将矩阵A分解成上三角矩阵和其转置的乘 积，即：A = R* R 其中R为上三角阵。方程 A*X=b 变成 R * R*X = b所以X = R(Rb)命令 R=chol(A)3.QR分解对于任何长方矩阵A,都可以进行QR分解，其中Q为正交矩阵，R为上三角矩阵的 初等变换形式，即： A=QR方程 A*X=b 变形成 QRX=b 所以X=R(Qb)命令 Q, R=qr(A) 上例中Q, R=

40、qr(A)X=R(QB)说明： 这三种分解，在求解大型方程组时很有用。其优点是运算速度快、可以节省磁 盘空间、节省内存。5.2 求线性齐次方程组的通解空间的一组基(基础解系)。格式：z = null% z的列向量为方程组的正交规范基，满足Zx Z = Iz = n u (A,t)% z的列向量是方程AX=0的有理基例 5-4 求解方程组的通解：xi + 2x2 + 2x3 + x4 = 0 2x + x 一 2x 一 2x = 0 1234I X1 一 X2 - 4X3 - 3X4 = 0 解：在Matlab编辑器中建立M文件：LX0719.m A=1 2 2 1;2 1 -2 -2;1 -1

41、 -4 -3; format rat%指定有理式格式输出B=null(A,r) %求解空间的有理基运行后显示结果如下：B =25/3-2-4/310 01写出通解： syms k1 k2X=k1*B(:,1)+k2*B(:,2)%写出方程组的通解pretty(X)%让通解表达式更加精美运行后结果如下：X = 2*k1+5/3*k2 -2*k1-4/3*k2k1k2%下面是其简化形式2k1 + 5/3k2 -2k1 - 4/3k2k1k25.3 求非齐次线性方程组的通解非齐次线性方程组需要先判断方程组是否有解，若有解，再去求通解。因此步骤为：第一步：判断AX=b是否有解，若有解则进行第二步第二步

42、：求 AX=b 的一个特解 第三步：求 AX=0 的通解 第四步：AX=b的通解=AX=O的通解+AX=b的一个特解。 例 5-5 求解方程组X1 - 2X2 + 3X3 - X4 = 1 3x x + 5x 一 3x = 212342x + x + 2x 一 2x = 31234解：在Matlab编辑器中建立M文件：LX0720.m A=1 -2 3 -1;3 -1 5 -3;2 1 2 -2; b=1 2 3;B=A b;n=4; R_A=rank(A) R_B=rank(B) format rat ifR_A=R_B&R_A=n%判断有唯一解X=Ab elseifR_A=R_B&R_An

43、%判断有无穷解X=Ab%求特解C=null(A,r)%求 AX=0 的基础解系else X=equition no solve%判断无解end 运行后结果显示：R_A =2R_B =3equition no solve说明该方程组无解例 5-6求解方程组的通解：x + x 一 3x 一 x = 11234K + 5x2 - 9x3 - 徨=0解法3x 一 x 一 3x + 4x = 41234在Matlab编辑器中建立M文件：LX07211.mA=1 1 -3 -1;3 -1 -3 4;1 5 -9 -8; b=1 4 0;B=A b;n=4;R_A=rank(A)R_B=rank(B)for

44、mat ratif R_A=R_B&R_A=nX=Abelseif R_A=R_B&R_A In D:Matlabpujunlx0723.m at line 11X =00-8/153/5C =-3/47/4013/23/210所以原方程组的通解为X=k1(3/23/210丿+k2(-3/47/401丿(0 )0-8/153/5 丿解法二：在Matlab编辑器中建立M文件：LX07212.mA=1 1 -3 -1;3 -1 -3 4;1 5 -9 -8;b=1 4 0;B=A b;C=rref(B) %求增广矩阵的行最简形，可得最简同解方程组 运行后结果显示为：C =1001-3/23/45/

45、40000对应齐次方程组的基础解系为：/ 3/2、-3/4、3/2，J =7/4101 0 J1 1丿非齐次方程组的特解为：/ 5/4、1/400丿所以，原方程组的通解为：x=ki i+k2.+n *六、特征值与二次型工程技术中的一些问题，如振动问题和稳定性问题，常归结为求一个方阵的特征值和 特征向量。6.1 方阵的特征值与特征向量设A为n阶方阵，如果数入和n维列向量x使得关系式Ax =九 x成立，则称九为方阵A的特征值，非零向量x称为A对应于特征值入的特征向量。在Matlab中，用如下几种调用格式来求A的特征值和特征向量。1. d=eig(A)% d为矩阵A的特征值排成的向量2. V, D=

46、eig(A)% D为A的特征值对角阵，V的列向量为对应特征值的特征向量(且为单位向量)3. V, D=eig(A, nobalance)%当A中有小到和截断误差相当的元素时，用nobalance 选项，其作用是减少计算误差(-2 11、例6-1求矩阵A = 02 0的特征值和特征向量4 1 3 丿解：在Matlab编辑器中建立M文件：LX0722.mA=-2 1 1;0 2 0;-4 1 3;V,D=eig(A)运行后结果显示：V =0 0 0.9045-0.7071 -0.9701 0.3015D =-1 0 0020002即：特征值-1 对应特征向量(-0.7071 0 -0.7071)T

47、特征值 2 对应特征向量(-0.2425 0 -0.9701)t和(-0.3015 0.9045 -0.3015)t1 1 0、例6-2求矩阵A = -4 3 0的特征值和特征向量I 1 0 2 丿解：在Matlab编辑器中建立M文件：LX0723.m A=-1 1 0;-4 3 0;1 0 2;V,D=eig(A) 运行后结果显示为V =00.4082-0.408200.8165-0.81651.0000-0.40820.4082D =2000 1 0001说明：当特征值为1 (二重根)时，对应特征向量都是k (0.4082 0.8165 -0.4082)Tk为 任意常数。6.2 正交矩阵及

48、二次型A为n阶方阵，且满足：AA = E (即A = A-i)，则A为正交矩阵A为正交矩阵o A的各列(行)向量的长度为1,而且A的各列(行)向量两两正交。6.2.1 向量的长度(范数)命令 norm格式 norm(X)% 求 X 的范数6.2.2 求矩阵的正交矩阵命令：orth格式： orth(A)%将矩阵 A 正交规范化(4 0 0、例如：求A = 0 3 1的正交矩阵1 3 丿在 Matlab 命令窗口键入A=4 0 0; 0 3 1; 01 3;P=orth(A)Q=P*P 则显示结果为 P =1.0000 0 0 0.7071 0 0.7071 Q =1.0000 00-0.7071

49、0.70710 1.000000001.00006.2.3 矩阵的 schur 分解格式 U, T=schur(A)%U为正交矩阵，使得A = UTU和UU = ET = schur(A)%生成schur矩阵T。当A为实对称阵时，T为特征值对角阵4 0 0、例6-3设A = 0 3 1，求一个正交矩阵P,使P-1AP=A为对角阵1 3 丿解法 1：在 Matlab 编辑器中建立 M 文件： LX07241.m A=4 0 0;0 3 1;0 1 3;V,D=eig(A)运行后结果如下：V =1.0000 0 00 0.7071 0.70710 0.7071 -0.7071D =40004000

50、2这里，V就是所求的正交矩阵P，D就是对角矩阵A 解法 2：在 Matlab 编辑器中建立 M 文件： LX07242.m A=4 0 0;0 3 1;0 1 3;U,T=schur(A)运行后结果显示如下U =1.00000000.70710.707100.7071-0.7071T =400040002这里，U就是所求的正交矩阵P, T就是对角矩阵A说明：对于实对称矩阵，用 eig 和 schur 分解效果一样例6-4求一个正交变换X=PY，把二次型f = 2x x + 2x x 一 2x x 一 2x x + 2x x + 2x x1 2 1 3 化成标准形。1 4232434解：先写出二

51、次型的实对称矩阵1011-A=i0 一 1ii一 1 0i(-1110 J在 Matlab 编辑器中建立 M 文件：LX0725.mA=0 1 1 -1;1 0 -1 1;1 -1 0 1;-1 1 1 0;P,D=schur(A) syms y1 y2 y3 y4 y=y1;y2;y3;y4；X=vpa(P,2)*y%vpa 表示可变精度计算，这里取 2 位精度f=y1 y2 y3 y4*D*y 运行后结果显示如下：P =780/989780/36911/2-390/1351780/3691780/989-1/2390/1351780/1351-780/1351-1/2390/1351001/21170/1351D =1000010000-300001X = .21*y1+.79*y2-.50*y3+.29*y4 .56*y1-.56*y2-.50*y3+.29*y4 .50*y3+.85*y4y1A2+y2A2-3*y3A2+y4A22 y+2 y 32y=f即