专题一MATLAB求解方程

1、关于专题一MATLAB求解方程第一张，PPT共四十一页，创作于2022年6月基本概念稀疏矩阵(sparse matrix)：矩阵中非零元素的个数远远小于矩阵元素的总数），并且非零元素的分布没有规律。与之相区别的是，如果非零元素的分布存在规律（如上三角矩阵、下三角矩阵、对称矩阵），则称该矩阵为特殊矩阵。稠密矩阵:非0元素占所有元素比例较大的矩阵。若n阶矩阵A的行列式不为零，即 |A|0，则称A为非奇异矩阵，否则称A为奇异矩阵。 把矩阵A的行换成相应的列，得到的新矩阵称为A的转置矩阵，记作A。AA=E（E为单位矩阵，）或AA=E，则n阶实矩阵 A称为正交矩阵 。第二张，PPT共四十一页，创作于20

2、22年6月基本概念求矩阵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) 第三张，PPT共四十一页，创作于2022年6月魔方矩阵魔方矩阵是每行、每列及两条对角线上的元素和都相等的矩阵。对于n阶魔方阵，其元素由1,2,3,n2共n2 个整数组成.magic(n) ：生成一个n阶魔方阵,各行各列及两条 对角线和为(1+2+3+.+n2 )/n第四张，PPT共四十一页，创作于2022年6月范得蒙矩阵范得蒙(Vandermonde)矩阵

3、最后一列全为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第五张，PPT共四十一页，创作于2022年6月希尔伯特矩阵Hilbert矩阵的每个元素hij=1/(i+j-1) hilb(n) 生成n阶希尔伯特矩阵 invhilb(n) 求n阶的希尔伯特矩阵的逆矩阵注意1：高阶Hilbert矩阵一般为病态矩阵，所以直接求逆可能出现错误结论，故应该采用invhilb(n) 注意2：由于Hilbert矩阵

4、本身接近奇异，所以建议处理该矩阵时建议尽量采用符号运算工具箱，若采用数值解时应该验证结果的正确性第六张，PPT共四十一页，创作于2022年6月托普利兹矩阵 (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

5、 1 2 3第七张，PPT共四十一页，创作于2022年6月帕斯卡矩阵二次项(x+y)n展开后的系数随n的增大组成一个三角形表，称为杨辉三角形。由杨辉三角形表组成的矩阵称为帕斯卡(Pascal)矩阵。 pascal(n) 生成一个n阶帕斯卡矩阵。第八张，PPT共四十一页，创作于2022年6月矩阵分解矩阵分解通过将复杂矩阵表示成形式简单或具有良好数学性质（统称为简单矩阵）的组合，以便于理论分析或数值计算。通常矩阵分解将复杂矩阵分解为几个简单矩阵的乘积。求解线性方程组不可避免地要用到矩阵分解的概念。MATLAB中，线性方程组的求解主要用到三种基本的矩阵分解，即对称正定矩阵的cholesky分解、一般

6、方程的gaussian消去法和矩阵的正交分解。这三种分解由函数chol、lu和qr完成。第九张，PPT共四十一页，创作于2022年6月矩阵的逆求方阵A的逆可用 inv(A)说明1:对于Hilbert求逆时，不建议用inv,可用 invhilb直接产生逆矩阵说明2：符号工具箱中也对符号矩阵定义了inv( )函数,即使对更高阶的奇异矩阵也可以精确的求解出逆矩阵说明3：对于奇异矩阵用数值方法的inv( )函数，会给出错误的结果，但符号工具箱中的inv( )会给出正确的结论第十张，PPT共四十一页，创作于2022年6月例 A=16,2,3,13;5,11,10,8; 9,7,6,12;4,14,15,

7、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)第十一张，PPT共四十一页，创作于2022年6月矩阵的伪逆 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)第十二张，PPT共四十一页，创作于2022年6月矩阵分解矩阵分解通过将复杂矩阵表示成形式简单或具有良好数学性质（统称为简单矩阵）的组合，

8、以便于理论分析或数值计算。通常矩阵分解将复杂矩阵分解为几个简单矩阵的乘积。求解线性方程组不可避免地要用到矩阵分解的概念。MATLAB中，线性方程组的求解主要用到三种基本的矩阵分解，即对称正定矩阵的cholesky分解、一般方程的gaussian消去法和矩阵的正交分解。这三种分解由函数chol、lu和qr完成。第十三张，PPT共四十一页，创作于2022年6月矩阵分解之LU分解矩阵的三角分解又称LU分解，它的目的是将一个矩阵分解成一个下三角矩阵L和一个上三角矩阵U的乘积，即A=LU。Matlab使用函数lu实现LU分解，其格式为：L,U = lu(A)其中U为上三角阵，L为下三角阵或其变换形式，满

9、足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)第十四张，PPT共四十一页，创作于2022年6月矩阵分解之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是有序的上

10、三角阵。例：A=pascal(4) ;R,p=chol(A)第十五张，PPT共四十一页，创作于2022年6月矩阵分解之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)第十六张，PPT共四十一页，创作于2022年6月代数方程的求解一般的代数方程包括：线性

11、方程非线性方程超越方程超越方程一般没有解析解，而只有数值解或近似解，只有特殊的超越方程才可以求出解析解来。matlab是获得数值解的一个最强大的工具。常用的命令有fsolve, fzero 等，但超越方程的解很难有精确的表达式，因此在matlab中常用eval（）函数得到近似数值解，再用vpa（）函数控制精度。第十七张，PPT共四十一页，创作于2022年6月求解线性方程第十八张，PPT共四十一页，创作于2022年6月求解线性方程对于Ax=b，A为m*n矩阵。1、m=n（恰定方程组）2、mn（超定方程），多用于曲线拟合。第十九张，PPT共四十一页，创作于2022年6月线性方程的组的求解若系数矩阵

12、的秩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 第三十一张，PPT共四十一页，创作于2022年6月 else %方程组不相容（无解） ，给出最小二乘解 disp(方程组的最小二乘法解是：); x=Ab; endelse %齐次方程组 if rank(A)=n %列满

13、秩 x=zero(n,1) %0解 else %非0解，无穷多个解 disp(方程组有无穷个解，基础解系为x); x=null(A,r); end endreturn第三十二张，PPT共四十一页，创作于2022年6月如在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)分别显示其求解结果。第三十三张，PPT共四十一页，创作于2022年6月迭代法求解线

14、性方程组迭代法适用于解大型稀疏方程组(万阶以上的方程组,系数矩阵中零元素占很大比例,而非零元按某种模式分布)背景: 电路分析、边值问题的数值解和数学物理方程常用迭代法：Jacobi迭代法Guass-Seidel迭代法SOR迭代法第三十四张，PPT共四十一页，创作于2022年6月编写实现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

15、nargin=err n=n+1 x0=s; s=B*x0+f; sendn第三十六张，PPT共四十一页，创作于2022年6月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例第三十七张，PPT共四十一页，创作于2022年6月编写实现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第三十九张，PPT共四十一页，创作于2022年6月A=9 -1 -1;-1 10 -1;-1 -1 15