matlab解线性方程组

MATLAB数学建模与仿真线性方程组目的了解线性方程组的数值解法掌握求解线性方程组迭代法的有关原理方法用迭代法的有关理论分析迭代法的收敛性和收敛速度预备知识背景工程中的许多问题最后都可以转化为求解线性方程组许多数值计算问题（如样条函数、常微分方程数值解、差分方程等）的研究也往往归结为此类问题线性方程组的数值法直接法——经过有限次算法运算求出精确解，最常用的是：高斯消元法、矩阵LU分解迭代法——从初值出发，用递推的方法，给出近似解序列。常用的方法：雅可比迭代法、高斯-赛德尔迭代法。预备知识适用情形直接法一般适用于系数矩阵A为低阶稠密矩阵(非零元素较多)的情形；在工程技术和科学计算中常会遇到大型稀疏矩阵(非零元素较少)形式的方程组；迭代法在计算和存储都适合这一情形。解线性方程组的迭代法直接法得到的解是理论上准确的，但是我们可以看得出，它们的计算量都是n3数量级，存储量为n2量级，这在n比较小的时候还比较合适（n1000,1G/s,15秒,8M），但是对于现在的很多实际问题，往往要我们求解很大的n的矩阵，而且这些矩阵往往是系数矩阵就是这些矩阵含有大量的0元素。对于这类的矩阵，在用直接法时就会耗费大量的时间和存储单元。因此我们有必要引入一类新的方法：迭代法。迭代法具有的特点是速度快。与非线性方程的迭代方法一样，需要我们构造一个等价的方程，从而构造一个收敛序列，序列的极限值就是方程组的根对方程组bAx做等价变换gGxxbMNxMxNxbMxbxNMbAx11)(如：令NMA，则则，我们可以构造序列gxGxkk)()1(若*)(xxkbAxgxGx***同时：*)(**)()()1(xxGGxGxxxkkk*)()0(1xxGk0kG所以，序列收敛与初值的选取无关定义：（收敛矩阵）0kG定理：矩阵G为收敛矩阵，当且仅当G的谱半径11)(0GGk由GG)(知，若有某种范数1pG则，迭代收敛Jacobi迭代nnnnnnnbxaxabxaxa1111111)(1)(1)(111112132312122211212111nnnnnnnnnnnnbxaxaaxbxaxaxaaxbxaxaax)(1)(1)(1)(11)(11)1(2)(1)(323)(12122)1(21)(1)(21211)1(1nknnnknnnknknnkkkknnkkbxaxaaxbxaxaxaaxbxaxaax格式很简单：)(11)(11)()1(inijkjijijkjijiikibxaxaaxJacobi迭代算法1、输入系数矩阵A和向量b，和误差控制eps2、x1={0,0,…..,0},x2={1,1,…..,1}//赋初值3、while(||A*x2-b||eps){x1=x2;for(i=0;in;i++){x2[i]=0;for(j=0;ji;j++){x2[i]+=A[i][j]*x1[j]}for(j=i+1;jn;j++){x2[i]+=A[i][j]*x1[j]}x2[i]=-(x2[i]-b[i])/A[i][i]}}4、输出解x2迭代矩阵ULDA记nnaaD0011000001121nnnaaaL000001112nnnaaaU易知，Jacobi迭代有bxULD)(bxULDx)(bDxULDx11)(bDgADIULDG111,)(收敛条件迭代格式收敛的充要条件是G的谱半径1。对于Jacobi迭代，我们有一些保证收敛的充分条件定理：若A满足下列条件之一，则Jacobi迭代收敛。①A为行对角占优阵②A为列对角占优阵③A满足ijijiiaajiijjjaa1jiiiijaa证明：)(1ULDGiiijijijiiijiaaaaG1max1max1jiiiijiaaG②A为列对角占优阵，则AT为行对角占优阵，有1)(1TADI1)()(11TADIADI＃证毕)(1)(1)(1)1(11)1(11)1(2)(1)(323)1(12122)1(21)(1)(21211)1(1nknnnknnnknknnkkkknnkkbxaxaaxbxaxaxaaxbxaxaaxGauss－Seidel迭代)(11)(11)1()1(inijkjijijkjijiikibxaxaax在Jacobi迭代中，使用最新计算出的分量值Gauss-Siedel迭代算法1、输入系数矩阵A和向量b，和误差控制eps2、x2={1,1,…..,1}//赋初值3、while(||A*x2-b||eps){for(i=0;in;i++){for(j=0;ji;j++){x2[i]+=A[i][j]*x2[j]}for(j=i+1;jn;j++){x2[i]+=A[i][j]*x2[j]}x2[i]=-(x2[i]-b[i])/A[i][i]}}4、输出解x2迭代矩阵)()()1(1)1(bUxLxDxkkkbDUxDxLDIkk1)(1)1(1)(bDLDIUxDLDIxkk111)(111)1()()(bLDUxLDxkk1)(1)1()()(bLDgULDG11)(,)(是否是原来的方程的解？A=(D-L)-U收敛条件迭代格式收敛的充要条件是G的谱半径1。我们看一些充分条件定理：若A满足下列条件之一，则Gauss-Seidel迭代收敛。①A为行或列对角占优阵②A对称正定阵证明：设G的特征多项式为)(sP，则ULDLDULDIGIPs)()()()(11ULDA为对角占优阵，则1时ULD)(为对角占优阵0)(ULD即0)(sP1即1)(G＃证毕注：二种方法都存在收敛性问题。有例子表明：Gauss-Seidel法收敛时，Jacobi法可能不收敛；而Jacobi法收敛时，Gauss-Seidel法也可能不收敛。(0)1807109,8,(0,0,0)'9117Abx1、预处理9117180,71098Ab2、格式(1)()()123(1)(1)11(1)(1)111/9(7)1/8(7)1/9(8)kkkkkkkxxxxxxx(1)(2)(3)(4)(0.7778,0.9722,0.9753)'(0.9942,0.9993,0.9994)'(0.9999,0.9999,0.9999)'(1.0000,1.0000,1.0000)'xxxx3、结果101/21/2()1011/21/20BDLU1、Jacobi迭代211111112A特征值为3504IB12,350,2i101/21/2()01/21/2001/2BDLU2、Gauss-Seidel迭代12,310,2松弛迭代记)()1()(kkkxxx则)()()1(kkkxxx可以看作在前一步上加一个修正量。，有)()()1(kkkxxx)()()1(1)1(bUxLxDxkkk对Gauss-Seidel迭代格式)()()1()()1(kkkkxxxx)()(1)(1)1(1)()1(kkkkkxbDUxDLxDxx)()1(1)(1)1(1)()1(bDUxDLxDxxkkkk若在修正量前乘以一个因子)(1)1()()1()(1)1()(11)(11)()1(2)(1)(323)(12122)(2)1(21)(1)(21211)(1)1(1nknnnknnnknknknnkkkkknnkkkbxaxaaxxbxaxaxaaxxbxaxaaxx写成分量形式，有松弛迭代算法1、输入系数矩阵A、向量b和松弛因子omega，和误差控制eps2、x2={1,1,…..,1}//赋初值3、while(||A*x2-b||eps){for(i=0;in;i++){temp-0for(j=0;ji;j++){temp+=A[i][j]*x2[j]}for(j=i+1;jn;j++){temp+=A[i][j]*x2[j]}temp=-(x2[i]-b[i])/A[i][i]x2[i]=(1-omega)*x2[i]+omega*temp}}4、输出解x2)()1(1)(1)1(1)()1(bDUxDLxDxxkkkk迭代矩阵bDxUDIxLDIkk1)(1)1(1))1(()(bDLDIxUDILDIxkk111)(111)1()())1(()(bDLDIgUDILDIG111111)())1(()(定理：松弛迭代收敛20定理：A对称正定，则松弛迭代收敛20是否是原来的方程的解？SOR方法收敛的快慢与松弛因子的选择有密切关系.但是如何选取最佳松弛因子,即选取=*,使(G)达到最小,是一个尚未很好解决的问题.实际上可采用试算的方法来确定较好的松弛因子.经验上可取1.41.6.LAB06线性方程组求根的迭代法1.编写Gauss-Seidel迭代和SOR迭代的通用程序2.用如上程序求根，并打印迭代步数和根。3113000100001335901100000931100000000107930000900030577050000074730000000030410000005002720009000229A152723020127710Ax3.取松弛因子为{i/50,i=1,2,…,99}，试给出一个最佳的值定理若SOR方法收敛,则02.证设SOR方法收敛,则(G)1,所以|det(G)|=|12…n|1而det(G)=det[(D-L)-1((1-)D+U)]=det[(E-D-1L)-1]det[(1-)E+D-1U)]=(1-)n于是|1-|1,或02定理设A是对称正定矩阵,则解方程组Ax=b的SOR方法,当02时收敛.证设是G的任一特征值,y是对应的特征向量,则[(1-)D+U]y=(D-L)y于是(1-)(Dy,y)+(Uy,y)=[(Dy,y)-(Ly,y)])()()())(1(yLy,yDy,yUy,yDy,由于A=D-L-U是对称正定的,所以D是正定矩阵,且L=UT.若记(Ly,y)=+i,则有(Dy,y)=0(Uy,y)=(y,Ly)=(Ly,y)=-i0(Ay,y)=(Dy,y)-(Ly,y)-(Uy,y)=