专题一MATLAB求解方程

关于专题一MATLAB求解方程第1页，讲稿共41页，2023年5月2日，星期三基本概念稀疏矩阵(sparsematrix)：矩阵中非零元素的个数远远小于矩阵元素的总数），并且非零元素的分布没有规律。与之相区别的是，如果非零元素的分布存在规律（如上三角矩阵、下三角矩阵、对称矩阵），则称该矩阵为特殊矩阵。稠密矩阵:非0元素占所有元素比例较大的矩阵。若n阶矩阵A的行列式不为零，即|A|≠0，则称A为非奇异矩阵，否则称A为奇异矩阵。 把矩阵A的行换成相应的列，得到的新矩阵称为A的转置矩阵，记作A‘。AA‘=E（E为单位矩阵，）或A’A=E，则n阶实矩阵A称为正交矩阵。第2页，讲稿共41页，2023年5月2日，星期三基本概念求矩阵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页，2023年5月2日，星期三魔方矩阵魔方矩阵是每行、每列及两条对角线上的元素和都相等的矩阵。对于n阶魔方阵，其元素由1,2,3,…,n2共n2

个整数组成.magic(n)：生成一个n阶魔方阵,各行各列及两条对角线和为(1+2+3+...+n2)/n第4页，讲稿共41页，2023年5月2日，星期三范得蒙矩阵范得蒙(Vandermonde)矩阵最后一列全为1，倒数第二列为一个指定的向量，其他各列是其后列与倒数第二列的点乘积。

vander(V)生成以向量V为基础向量的范得蒙矩阵。例>>A=vander([1;2;3;5])A=11118421279311252551第5页，讲稿共41页，2023年5月2日，星期三希尔伯特矩阵Hilbert矩阵的每个元素hij=1/(i+j-1)hilb(n)生成n阶希尔伯特矩阵

invhilb(n)求n阶的希尔伯特矩阵的逆矩阵注意1：高阶Hilbert矩阵一般为病态矩阵，所以直接求逆可能出现错误结论，故应该采用invhilb(n)注意2：由于Hilbert矩阵本身接近奇异，所以建议处理该矩阵时建议尽量采用符号运算工具箱，若采用数值解时应该验证结果的正确性第6页，讲稿共41页，2023年5月2日，星期三托普利兹矩阵(toeplitz)toeplitz矩阵除第一行第一列外，其他每个元素都与左上角的元素相同。toeplitz(x,y)生成一个以x为第一列，y为第一行的托普利兹矩阵。这里x,y均为向量，两者不必等长。toeplitz(x)用向量x生成一个对称的托普利兹矩阵。例>>T=toeplitz([1:5],[1:7])T=12345672123456321234543212345432123第7页，讲稿共41页，2023年5月2日，星期三帕斯卡矩阵二次项(x+y)n展开后的系数随n的增大组成一个三角形表，称为杨辉三角形。由杨辉三角形表组成的矩阵称为帕斯卡(Pascal)矩阵。

pascal(n)生成一个n阶帕斯卡矩阵。第8页，讲稿共41页，2023年5月2日，星期三矩阵分解矩阵分解通过将复杂矩阵表示成形式简单或具有良好数学性质（统称为简单矩阵）的组合，以便于理论分析或数值计算。通常矩阵分解将复杂矩阵分解为几个简单矩阵的乘积。求解线性方程组不可避免地要用到矩阵分解的概念。MATLAB中，线性方程组的求解主要用到三种基本的矩阵分解，即对称正定矩阵的cholesky分解、一般方程的gaussian消去法和矩阵的正交分解。这三种分解由函数chol、lu和qr完成。第9页，讲稿共41页，2023年5月2日，星期三矩阵的逆求方阵A的逆可用inv(A)说明1:对于Hilbert求逆时，不建议用inv,可用invhilb直接产生逆矩阵说明2：符号工具箱中也对符号矩阵定义了inv()函数,即使对更高阶的奇异矩阵也可以精确的求解出逆矩阵说明3：对于奇异矩阵用数值方法的inv()函数，会给出错误的结果，但符号工具箱中的inv()会给出正确的结论第10页，讲稿共41页，2023年5月2日，星期三例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页，2023年5月2日，星期三矩阵的伪逆pinv(A)若A不是一个方阵，或A是一个非满秩的方阵时，矩阵A没有逆矩阵，但可以找到一个与A的转置矩阵A’同型的矩阵B，使得：

A·B·A=A,B·A·B=B此时称矩阵B为矩阵A的伪逆。例求矩阵A的伪逆

A=[0,0,0;0,1,0;0,0,1];pinv(A)第12页，讲稿共41页，2023年5月2日，星期三矩阵分解矩阵分解通过将复杂矩阵表示成形式简单或具有良好数学性质（统称为简单矩阵）的组合，以便于理论分析或数值计算。通常矩阵分解将复杂矩阵分解为几个简单矩阵的乘积。求解线性方程组不可避免地要用到矩阵分解的概念。MATLAB中，线性方程组的求解主要用到三种基本的矩阵分解，即对称正定矩阵的cholesky分解、一般方程的gaussian消去法和矩阵的正交分解。这三种分解由函数chol、lu和qr完成。第13页，讲稿共41页，2023年5月2日，星期三矩阵分解之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=[123;456;789];[L,U]=lu(A)[L,U,P]=lu(A)第14页，讲稿共41页，2023年5月2日，星期三矩阵分解之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页，2023年5月2日，星期三矩阵分解之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=[123;456;789;101112];[Q,R]=qr(A)第16页，讲稿共41页，2023年5月2日，星期三代数方程的求解一般的代数方程包括：线性方程非线性方程超越方程超越方程一般没有解析解，而只有数值解或近似解，只有特殊的超越方程才可以求出解析解来。matlab是获得数值解的一个最强大的工具。常用的命令有fsolve,fzero等，但超越方程的解很难有精确的表达式，因此在matlab中常用eval（）函数得到近似数值解，再用vpa（）函数控制精度。第17页，讲稿共41页，2023年5月2日，星期三求解线性方程→第18页，讲稿共41页，2023年5月2日，星期三求解线性方程对于Ax=b，A为m*n矩阵。1、m=n（恰定方程组）2、m<n（不定方程）3、m>n（超定方程），多用于曲线拟合。第19页，讲稿共41页，2023年5月2日，星期三线性方程的组的求解若系数矩阵的秩r=n(n为方程组中未知变量的个数)，则有唯一解；若系数矩阵的秩r<n，则可能有无穷解；求线性方程组的唯一解或特解（直接法）求线性方程组的通解线性方程组的无穷解=对应齐次方程组的通解+非齐次方程组的一个特解齐次线性方程组：常数项全部为零的线性方程组第20页，讲稿共41页，2023年5月2日，星期三直接法求线性方程组求逆矩阵的思想:x=inv(A)*b或x=A^(-1)*b左除法（原理上是运用高斯消元法求解,但Matlab在实际执行过程中是通过LU分解法进行的）：x=A\b符号矩阵法（此法最接近精确值，但运算速度慢）x=sym(A)\sym(b)化为行最简形：C=[A,b]R=rref(C)则R的最后一列元素就是所求之解。如果无解，那么返回的是一个单位矩阵。第21页，讲稿共41页，2023年5月2日，星期三例1：A=[-0.040.040.12;0.56-1.560.32;-0.241.24-0.28];b=[3;1;0];x11=inv(A)*b%法1：求逆思想x12=A^(-1)*b%法1：求逆思想x3=A\b%法2：左除法x4=sym(A)\sym(b)%法3：符号矩阵法C=[A,b];rref(C)%法4：化为行最简形第22页，讲稿共41页，2023年5月2日，星期三例:2：求下列方程组的解。A=[5600015600015600015600015];B=[10001]';R_A=rank(A) %求秩X=A\B %求解第23页，讲稿共41页，2023年5月2日，星期三例3：求下列方程组的解。A=[11-3-1;3-1-34;15-9-8];B=[140]'; R_A=rank(A) %求秩X=A\B %由于系数矩阵不满秩,该解法可能存在误差X=[00-0.53330.6000]‘%一个特解近似值A=[11-3-1;3-1-34;15-9-8];B=[140]';C=[A,B];%构成增广矩阵R=rref(C)若用rref求解，则比较精确：第24页，讲稿共41页，2023年5月2日，星期三LU分解求方程组的解LU分解又称Gauss消去分解，可把任意方阵分解为下三角矩阵和上三角矩阵的乘积。即

A=LU，L为下三角阵，U为上三角阵。则：

A*X=bL*U*X=b所以

X=U\(L\b)这样可以大大提高运算速度。若采用[l，u，p]=lu(A)，则Ax=b的解为x=u(l\(p*b))第25页，讲稿共41页，2023年5月2日，星期三利用LU分解求方程组的解例4:求下列方程组的一个特解。

解：A=[42-1;3-12;1130];B=[2108]';[L,U]=lu(A)X=U\(L\B)第26页，讲稿共41页，2023年5月2日，星期三

例5：利用LU分解法求解方程组A=[2426;49615;26918;6151840]%先录入系数矩阵b=[9232247]’[L,U]=lu(A)%对A矩阵做LU分解y=L\b%求解方程组Ly=b

x=U\y%求解方程组Ux=y得到方程组的最终解故方程组的最终解为:

x=(0.5000,2.0000,3.0000,-1.0000)’第27页，讲稿共41页，2023年5月2日，星期三或A=[2426;49615;26918;6151840]%先录入系数矩阵b=[9232247]’[L,U,P]=lu(A)%对A矩阵做LU分解y=L\P*b%求解方程组Ly=Pb

x=U\y%求解方程组Ux=y得到方程组的最终解第28页，讲稿共41页，2023年5月2日，星期三Cholesky分解求方程组的解若A为对称正定矩阵，则Cholesky分解可将矩阵A分解成上三角矩阵和其转置的乘积，即：

A=R'*R其中R为上三角阵。方程A*X=b变成R‘*R*X=b所以 X=R\(R'\b)第29页，讲稿共41页，2023年5月2日，星期三QR分解求方程组的解对于任何长方矩阵A，都可以进行QR分解，其中Q为正交矩阵，R为上三角矩阵的初等变换形式，即：

A=QR方程A*X=b变形成

QRX=b所以

X=R\(Q\b)若采用[Q，R，E]=qr(X)，则Ax=b的解为 x=E(R\(Q\b)) 第30页，讲稿共41页，2023年5月2日，星期三解线性方程组的一般函数文件如下：function[x,y]=line_solution(A,b)[m,n]=size(A);y=[];

ifnorm(b)>0%b非零，非齐次方程组

ifrank(A)==rank([A,b])%方程组相容（有解）

ifrank(A)==n%有唯一解

x=A\b;

else%方程组有无穷多个解,基础解系

disp('原方程组有有无穷个解，其齐次方程组的基础解系为y，特解为x');y=null(A,'r');x=A\b;

end

第31页，讲稿共41页，2023年5月2日，星期三

else%方程组不相容（无解），给出最小二乘解

disp('方程组的最小二乘法解是：');x=A\b;

endelse%齐次方程组

ifrank(A)==n%列满秩

x=zero(n,1)%0解

else%非0解，无穷多个解

disp('方程组有无穷个解，基础解系为x');x=null(A,'r');

end

endreturn第32页，讲稿共41页，2023年5月2日，星期三如在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页，2023年5月2日，星期三迭代法求解线性方程组迭代法适用于解大型稀疏方程组(万阶以上的方程组,系数矩阵中零元素占很大比例,而非零元按某种模式分布)背景:电路分析、边值问题的数值解和数学物理方程常用迭代法：Jacobi迭代法Guass-Seidel迭代法SOR迭代法第34页，讲稿共41页，2023年5月2日，星期三编写实现Jacobo迭代法的函数jacobi.m如下:functions=jacobi(a,b,x0,err)%jacobi迭代法求解线性方程组，a为系数矩阵，b为%ax=b右边的%矩阵b,x0为迭代初值，err为迭代误差ifnargin==3 %nargin判断输入变量个数的函数

err=1.0e-6;elseifnargin<3errorreturnendD=diag(diag(a)); %构造对角阵DD1=inv(D); %求对角阵D的逆矩阵L=tril(a,-1); %构造严格下三角阵U=triu(a,1); %构造严格上三角阵B=-D1*(L+U); %求出迭代矩阵第35页，讲稿共41页，2023年5月2日，星期三f=D1*b;s=B*x0+f;n=1;%n为迭代次数whilenorm(s-x0)>=errn=n+1x0=s;s=B*x0+f;s'endn第36页，讲稿共41页，2023年5月2日，星期三A=[9-1-1;-110-1;-1-115];b=[7;8;13];x0=[0;0;0];err=0.00005;s=jacobi(A,b,x0,err)结果：n

1

0.77780.80000.866720.96300.96440.971930.99290.99350.995240.99870.99880.999150.99980.99980.9998