线性方程组求解的数值方法课件

Ch3线性方程组求解的数值方法NumericalSolutionsofSystemsofLinearEquationsCh3线性方程组求解的数值方法NumericalSolu1§3.0Introduction线性方程组工程问题（电子网络，船体放样等）自然科学（实验数据的最小二乘法曲线拟合）解非线性方程组解常/偏微分方程组（差分法，有限元法）§3.0Introduction线性方程组工程问题（电子网23.0Introduction线性方程组的系数矩阵：（1）低价稠密矩阵（阶数<150）；（2）大型稀疏矩阵（阶数高且零元素较多）。求解线性方程组的方法：（1）直接法:经过有限步算术运算，求得精确解（假设计算过程没有舍入误差）。如Gauss消去法，三角分解法。（2）间接法(迭代法)：通过迭代序列，逐步逼近方程组的解。如Jacobi迭代法，Gauss-Seidel迭代法。3.0Introduction线性方程组的系数矩阵：（1）33.0Introduction【最简单的情形】三角形线性方程组：简记为，A为上三角矩阵。若所有aii0，可用“回代”过程得到方程组的解。3.0Introduction【最简单的情形】三角形线性方4functionX=backsub(U,b)%解上三角方程组--回代过程%BacksubstitutioninGausselimination%Input--Uisanxnupper-trianglularmatrix%--bisanx1constantvctor%Output--XisthesolutionvectorofUX=b%Usage:X=backsub(U,b)%FindthedimensionofbandinitializeXn=length(b);X=zeros(n,1);X(n)=b(n)/U(n,n);forp=n-1:-1:1X(p)=(b(p)-U(p,p+1:n)*X(p+1:n))/U(p,p);endBacksub.m3.0IntroductionfunctionX=backsub(U,b)Backsub53.0Introduction【作业】（1）编写Matlab程序解下列上三角形方程组：（2）编写Matlab程序解下列下三角形方程组：3.0Introduction【作业】（1）编写Matla6§3.1Gauss消去法与LU分解法直接法的基本思想：将线性方程组化成与之等价的上三角形或下三角形，再用回代法求解。它的核心是矩阵分解。核心：矩阵分解。§3.1Gauss消去法与LU分解法直接法的基本思想：将线7Gauss消去法（Gaussianelimination）：（对方程组的三种变换）（1）交换两个方程的次序；（2）用一个非零常数乘一个方程；（3）将一个方程的非零倍数加到另一个方程上去。这等价与对增广矩阵进行三种行初等变换，将它化为行阶梯形。Gauss消去法（Gaussianelimination）8【例3.1】解下列线性方程组：解：用行初等变换将方程组的增广矩阵化为行阶梯形：解得

x=(0,-1,1)’.§3.1GaussianEliminationandLUDecomposition【例3.1】解下列线性方程组：解：用行初等变换将方程组的增广9【注】上述过程可用矩阵表示为：LU=PA其中§3.1GaussianEliminationandLUDecomposition【注】上述过程可用矩阵表示为：LU=PA其中§3.1Gau10【高斯消元法思路】(1)用消元法将A化为上三角阵

upper-triangularmatrix;(2)回代求解backwardsubstitution。§3.1GaussianEliminationandLUDecomposition【高斯消元法思路】(1)用消元法将A化为上三角阵upp11§3.1GaussianEliminationandLUDecomposition消元记Step1：设

，计算因子将增广矩阵/*augmentedmatrix*/第i行mi1

第1行，得到其中Stepk：设，计算因子且计算共进行？步n

1§3.1GaussianEliminationand12回代Whatif?Nouniquesolutionexists.Whatif?Thenwemustfindthesmallestintegerkiwith,andinterchangethek-throwwiththei-throw.Whatifwecan’tfindsuchk

?Nouniquesolutionexists.定理若A的所有顺序主子式

/*determinantofleadingprincipalsubmatrices*/

均不为0，则高斯消元无需换行即可进行到底，得到唯一解。注：事实上，只要A非奇异，即A1存在，则可通过逐次消元及行交换，将方程组化为三角形方程组，求出唯一解。§3.1GaussianEliminationandLUDecomposition回代Whatif?Noun13高斯消元法的矩阵表达----矩阵分解课本P40-42§3.1GaussianEliminationandLUDecomposition高斯消元法的矩阵表达----矩阵分解课本P40-42§3.1143.1.2矩阵的LU分解（LUDecomposition）若在Gauss消去法过程中，没有进行交换行的变换，即P=I,则分解结果为§3.1GaussianEliminationandLUDecomposition这样的分解式称为矩阵的LU分解。【Matlab】>>[L,U,P]=lu(A)其中【作业】P.67题1,23.1.2矩阵的LU分解（LUDecomposition153.1.3选主元消去法

/*PivotingStrategies*/例：单精度解方程组/*精确解为和*/8个8个用GaussianElimination计算：8个小主元/*Smallpivotelement*/可能导致计算失败。§3.1GaussianEliminationandLUDecomposition3.1.3选主元消去法/*PivotingStra16

全主元消去法

/*CompletePivoting*/每一步选绝对值最大的元素为主元素，保证。Stepk:

①选取②Ifik

k

then交换第k行与第ik

行;Ifjk

k

then交换第k列与第jk

列;③消元注：列交换改变了xi的顺序，须记录交换次序，解完后再换回来。列主元消去法

/*PartialPivoting,ormaximalcolumnpivoting*/省去换列的步骤，每次仅选一列中最大的元。§3.1GaussianEliminationandLUDecomposition全主元消去法/*CompletePivoting17例：注：列主元法没有全主元法稳定。例：注意：这两个方程组在数学上严格等价。标度化列主元消去法/*ScaledPartialPivoting*/对每一行计算。为省时间，si只在初始时计算一次。以后每一步考虑子列中最大的aik为主元。注：稳定性介于列主元法和全主元法之间。§3.1GaussianEliminationandLUDecomposition例：注：列主元法没有全主元法稳定。例：注意：这两个方程组在18实际应用中直接调用GaussElimination解3阶线性方程组的结果：结合全主元消去后的结果：§3.1GaussianEliminationandLUDecomposition实际应用中直接调用GaussElimination解3阶19§3.2Cholesky分解【Matlab】设A是对称正定矩阵，则A有以下形式的LU分解：A=PTP（3.2.1）其中P是一个上三角矩阵,上式称为A的Cholesky分解(Choleskiandecomposition).>>P=chol(A)定理【作业】P.68题4§3.2Cholesky分解【Matlab】设A是对称正定20§3.3向量范数与矩阵范数—误差的度量3.3.1

向量范数(vectornorms)常用向量范数：这些范数满足：一般范数：§3.3向量范数与矩阵范数—误差的度量3.3.1向量范21如，设向量x=(2,-4,-1)T,则§3.3向量范数与矩阵范数【Matlab】x=(1:4)/5N1=norm(x,1)N2=norm(x)N7=norm(x,7)Ninf=norm(x,inf)x=0.20000.40000.60000.8000N1=2N2=1.0954N7=0.8153Ninf=0.8000如，设向量x=(2,-4,-1)T,则§3.3向22定义

矩阵A的范数定义为§3.3向量范数与矩阵范数3.3.2

矩阵范数(matrixnorms)

命题：矩阵的常用范数：其中，λ1,…,λ

n是的特征值，称为的谱半径。定义矩阵A的范数定义为§3.3向量范数与矩阵范数3.3.23例

设矩阵§3.3向量范数与矩阵范数则可以证明，矩阵范数满足：例设矩阵§3.3向量范数与矩阵范数则可以证明，矩阵范数满24§3.4经典迭代法－Jacobi迭代法与Gauss-Seidel迭代法3.4.1．Jacobi迭代法设有n阶方程组(3.4.1)

§3.4经典迭代法－Jacobi迭代法与Gauss-Sei25若系数矩阵非奇异，且(i=1,2,…,n)，将方程组（4.1）改写成§3.4经典迭代法若系数矩阵非奇异，且(26然后写成迭代格式（3.4.2）（3.4.2）式也可以简单地写为（3.4.3）§3.4经典迭代法然后写成迭代格式（3.4.2）（3.4.2）式也可以简单地写27写成矩阵形式：A=LUDBJacobi迭代阵(3.4.6)§3.4经典迭代法写成矩阵形式：A=LUDBJacobi迭代阵(3.4.628…………只存一组向量即可。写成矩阵形式：BGauss-Seidel迭代阵3.4.2Gauss-Seidel迭代法(3.4.7)§3.4经典迭代法…………只存一组向量即可。写成矩阵形29定理3.1

对任意初始向量X(0)及常向量F，迭代格式(3.5.1)(3.5.1)推论：若迭代矩阵Ｂ的某种范数，则（3.5.1）收敛的充分必要条件是迭代矩阵B的谱半径(B)<1。确定的迭代法对任意初值X(0)均收敛于方程组X=BX+F§3.5迭代法的分析——迭代法的收敛性、收敛速度与数值稳定性3.5.1收敛性的唯一解x*。定理3.1对任意初始向量X(0)及常向量F，迭代格式(3.30定义1：如果矩阵的每一行中，不在主对角线上的所有元素绝对值之和小于主对角线上元素的绝对值，即

则称矩阵A(按行)严格对角占优。定理3.2：若A是严格对角占优，则A是非奇异的。§3.5迭代法的分析定理3：若线性方程组AX=b的系数矩阵A按行严格对角占优，则Jacobi迭代法和Gauss-Seidel迭代法对任意给定初值均收敛。定义1：如果矩阵的每一行中，不在主对角线上的313.5.2迭代法的误差估计与收敛速度定理3.4

设X*是方程组AX=b的同解方程X=BX+F的准确解，若迭代公式中迭代矩阵B的某种范数，(1)(2)则有§3.5迭代法的分析3.5.2迭代法的误差估计与收敛速度定理32迭代加速（3.6.2）§3.6超松驰法(SOR)及分块迭代法定理3.５

SOR方法收敛的必要条件是其中ω称为松弛因子．（3.6.1）迭代加速（3.6.2）§3.6超松驰法(SOR)及分块迭代33§3.7条件数、病态方程组与算法的稳定性3.7.1

病态方程组与数值稳定的算法定义对线性方程组Ax=b,如果解x关于问题（矩阵A和b）的“微小”变化（即舍入误差）不敏感，则称此方程组是“好的”，或称为“良态的”，否则，称之为“坏的”，或称为“病态的”。定义对求解线性方程组Ax=b的一个(迭代)算法，如果关于初始值的误差，不随计算逐步放大，则称算法是数值稳定的，否则，称之为是数值不稳定的。数值不稳定的算法是不可用的！§3.7条件数、病态方程组与算法的稳定性3.7.1病态方34例（病态方程组）考虑线性方程组§3.7条件数、病态方程组与算法的稳定性已知其精确解是u=(1,1,1,1)T.(1)对右边常数向量b作“微小”扰动，得方程(2)对系数矩阵A作“微小”扰动，得方程组它的解是v=(9.2,-12.6,4.5,-1.1)T.它的解是w=(-81,137,-34,22)T.注意：上述两种情况的解都是精确解！例（病态方程组）考虑线性方程组§3.7条件数、病态方程组与35一般地，对线性方程组Au=b（1）两方程相减得于是(2)由此得到设其解为v=u+u,即设(1)考虑对此线性方程组的第一类扰动方程组（2）式右端给出了解的相对误差界。§3.7条件数、病态方程组与算法的稳定性一般地，对线性方程组两方程相减得于是(2)由此得到设其解36两方程相减得于是(3)由此得到设其解为v=u+u,即设(2)考虑对线性方程组（1）的第二类扰动方程组（4）式右端给出了解的相对误差界。进一步得到(4)§3.7条件数、病态方程组与算法的稳定性两方程相减得于是(3)由此得到设其解为v=u+u,即设373.7.2

条件数与方程组病态性定义设A是一个非奇异矩阵，||.||是与某种向量范数相容的矩阵范数，正数cond(A):=||A-1||||A||称为矩阵A的

条件数（conditionalnumber）。Remark:由上面分析可知，系数矩阵的条件数很大(>>1)（此时矩阵也称为“病态矩阵”）的线性方程组Ax=b是病态的；反之，系数矩阵的条件数很小的线性方程组Ax=b是良态的。>>A=hilb(n);>>c=cond(A)

Example:

著名的Hilbert矩阵是病态的。【Matlab】§3.7条件数、病态方程组与算法的稳定性3.7.2条件数与方程组病态性定义设A是一个非38定理3.6

(条件数的性质)

设A是一个非奇异矩阵.(1)cond(A)≥1,cond(A)=cond(A-1),con(aA)=cond(A),a≠

0.(2)如果Q是正交矩阵，则cond2(Q)=1,cond2(A)=cond2(QA)=cond2(AQ),其中cond2(A)=||A-1||2||A||2.

§3.7条件数、病态方程组与算法的稳定性定理3.6(条件数的性质)设A是一个非奇异矩阵39§3.8稀疏矩阵的计算例（比较稀疏矩阵算法与普通矩阵算法：机器时间与占用空间）考虑线性方程组§3.8稀疏矩阵的计算例（比较稀疏矩阵算法与普通矩阵算法40%sparseandfullmatrixn=100;row=2:n;col=1:n-1;val=ones(1,n-1);offd=sparse(row,col,val,n,n);a=sparse(1:n,1:n.4*ones(1,n),n,n);b=full(a);rh=[1:n];tic;x=a/rh;t1=toctic;y=b/rh;t2=toc【实验要求】逐步增大n的值，观察发生的现象。§3.8稀疏矩阵的计算%sparseandfullmatrix【实验要求】41Ch3线性方程组求解的数值方法NumericalSolutionsofSystemsofLinearEquationsCh3线性方程组求解的数值方法NumericalSolu42§3.0Introduction线性方程组工程问题（电子网络，船体放样等）自然科学（实验数据的最小二乘法曲线拟合）解非线性方程组解常/偏微分方程组（差分法，有限元法）§3.0Introduction线性方程组工程问题（电子网433.0Introduction线性方程组的系数矩阵：（1）低价稠密矩阵（阶数<150）；（2）大型稀疏矩阵（阶数高且零元素较多）。求解线性方程组的方法：（1）直接法:经过有限步算术运算，求得精确解（假设计算过程没有舍入误差）。如Gauss消去法，三角分解法。（2）间接法(迭代法)：通过迭代序列，逐步逼近方程组的解。如Jacobi迭代法，Gauss-Seidel迭代法。3.0Introduction线性方程组的系数矩阵：（1）443.0Introduction【最简单的情形】三角形线性方程组：简记为，A为上三角矩阵。若所有aii0，可用“回代”过程得到方程组的解。3.0Introduction【最简单的情形】三角形线性方45functionX=backsub(U,b)%解上三角方程组--回代过程%BacksubstitutioninGausselimination%Input--Uisanxnupper-trianglularmatrix%--bisanx1constantvctor%Output--XisthesolutionvectorofUX=b%Usage:X=backsub(U,b)%FindthedimensionofbandinitializeXn=length(b);X=zeros(n,1);X(n)=b(n)/U(n,n);forp=n-1:-1:1X(p)=(b(p)-U(p,p+1:n)*X(p+1:n))/U(p,p);endBacksub.m3.0IntroductionfunctionX=backsub(U,b)Backsub463.0Introduction【作业】（1）编写Matlab程序解下列上三角形方程组：（2）编写Matlab程序解下列下三角形方程组：3.0Introduction【作业】（1）编写Matla47§3.1Gauss消去法与LU分解法直接法的基本思想：将线性方程组化成与之等价的上三角形或下三角形，再用回代法求解。它的核心是矩阵分解。核心：矩阵分解。§3.1Gauss消去法与LU分解法直接法的基本思想：将线48Gauss消去法（Gaussianelimination）：（对方程组的三种变换）（1）交换两个方程的次序；（2）用一个非零常数乘一个方程；（3）将一个方程的非零倍数加到另一个方程上去。这等价与对增广矩阵进行三种行初等变换，将它化为行阶梯形。Gauss消去法（Gaussianelimination）49【例3.1】解下列线性方程组：解：用行初等变换将方程组的增广矩阵化为行阶梯形：解得

x=(0,-1,1)’.§3.1GaussianEliminationandLUDecomposition【例3.1】解下列线性方程组：解：用行初等变换将方程组的增广50【注】上述过程可用矩阵表示为：LU=PA其中§3.1GaussianEliminationandLUDecomposition【注】上述过程可用矩阵表示为：LU=PA其中§3.1Gau51【高斯消元法思路】(1)用消元法将A化为上三角阵

upper-triangularmatrix;(2)回代求解backwardsubstitution。§3.1GaussianEliminationandLUDecomposition【高斯消元法思路】(1)用消元法将A化为上三角阵upp52§3.1GaussianEliminationandLUDecomposition消元记Step1：设

，计算因子将增广矩阵/*augmentedmatrix*/第i行mi1

第1行，得到其中Stepk：设，计算因子且计算共进行？步n

1§3.1GaussianEliminationand53回代Whatif?Nouniquesolutionexists.Whatif?Thenwemustfindthesmallestintegerkiwith,andinterchangethek-throwwiththei-throw.Whatifwecan’tfindsuchk

?Nouniquesolutionexists.定理若A的所有顺序主子式

/*determinantofleadingprincipalsubmatrices*/

均不为0，则高斯消元无需换行即可进行到底，得到唯一解。注：事实上，只要A非奇异，即A1存在，则可通过逐次消元及行交换，将方程组化为三角形方程组，求出唯一解。§3.1GaussianEliminationandLUDecomposition回代Whatif?Noun54高斯消元法的矩阵表达----矩阵分解课本P40-42§3.1GaussianEliminationandLUDecomposition高斯消元法的矩阵表达----矩阵分解课本P40-42§3.1553.1.2矩阵的LU分解（LUDecomposition）若在Gauss消去法过程中，没有进行交换行的变换，即P=I,则分解结果为§3.1GaussianEliminationandLUDecomposition这样的分解式称为矩阵的LU分解。【Matlab】>>[L,U,P]=lu(A)其中【作业】P.67题1,23.1.2矩阵的LU分解（LUDecomposition563.1.3选主元消去法

/*PivotingStrategies*/例：单精度解方程组/*精确解为和*/8个8个用GaussianElimination计算：8个小主元/*Smallpivotelement*/可能导致计算失败。§3.1GaussianEliminationandLUDecomposition3.1.3选主元消去法/*PivotingStra57

全主元消去法

/*CompletePivoting*/每一步选绝对值最大的元素为主元素，保证。Stepk:

①选取②Ifik

k

then交换第k行与第ik

行;Ifjk

k

then交换第k列与第jk

列;③消元注：列交换改变了xi的顺序，须记录交换次序，解完后再换回来。列主元消去法

/*PartialPivoting,ormaximalcolumnpivoting*/省去换列的步骤，每次仅选一列中最大的元。§3.1GaussianEliminationandLUDecomposition全主元消去法/*CompletePivoting58例：注：列主元法没有全主元法稳定。例：注意：这两个方程组在数学上严格等价。标度化列主元消去法/*ScaledPartialPivoting*/对每一行计算。为省时间，si只在初始时计算一次。以后每一步考虑子列中最大的aik为主元。注：稳定性介于列主元法和全主元法之间。§3.1GaussianEliminationandLUDecomposition例：注：列主元法没有全主元法稳定。例：注意：这两个方程组在59实际应用中直接调用GaussElimination解3阶线性方程组的结果：结合全主元消去后的结果：§3.1GaussianEliminationandLUDecomposition实际应用中直接调用GaussElimination解3阶60§3.2Cholesky分解【Matlab】设A是对称正定矩阵，则A有以下形式的LU分解：A=PTP（3.2.1）其中P是一个上三角矩阵,上式称为A的Cholesky分解(Choleskiandecomposition).>>P=chol(A)定理【作业】P.68题4§3.2Cholesky分解【Matlab】设A是对称正定61§3.3向量范数与矩阵范数—误差的度量3.3.1

向量范数(vectornorms)常用向量范数：这些范数满足：一般范数：§3.3向量范数与矩阵范数—误差的度量3.3.1向量范62如，设向量x=(2,-4,-1)T,则§3.3向量范数与矩阵范数【Matlab】x=(1:4)/5N1=norm(x,1)N2=norm(x)N7=norm(x,7)Ninf=norm(x,inf)x=0.20000.40000.60000.8000N1=2N2=1.0954N7=0.8153Ninf=0.8000如，设向量x=(2,-4,-1)T,则§3.3向63定义

矩阵A的范数定义为§3.3向量范数与矩阵范数3.3.2

矩阵范数(matrixnorms)

命题：矩阵的常用范数：其中，λ1,…,λ

n是的特征值，称为的谱半径。定义矩阵A的范数定义为§3.3向量范数与矩阵范数3.3.64例

设矩阵§3.3向量范数与矩阵范数则可以证明，矩阵范数满足：例设矩阵§3.3向量范数与矩阵范数则可以证明，矩阵范数满65§3.4经典迭代法－Jacobi迭代法与Gauss-Seidel迭代法3.4.1．Jacobi迭代法设有n阶方程组(3.4.1)

§3.4经典迭代法－Jacobi迭代法与Gauss-Sei66若系数矩阵非奇异，且(i=1,2,…,n)，将方程组（4.1）改写成§3.4经典迭代法若系数矩阵非奇异，且(67然后写成迭代格式（3.4.2）（3.4.2）式也可以简单地写为（3.4.3）§3.4经典迭代法然后写成迭代格式（3.4.2）（3.4.2）式也可以简单地写68写成矩阵形式：A=LUDBJacobi迭代阵(3.4.6)§3.4经典迭代法写成矩阵形式：A=LUDBJacobi迭代阵(3.4.669…………只存一组向量即可。写成矩阵形式：BGauss-Seidel迭代阵3.4.2Gauss-Seidel迭代法(3.4.7)§3.4经典迭代法…………只存一组向量即可。写成矩阵形70定理3.1

对任意初始向量X(0)及常向量F，迭代格式(3.5.1)(3.5.1)推论：若迭代矩阵Ｂ的某种范数，则（3.5.1）收敛的充分必要条件是迭代矩阵B的谱半径(B)<1。确定的迭代法对任意初值X(0)均收敛于方程组X=BX+F§3.5迭代法的分析——迭代法的收敛性、收敛速度与数值稳定性3.5.1收敛性的唯一解x*。定理3.1对任意初始向量X(0)及常向量F，迭代格式(3.71定义1：如果矩阵的每一行中，不在主对角线上的所有元素绝对值之和小于主对角线上元素的绝对值，即

则称矩阵A(按行)严格对角占优。定理3.2：若A是严格对角占优，则A是非奇异的。§3.5迭代法的分析定理3：若线性方程组AX=b的系数矩阵A按行严格对角占优，则Jacobi迭代法和Gauss-Seidel迭代法对任意给定初值均收敛。定义1：如果矩阵的每一行中，不在主对角线上的723.5.2迭代法的误差估计与收敛速度定理3.4

设X*是方程组AX=b的同解方程X=BX+F的准确解，若迭代公式中迭代矩阵B的某种范数，(1)(2)则有§3.5迭代法的分析3.5.2迭代法的误差估计与收敛速度定理73迭代加速（3.6.2）§3.6超松驰法(SOR)及分块迭代法定理3.５

SOR方法收敛的必要条件是其中ω称为松弛因子．（3.6.1）迭代加速（3.6.2）§3.6超松驰法(SOR)及分块迭代74§3.7条件数、病态方程组与算法的稳定性3.7.1

病态方程组与数值稳定的算法定义对线性方程组Ax=b,如果解x关于问题（矩阵A和b）的“微小”变化（即舍入误差）不敏感，则称此方程组是“好的”，或称为“良态的”，否则，称之为“坏的”，或称为“病态的”。定义对求解线性方程组Ax=b的一个(迭代)算法，如果关于初始值的误差，不随计算逐步放大，则称算法是数值稳定的，否则，称之为是数值不稳定的。数值不稳定的算法是不可用的！§3.7条件数、病态方程组与算法的稳定性3.7.1病态方75例（病态方程组）考虑线性方程组§3.7条件数、病态方程组与算法的稳定性已知其精确解是u