专题一matlab求解方程.ppt
关于专题一MATLAB求解方程现在学习的是第1页，共41页基本概念 稀疏矩阵(sparse matrix)：矩阵中非零元素的个数远远小于矩阵元素的总数），并且非零元素的分布没有规律。与之相区别的是，如果非零元素的分布存在规律（如上三角矩阵、下三角矩阵、对称矩阵），则称该矩阵为特殊矩阵。 稠密矩阵:非0元素占所有元素比例较大的矩阵。 若n阶矩阵A的行列式不为零，即 |A|0，则称A为非奇异矩阵，否则称A为奇异矩阵。 把矩阵A的行换成相应的列，得到的新矩阵称为A的转置矩阵，记作A。 AA=E（E为单位矩阵，）或AA=E，则n阶实矩阵 A称为正交矩阵 。现在学习的是第2页，共41页基本概念 求矩阵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) 现在学习的是第3页，共41页魔方矩阵 魔方矩阵是每行、每列及两条对角线上的元素和都相等的矩阵。对于n阶魔方阵，其元素由1,2,3,n2共n2 个整数组成. magic(n) ：生成一个n阶魔方阵,各行各列及两条 对角线和为(1+2+3+.+n2 )/n现在学习的是第4页，共41页范得蒙矩阵 范得蒙(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现在学习的是第5页，共41页希尔伯特矩阵 Hilbert矩阵的每个元素hij=1/(i+j-1) hilb(n) 生成n阶希尔伯特矩阵 invhilb(n) 求n阶的希尔伯特矩阵的逆矩阵 注意1：高阶Hilbert矩阵一般为病态矩阵，所以直接求逆可能出现错误结论，故应该采用invhilb(n) 注意2：由于Hilbert矩阵本身接近奇异，所以建议处理该矩阵时建议尽量采用符号运算工具箱，若采用数值解时应该验证结果的正确性现在学习的是第6页，共41页托普利兹矩阵 (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现在学习的是第7页，共41页帕斯卡矩阵 二次项(x+y)n展开后的系数随n的增大组成一个三角形表，称为杨辉三角形。由杨辉三角形表组成的矩阵称为帕斯卡(Pascal)矩阵。 pascal(n) 生成一个n阶帕斯卡矩阵。现在学习的是第8页，共41页矩阵分解 矩阵分解通过将复杂矩阵表示成形式简单或具有良好数学性质（统称为简单矩阵）的组合，以便于理论分析或数值计算。通常矩阵分解将复杂矩阵分解为几个简单矩阵的乘积。求解线性方程组不可避免地要用到矩阵分解的概念。 MATLAB中，线性方程组的求解主要用到三种基本的矩阵分解，即对称正定矩阵的cholesky分解、一般方程的gaussian消去法和矩阵的正交分解。这三种分解由函数chol、lu和qr完成。现在学习的是第9页，共41页矩阵的逆 求方阵A的逆可用 inv(A) 说明1:对于Hilbert求逆时，不建议用inv,可用 invhilb直接产生逆矩阵 说明2：符号工具箱中也对符号矩阵定义了inv( )函数,即使对更高阶的奇异矩阵也可以精确的求解出逆矩阵 说明3：对于奇异矩阵用数值方法的inv( )函数，会给出错误的结果，但符号工具箱中的inv( )会给出正确的结论现在学习的是第10页，共41页 例例 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)现在学习的是第11页，共41页矩阵的伪逆 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)现在学习的是第12页，共41页矩阵分解 矩阵分解通过将复杂矩阵表示成形式简单或具有良好数学性质（统称为简单矩阵）的组合，以便于理论分析或数值计算。通常矩阵分解将复杂矩阵分解为几个简单矩阵的乘积。求解线性方程组不可避免地要用到矩阵分解的概念。 MATLAB中，线性方程组的求解主要用到三种基本的矩阵分解，即对称正定矩阵的cholesky分解、一般方程的gaussian消去法和矩阵的正交分解。这三种分解由函数chol、lu和qr完成。现在学习的是第13页，共41页矩阵分解之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)现在学习的是第14页，共41页矩阵分解之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)现在学习的是第15页，共41页矩阵分解之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)现在学习的是第16页，共41页代数方程的求解 一般的代数方程包括： 线性方程 非线性方程 超越方程 超越方程一般没有解析解，而只有数值解或近似解，只有特殊的超越方程才可以求出解析解来。 matlab是获得数值解的一个最强大的工具。常用的命令有fsolve, fzero 等，但超越方程的解很难有精确的表达式，因此在matlab中常用eval（）函数得到近似数值解，再用vpa（）函数控制精度。现在学习的是第17页，共41页求解线性方程11 11221121 1222221 122mnnnnmmmnna xa xa xba xa xa xba xaxaxbL L L L L L L L L L L L L L L L LLLLAxb现在学习的是第18页，共41页求解线性方程 对于Ax=b，A为m*n矩阵。 1、m=n（恰定方程组）（恰定方程组） 2、mn（超定方程），多用于曲线拟合（超定方程），多用于曲线拟合。x直接法：理论，无舍入误差，有限步，精确解迭代法：格式，无穷序列解向量现在学习的是第19页，共41页线性方程的组的求解 若系数矩阵的秩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 现在学习的是第31页，共41页 else %方程组不相容（无解）方程组不相容（无解） ，给出最小二乘解，给出最小二乘解 disp(方程组的最小二乘法解是：方程组的最小二乘法解是：); x=Ab; endelse %齐次方程组齐次方程组 if rank(A)=n %列满秩列满秩 x=zero(n,1) %0解解 else %非非0解，无穷多个解解，无穷多个解 disp(方程组有无穷个解，基础解系为方程组有无穷个解，基础解系为x); x=null(A,r); end endreturn现在学习的是第32页，共41页如在如在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)分别显示其求解结果。分别显示其求解结果。现在学习的是第33页，共41页迭代法求解线性方程组 迭代法适用于解大型稀疏方程组迭代法适用于解大型稀疏方程组(万阶以上的万阶以上的方程组方程组,系数矩阵中零元素占很大比例系数矩阵中零元素占很大比例,而非零而非零元按某种模式分布元按某种模式分布) 背景背景: 电路分析、边值问题的数值解和数学物电路分析、边值问题的数值解和数学物理方程理方程 常用迭代法：常用迭代法： Jacobi迭代法迭代法 Guass-Seidel迭代法迭代法 SOR迭代法迭代法现在学习的是第34页，共41页编写实现编写实现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; sendn现在学习的是第36页，共41页123123123971081513xxxxxxxxx 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例例现在学习的是第37页，共41页编写实现编写实现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; sendn现在学习的是第39页，共41页123123123971081513xxxxxxxxx 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快快 现在学习的是第40页，共41页感谢大家观看现在学习的是第41页，共41页
