解线性方程组

解线性方程组的直接法对方程AXb,矩阵A的维数是mn,b为一给定的m维向量.则有(1)mn恰定方程,寻求精确解;(2)mn超定方程,寻求最小二乘解;(3)mn欠定方程，寻求基本解,其中至多有m个非零元素.针对不同的况,MATLAB将采用不同的算法来求解.一、超定方程组当mn时,方程组个数比未知数多,上式为一超定方程组.一般而言,满足方程AXb的解将不存在,方程没有精确解,在MATLAB中,利用左除命令\XAb来寻求它的最小二乘解,即使2AXb最小.还可以用广义逆来求,即()*XpinvAb,所得的解不一定满足AXb,X只是最小二乘意义上的解.左除的方法是建立在奇异值分解基础之上,由此获得的解最可靠;广义逆法是建立在对原超定方程直接进行householder变换的基础上,其算法可靠性稍逊于奇异值分解,但速度快.二、欠定方程组当mn时,未知数比方程个数多,上式为一欠定方程组.这时,满足方程AXb的解有无穷多个,MATLAB只给出其中一个基本解\XAb,要得到它的通解,可用函数null来实现三、恰定方程组恰定方程组由n个未知数,n个方程构成,对AXb在线性代数教科书中,最常用的方程解法有:(1)利用cramer公式求解法;(2)利用矩阵求逆解法,即1XAb;(3)利用gauss消去法;(4)利用LU分解法求解.一般来说,对于维数不高、条件数不大的矩阵,上面4种解法所得的结果差别不大.前两种解法的真正意义是在其理论上,而不是实际的数值计算.而Gauss消去法,其本质上利用LU分解,在MATLAB中,出于对算法稳定性的考虑,行列式及逆矩阵的计算大都在LU分解的基础上进行.因此,在MATLAB中,求解这类方程组时可直接采用表达式:\XAb.(LU分解)若n阶矩阵A可逆且顺序主子式不为零,则A可以分解为一个单位下三角阵L和一个上三角阵U的积ALU,并且这种分解是唯一的.由于AXLUXb记UXZ,则LZb,从而由LZb求得\ZLb,再由UXZ求\XUZ,\(\)XULbMATLAB中,用[,]()LUluA函数求得L,U,再用\(\)XULb求得解.如果矩阵A是对称正定矩阵,可采用Cholesky分解法,矩阵Cholesky分解定理为:如果A是对称正定矩阵,则(至少)存在一个实的下三角矩阵L使得TALL此外,我们可以限定矩阵L的对角元素全部为正,那么,对应的分解TALL是唯一的.在MATLAB中,用()LcholA函数求得L,再用\('\)XLLb求得解.Guass法(LU分解)和cholesky法的基础都是把线性方程组的矩阵分解为下三角矩阵和上三角矩阵的乘积,但对于对称正定矩阵的情形,cholesky比Gauss法更加简便.我们也可以使用Householder法则把矩阵分解为正交矩阵Q和上三角矩阵R的乘积,称为QR因子分解法.给定n阶矩阵A,则存在一个酉矩阵Q,以及一个上三角矩阵R,使得AQR此外,我们可以设法使矩阵R的对角元素都为正.如果A是可逆的,则这时所对应的分解AQR是唯一的.在MATLAB中,用[,]()QRqrA函数求得Q,R,再用\(\)XRQb求得解.线性方程组的最小二乘解线性最小二乘问题的解法在数据拟合、测量平差、控制理论等方面均得到广泛的应用。例如,已知m对数据(,)iity（其中1,2,,im）和n个已知函数()jht（其中1,2,,jn）试构造线性组合1()njjjxht。用该线性组合最佳地拟合这m对数据(,)iity（其中1,2,,im）即我们希望适当地选取组合系数jx（1,2,,jn）,使得在某种范数意义下,误差1()(),(1,2,,)niijjijrxyxhtim能够达到最小。令(),,(1,2,,;1,2,,)jiijiihtaybimjn。则上式可用矩阵向量形式把误差向量表为()rxbAx，其中11112122122212()()(),()nnmnnnnrxaaarxaaarxArxaaa1212(,,,),(,,,)TTmnbyyyxxxx当mn时,在上式中可要求()0rx，则估计12(,,,)Tnxxxx的问题就转化为求解线性方程组。当mn时,一般()0rx，最小二乘问题就是适当选取x使误差()rx在2范数意义下等于最小。给定矩阵mnAC及向量mbC,寻找nxC满足222()minnyCrxbAxAyb这就是线性最小二乘问题,其解也称为线性方程组,mnAxbAC的最小二乘解。在方程组不相容情形,它可视为方程组在最小二乘意义下的最优近似解。在方程组相容时,则最小二乘解与通常意义下的解是一致的。奇异值分解1、奇异值与奇异值分解定理奇异值定理：设mnAC，r=rank(A)，则一定存在m阶酉矩阵U和n阶酉矩阵V和对角矩阵1212(,,,)(1,2,,)rridiagir，且，而，使得HAUV，称为A的奇异值分解。复数域内的奇异值：设(0)mnHrACrAA，的特征值为1210rrn则称(1,2,,)iiin为A的正奇异值；当A为零矩阵时，它的奇异值都是零。易见，矩阵A的奇异值的个数等于A的列数，A的非零奇异值的个数等于rank(A)。2、奇异值分解的理解奇异值σ跟特征值类似，在矩阵Σ中也是从大到小排列，而且σ的减少特别的快，在很多情况下，前10%甚至1%的奇异值的和就占了全部的奇异值之和的99%以上了。也就是说，我们也可以用前r大的奇异值来近似描述矩阵。rrrTrrrTvuvuvuVUVUA222111即A可表示为r个秩为一的举证的和，这是A得奇异值分解的紧凑格式。3、奇异值分解特征奇异值分解的第一个特征是可以降维。A表示n个m维向量,通过奇异值分解可表示成m+n个r维向量,若A的秩r远远小于m和n,则通过奇异值分解可以大大降低A的维数。可以计算出,当rmn/(m+n+1)时,可以达到降维的目的,同时可以降低计算机对存贮器的要求。奇异值分解的第二个特征是奇异值对矩阵的扰动不敏感,而特征值对矩阵的扰动敏感。在数学上可以证明,奇异值的变化不会超过相应矩阵的变化,即对任何相同阶数的实矩阵A、B,按从大到小排列的奇异值i和i,有2BAiii.奇异值的第三个特征是奇异值的旋转不变性。即若P是正交阵,PA的奇异值与A的奇异值相同。奇异值的比例和旋转不变性特征在数字图像的旋转、镜像、平移、放大、缩小等几何变化方面有很好的应用。奇异值的第四个特征是容易得到矩阵A的秩为k(k≤r)的一个最佳逼近矩阵。奇异值的这个特征可以应用于信号的分解和重构,提取有用信息,消除信号噪声。讨论理解判别收敛条件1．收敛性问题现在来研究与方程组fBxx对应的基本型迭代公式),2,1,0()()1(kfBxxkk设*x是方程组fBxx（也即bAx）的就解，即有fBxx**要研究由迭代公式fBxxkk)()1(产生的序列)(kx当k时是否收敛于*x,)()(*)1((*)kkxxBxx=)()1(*2kxxB=)()0(*1xxBk可见当k时，是否有*)(xxk，等价于是否有0kB（零矩阵，即kB的每一个元素趋于零）根据线性代数，任何n阶矩阵B都存在非奇异矩阵P，使得JPPB1PJPBkk1其中J为B的Jordan标准形1JJ2JnnrJkkJJ1kJ2nnkrJiiJ1iiinni1nnrii1于是，可得),,2,1;(0)(0)(0rikJkJkBkikk这时，若设2J13J11mJ111kkJ2kkk1=kkkkc11kkJ3kkk1kkkkkk122/)1(kkkkc11kkkkkcc1122kkmJkkkc111122kkkkcckmkkmmkkmcc2211其中)!(!!mkmkckm于是又可得到),,2,1(1)(0rikJiki注意到),,2,1(1ri，即1)(B，于是最后可得1)()(0BkBk由上面就得到下面迭代法的收敛定理。2．收敛性基本定理定理1（迭代法收敛性基本定理）设方程组为fBxx对任意的初始向量)0(x，解此方程组的迭代法),2,1,0()()1(kfBxxkk收敛的充分必要条件是迭代矩阵B的谱半径1)(B注意：此定理为判断迭代法的敛散性提供了一个强有力的手段（充分必要条件）。然而，定理的条件往往不容易验证。因此，利用特征值上界性质BB)(,可以给出另一个条件较弱的结果。定理2（迭代法收敛性充分条件）如果迭代法)2,1,0()()1(kfBxxkk的迭代矩阵B的某一种算子范数1B，则（1）对任意的初始向量)0(x,迭代法收敛；（2）迭代序列与方程组的解*x存在误差估计式)1()()(*1kkkxxBBxx或)0()1()(*1xxBBxxkk证明（1）由条件1B及BB)(从而1)(BB，按定理3.3.2可知迭代法收敛。（2）)()()1(*)1()(*kkkkxxxxxx=)()()(*)1()(kkkxxBxxB)(*)1()(kkkxxBxxB整理得)1()()(*1kkkxxBBxx由)()2()1()1(kkkkxxBxx得)2()1()1()(kkkkxxBxx则可得(2).3．其他定理设nnijRaA)(，若),,2,1(||||1niaanijjijii则称A为严格对角占优矩阵。容易证明，若nnijRaA)(严格对角占优，则),,2,1(0niaii，且A为非奇异。事实上，前者由定义可得到。后者用反证法，设A奇异，则有0x满足0Ax。又设其中kx为0||xxk，则0Ax的第k个方程为nkjjjkjkkkxaxa1由此有nkjjkjnkjjkjkjkkaxxaa11||||||||定理3若方程组bAx中nnijRaA)(为严格对角占优矩阵，则解此方程组的J法和GS法收敛。定理3.3.4若方程组bAx中nnijRaA)(为对称正定矩阵，则解此方程组的GS法收敛。4．收敛速度由定理2（2）式可知（1）)0(x越接近*x，迭代法收敛速度越快，即收敛速度与初值选取有关；(2)迭代矩阵的范数B越小，(当然应有1B)，迭代法收敛也越快；（3）由于BB)(，所以迭代矩阵的谱半径)(B越小，迭代收敛越快。讨论设计算法及编程实现共轭斜量法根据书P175的算法表述可得算法流程图如下：开始输入A,bx=rand(3,1)r=b-Axp=rx=x+αpr=r-αpr==0?p=r-βp结束NY输出x利用Matlab编程如下：A=[1-.1-.2;-.11-.2;-.2-.21];b=[.72.83.84]';x=F(A,b)F.m:functionx=F(A,b)N=size(A,2);x=rand(N,1);r=b-A*x;p=r;fork=1:Naerf=sum(r.*p)/sum((A*p).*p);x=x+aerf*p;r=r-aerf*A*p;ifr==0break;endbeita=sum(r.*(A*p))/sum(A*p.*p);p=r