《专题一：MATLAB 求解方程.ppt》由会员分享，可在线阅读，更多相关《专题一：MATLAB 求解方程.ppt（40页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、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的行列式
2、det(A)求矩阵V 的1范数 norm(V,1)求矩阵V 的2范数 norm(V)或norm(V,2)求矩阵V 的范数 norm(V,inf)魔方矩阵魔方矩阵是每行、每列及两条对角线上的元素和都相等的矩阵。对于n阶魔方阵，其元素由1,2,3,n2共n2 个整数组成.magic(n)：生成一个n阶魔方阵,各行各列及两条 对角线和为(1+2+3+.+n2)/n范得蒙矩阵范得蒙(Vandermonde)矩阵最后一列全为1，倒数第二列为一个指定的向量，其他各列是其后列与倒数第二列的点乘积。vander(V)生成以向量V为基础向量的范得蒙矩阵。例 A=vander(1;2;3;5)A=1 1 1 1
3、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为第一行的
4、托普利兹矩阵。这里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阶帕斯卡矩阵。矩阵分解矩阵分解通过将复杂矩阵表示成形式简单或具有良好数学性质（统称为简单矩阵）的组合，以便于理论分析或数值计算。通常矩阵分解将
5、复杂矩阵分解为几个简单矩阵的乘积。求解线性方程组不可避免地要用到矩阵分解的概念。MATLAB中，线性方程组的求解主要用到三种基本的矩阵分解，即对称正定矩阵的cholesky分解、一般方程的gaussian消去法和矩阵的正交分解。这三种分解由函数chol、lu和qr完成。矩阵的逆求方阵A的逆可用 inv(A)说明1:对于Hilbert求逆时，不建议用inv,可用 invhilb直接产生逆矩阵说明2：符号工具箱中也对符号矩阵定义了inv()函数,即使对更高阶的奇异矩阵也可以精确的求解出逆矩阵说明3：对于奇异矩阵用数值方法的inv()函数，会给出错误的结果，但符号工具箱中的inv()会给出正确的结论
6、例例 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)矩阵分解矩阵分解通过将复杂矩阵表示成形式简单或具有良好数学性质（统称为简单矩阵）的组合，以便于理论分析或数值计算
7、。通常矩阵分解将复杂矩阵分解为几个简单矩阵的乘积。求解线性方程组不可避免地要用到矩阵分解的概念。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为单位矩阵的
8、行变换矩阵，满足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，称为
9、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（）函
10、数得到近似数值解，再用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(原方程组有有无穷个解，其齐次方程组的基础原方程
11、组有有无穷个解，其齐次方程组的基础 解系为解系为y，特解为特解为x);y=null(A,r);x=Ab;end else%方程组不相容（无解）方程组不相容（无解），给出最小二乘解，给出最小二乘解 disp(方程组的最小二乘法解是：方程组的最小二乘法解是：);x=Ab;endelse%齐次方程组齐次方程组 if rank(A)=n%列满秩列满秩 x=zero(n,1)%0解解 else%非非0解，无穷多个解解，无穷多个解 disp(方程组有无穷个解，基础解系为方程组有无穷个解，基础解系为x);x=null(A,r);end endreturn如在如在MATLAB命令窗口，输入命令命令窗口，输入命
12、令 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)分别显示其求解结果。分别显示其求解结果。迭代法求解线性方程组迭代法适用于解大型稀疏方程组迭代法适用于解大型稀疏方程组(万阶以上的万阶以上的方程组方程组,系数矩阵中零元素占很大比例系数矩阵中零元素占很大比例,而非零而非零元按某种模式分布元按某种模式分布)背景背景:电路分析、边值问题的数值解和数学物电路分析、边值问题的数值解和数学物理方
13、程理方程常用迭代法：常用迭代法：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
14、+f;sendnA=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
15、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;sendnA=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快快