专题一：MATLAB 求解方程

单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,MATLAB,求解方程,制作,:,陈学明,基本概念,稀疏矩阵,(sparse matrix),：矩阵中非零元素的个数远远小于矩阵元素的总数），并且非零元素的分布没有规律。与之相区别的是，如果非零元素的分布存在规律（如上三角矩阵、下三角矩阵、对称矩阵），则称该矩阵为特殊矩阵。,稠密矩阵,:,非,0,元素占所有元素比例较大的矩阵。,若,n,阶矩阵,A,的行列式不为零，即,|A|0,，则称,A,为非奇异矩阵，否则称,A,为奇异矩阵。,把矩阵,A,的行换成相应的列，得到的新矩阵称为,A,的转置矩阵，记作,A,。,AA=E,（,E,为单位矩阵，）或,AA=E,，则,n,阶实矩阵,A,称为正交矩阵。,基本概念,求矩阵,A,的秩,rank(A,),求矩阵,A,的迹,trace(A,),求矩阵,A,的行列式,det(A,),求矩阵V,的1范数,norm(V,1),求矩阵V 的2范数,norm(V)或norm(V,2),求矩阵V,的范数,norm(V,inf,),魔方矩阵,魔方矩阵是每行、每列及两条对角线上的元素和都相等的矩阵。对于,n,阶魔方阵，其元素由,1,2,3,n,2,共,n,2,个整数组成,.,magic(n,),：生成一个,n,阶魔方阵,各行各列及两条 对角线和为,(1+2+3+.+n,2,)/n,范得蒙矩阵,范得蒙,(,Vandermonde,),矩阵最后一列全为,1,，倒数第二列为一个指定的向量，其他各列是其后列与倒数第二列的点乘积。,vander(V,),生成以向量,V,为基础向量的范得蒙矩阵。,例,A=vander(1;2;3;5),A=,1 1 1 1,8 4 2 1,27 9 3 1,125 25 5 1,希尔伯特矩阵,Hilbert,矩阵的每个元素,hij,=1/(i+j-1),hilb(n,),生成,n,阶希尔伯特矩阵,invhilb(n,),求,n,阶的希尔伯特矩阵的逆矩阵,注意,1,：高阶,Hilbert,矩阵一般为病态矩阵，所以直接求逆可能出现错误结论，故应该采用,invhilb(n,),注意,2,：由于,Hilbert,矩阵本身接近奇异，所以建议处理该矩阵时建议尽量采用符号运算工具箱，若采用数值解时应该验证结果的正确性,托普利兹矩阵,(,toeplitz,),toeplitz,矩阵除第一行第一列外，其他每个元素都与左上角的元素相同。,toeplitz(x,y,),生成一个以,x,为第一列，,y,为第一行的托普利兹矩阵。这里,x,y,均为向量，两者不必等长。,toeplitz(x,),用向量,x,生成一个对称的托普利兹矩阵。,例,T=toeplitz(1:5,1:7),T=,1 2 3 4 5 6 7,2 1 2 3 4 5 6,3 2 1 2 3 4 5,4 3 2 1 2 3 4,5 4 3 2 1 2 3,帕斯卡矩阵,二次项,(,x+y),n,展开后的系数随,n,的增大组成一个三角形表，称为杨辉三角形。由杨辉三角形表组成的矩阵称为帕斯卡,(Pascal),矩阵。,pascal(n,),生成一个,n,阶帕斯卡矩阵。,矩阵分解,矩阵分解通过将复杂矩阵表示成形式简单或具有良好数学性质（统称为简单矩阵）的组合，以便于理论分析或数值计算。通常矩阵分解将复杂矩阵分解为几个简单矩阵的乘积。求解线性方程组不可避免地要用到矩阵分解的概念。,MATLAB,中，线性方程组的求解主要用到三种基本的矩阵分解，即对称正定矩阵的,cholesky,分解、一般方程的,gaussian,消去法和矩阵的正交分解。这三种分解由函数,chol,、,lu,和,qr,完成。,矩阵的逆,求方阵,A,的逆可用,inv(A,),说明,1:,对于,Hilbert,求逆时，不建议用,inv,可用,invhilb,直接产生逆矩阵,说明,2,：符号工具箱中也对符号矩阵定义了,inv(),函数,即使对更高阶的奇异矩阵也可以精确的求解出逆矩阵,说明,3,：对于奇异矩阵用数值方法的,inv(),函数，会给出错误的结果，但符号工具箱中的,inv(),会给出正确的结论,例,A=16,2,3,13;5,11,10,8;9,7,6,12;4,14,15,1;,det(A,),B=,inv(A,),A=16,2,3,13;5,11,10,8;9,7,6,12;4,14,15,1;,A=sym(A),inv(A),矩阵的伪逆,pinv(A,),若,A,不是一个方阵，或,A,是一个非满秩的方阵时，矩阵,A,没有逆矩阵，但可以找到一个与,A,的转置矩阵,A,同型的矩阵,B,，使得：,ABA=A,BAB=B,此时称矩阵,B,为矩阵,A,的伪逆。,例 求矩阵,A,的伪逆,A=0,0,0;0,1,0;0,0,1;,pinv(A,),矩阵分解,矩阵分解通过将复杂矩阵表示成形式简单或具有良好数学性质（统称为简单矩阵）的组合，以便于理论分析或数值计算。通常矩阵分解将复杂矩阵分解为几个简单矩阵的乘积。求解线性方程组不可避免地要用到矩阵分解的概念。,MATLAB,中，线性方程组的求解主要用到三种基本的矩阵分解，即对称正定矩阵的,cholesky,分解、一般方程的,gaussian,消去法和矩阵的正交分解。这三种分解由函数,chol,、,lu,和,qr,完成。,矩阵分解之,LU,分解,矩阵的三角分解又称,LU,分解，它的目的是将一个矩阵分解成一个下三角矩阵,L,和一个上三角矩阵,U,的乘积，即,A=LU,。,Matlab,使用函数,lu,实现,LU,分解，其格式为：,L,U=,lu(A,),其中,U,为上三角阵，,L,为下三角阵或其变换形式，满足,LU=A,。,L,U,P=,lu(A,),U,为上三角阵，,L,为下三角阵，,P,为单位矩阵的行变换矩阵，满足,LU=PA,。,例,:,A=1 2 3;4 5 6;7 8 9;L,U=lu(A),L,U,P=lu(A),矩阵分解之,Cholesky分解,对于正定矩阵,A,，则存在一个实的非奇异上三角阵,R,，满足,R*R=A,，称为,Cholesky,分解。,Matlab,使用函数,chol,实现,Cholesky,分解，其格式为：,R=,chol(A,),若,A,非正定，则产生错误信息。,R,p,=,chol(A,),不产生任何错误信息，若,A,为正定阵，则,p=0,，,R,与上相同；若,A,非正定，则,p,为正整数，,R,是有序的上三角阵。,例：,A=pascal(4);,R,p,=,chol(A,),矩阵分解之,QR分解,将矩阵,A,分解成一个正交矩阵,Q,与一个上三角矩阵,R,的乘积,A=QR,，称为,QR,分解。,Matlab,使用函数,qr,实现,QR,分解，其格式为：,Q,R=,qr(A,),Q,R,E=,qr(A,),求得正交矩阵,Q,和上三角阵,R,，,E,为单位矩阵的变换形式，,R,的对角线元素按大小降序排列，满足,AE=QR,。,Q,R=qr(A,0),例：,A=1 2 3;4 5 6;7 8 9;10 11 12;Q,R=,qr(A,),代数方程的求解,一般的代数方程包括：,线性方程,非线性方程,超越方程,超越方程一般没有解析解，而只有数值解或近似解，只有特殊的超越方程才可以求出解析解来。,matlab,是获得数值解的一个最强大的工具。常用的命令有,fsolve,fzero,等，但超越方程的解很难有精确的表达式，因此在,matlab,中常用,eval,（）函数得到近似数值解，再用,vpa,（）函数控制精度。,求解线性方程,求解线性方程,对于,Ax=b,，,A,为,m*n,矩阵。,1,、,m=n,（恰定方程组）,2,、,mn,（超定方程），多用于曲线拟合,。,线性方程的组的求解,若系数矩阵的秩,r=n(n,为方程组中未知变量的个数,),，则有唯一解；,若系数矩阵的秩,r0%b,非零，,非齐次方程组,if,rank(A,)=,rank(A,b,)%,方程组相容（有解）,if,rank(A,)=n%,有唯一解,x=Ab;,else,%,方程组有无穷多个解,基础解系,disp,(,原方程组有有无穷个解，其齐次方程组的基础,解系为,y,，,特解为,x);,y=,null(A,r,);,x=Ab;,end,else,%,方程组不相容（无解），给出最小二乘解,disp,(,方程组的最小二乘法解是：,);,x=Ab;,end,else,%,齐次方程组,if,rank(A,)=n%,列满秩,x=zero(n,1)%0,解,else,%,非,0,解，无穷多个解,disp,(,方程组有无穷个解，基础解系为,x);,x=,null(A,r,);,end,end,return,如在,MATLAB,命令窗口，输入命令,A=2,2,-1,1;4,3,-1,2;8,5,-3,4;3,3,-2,2;b=4,6,12,6;,x,y=line_solution(A,b),及：,A=2,7,3,1;3,5,2,2;9,4,1,7;b=6,4,2;,x,y=line_solution(A,b),分别显示其求解结果。,迭代法求解线性方程组,迭代法适用于解大型稀疏方程组,(,万阶以上的方程组,系数矩阵中零元素占很大比例,而非零元按某种模式分布,),背景,:,电路分析、边值问题的数值解和数学物理方程,常用迭代法：,Jacobi,迭代法,Guass,-Seidel,迭代法,SOR,迭代法,编写实现,Jacobo,迭代法的函数,jacobi.m,如下,:,function s=jacobi(a,b,x0,err),%,jacobi,迭代法求解线性方程组，,a,为系数矩阵，,b,为,%ax=b,右边的,%,矩阵,b,x0,为迭代初值，,err,为迭代误差,if,nargin,=3%,nargin,判断输入变量个数的函数,err=1.0e-6;,elseif,nargin,=err,n=n+1,x0=s;,s=B*x0+f;,s,end,n,A=9-1-1;-1 10-1;-1-1 15;,b=7;8;13;,x0=0;0;0;,err=0.00005;,s=jacobi(A,b,x0,err),结果：,n,1,0.7778 0.8000 0.8667,2 0.9630 0.9644 0.9719,3 0.9929 0.9935 0.9952,4 0.9987 0.9988 0.9991,5 0.9998 0.9998 0.9998,6 1.0000 1.0000 1.0000,7 1.0000 1.0000 1.0000,例,编写实现,Seidel,迭代法的函数,seidel.m,如下,:,function s=seidel(a,b,x0,err),%,seidel,迭代法求解线性方程组，,a,为系数矩阵，,b,为,%ax=b,右边的,%,矩阵,b,x0,为迭代初值，,err,为迭代误差,if,nargin,=3,err=1.0e-6;,elseif,nargin,=err,n=n+1,x0=s;,s=B*x0+f;,s,end,n,A=9-1-1;-1 10-1;-1-1 15;,b=7;8;13;,x0=0;0;0;,err=0.00005;,s=seidel(A,b,x0,err),结果：,n,1 0.7778 0.8778 0.9770,2 0.9839 0.9961 0.9987,3 0.9994 0.9998 0.9999,4 1.0000 1.0000 1.0000,5 1.0000 1.0000 1.0000,例,发现：,seidel,收敛速度比,jacobi,快,