线性方程组AX=B的数值计算方法实验

《数值方法》实验报告1线性方程组AX=B的数值计算方法实验【摘要】在自然科学与工程技术中很多问题的解决常常归结为解线性代数方程组。例如电学中的网络问题，船体数学放样中建立三次样条函数问题，用最小二乘法求实验数据的曲线拟合问题，解非线性方程组的问题，用差分法或者有限元法解常微分方程，偏微分方程边值问题等都导致求解线性方程组。线性代数方面的计算方法就是研究求解线性方程组的一些数值解法与研究计算矩阵的特征值及特征向量的数值方法。关于线性方程组的数值解法一般有两类：直接法和迭代法。关键字高斯消元法、三角分解法、高斯-赛德尔迭代、稀疏矩阵一、实验目的1.掌握高斯消元法、三角分解法、高斯—赛德尔迭代发的编程技巧。2.掌握线性方程组AX=B的数值计算方法。3.掌握矩阵的基本编程技巧。二、实验原理1.高斯消元法数学上，高斯消元法是线性代数规划中的一个算法，可用来为线性方程组求解。高斯（Gauss）夏鸥按法其实是将一般的线性方程组变换为三角形（上三角）方程组求解问题（消元法），只是步骤规范，便于编写计算机程序。一般高斯消元法包括两过程：先把方程组化为同解的上三角形方程组，再按相反顺序求解上三角方程组。前者称为消去或消元过程，后者称回代过程。消去过程实际上是对增广矩阵作行初等变换。对一般的n阶方程组，消去过程分n-1步：第一步消去11a下方元素。第二步消去22a下方元素，......，第n-1步消去1-n1-na，下方元素。即第k步将第k行的适当倍数加于其后各行，或可说是从k+1~n行减去第k行的适当倍数，使它们第k列元素变为零，而其余列元素减去第k行对应列元素的倍数。2.三角分解法《数值方法》实验报告2三角分解法是将原正方(square)矩阵分解成一个上三角形矩阵或是排列(permuted)的上三角形矩阵和一个下三角形矩阵，这样的分解法又称为LU分解法。它的用途主要在简化一个大矩阵的行列式值的计算过程，求反矩阵，和求解联立方程组。不过要注意这种分解法所得到的上下三角形矩阵并非唯一，还可找到数个不同的一对上下三角形矩阵，此两三角形矩阵相乘也会得到原矩阵。3.高斯—赛德尔迭代高斯－赛德尔迭代（Gauss–Seidelmethod）是数值线性代数中的一个迭代法，可用来求出线性方程组解的近似值。研究雅可比迭代法，我们发现在逐个求）（1kX的分量时，当计算到）（1kiX时,分量）（1k1X，......，）（1k1-iX都已经求得，而仍用旧分量k1X，......，）（k1-iX计算）（1kiX。由于新计算出的分量比旧分量准确些，因此设想一旦新分量）（1k1X，......，）（1k1-iX求出，马上就用新分量）（1k1X，......，）（1k1-iX代替雅可比迭代法中k1X，......，）（k1-iX来求）（1kiX这就是高斯-赛德尔(Gauss-Seidel)迭代法。把矩阵A分解成ULDA(6)其中nna,...,a,adiagD2211,U,L分别为A的主对角元除外的下三角和上三角部分,于是,方程组(1)便可以写成bUxxLD即22fxBx其中bLDf,ULDB1212(7)以2B为迭代矩阵构成的迭代法(公式)221fxBxkk(8)称为高斯—塞德尔迭代法(公式),用变量表示的形式为,...,,k,n,,ixaxabaxijnij)k(jij)k(jijiii)k(i21021111111(9)4.稀疏矩阵《数值方法》实验报告3矩阵中非零元素的个数远远小于矩阵元素的总数，并且非零元素的分布没有规律，则称该矩阵为稀疏矩阵(sparsematrix)；与之相区别的是，如果非零元素的分布存在规律（如上三角矩阵、下三角矩阵、对角矩阵），则称该矩阵为特殊矩阵。常见于进行大量数据计算。三、实验内容1.P1081许多科学应用包含的矩阵带有很多的零。在实际情况中很重要的三角形线性方程组有如下形式：d1x1+c1x2=b1a1x1+d2x2+c2x3=b2a2x2+d3x3+c3x4=b3············aN-2xN-2+dN-1xN-1+cN-1xN=bN-1aN-1xN-1+dNxN=bN构造一个程序求解三角形线性方程组。可假定不需要行变换，而且可用第k行消去第k+1行的xk。2.P1201求解线性方程组AX=B，其中A=1357213500252631B=1234使用三角分解法求解X。3.P1202求解线性方程组AX=B，其中A=[aij]N×N，,aij=ij-1；而且B=[bij]N×1，b11=N，当i≥时，bij=(iN-1)/(i-1)。对N=3，7，11的情况分别求解。精确解为X=[11…11]’。对得到的结果与精确解的差异进行解释。4.P1203《数值方法》实验报告4通过重复求解N各线性方程组ACJ=EJ，其中J=1，2，…，N来得到A-1,则A[C1C2…CN]=[E1E2…EN]而且A-1=[C1C2…CN]保证对LU分解只计算一次!5.P1293设有如下三角线性方程组，而且系数矩阵具有严格对角优势：d1x1+c1x2=b1a1x1+d2x2+c2x3=b2a2x2+d3x3+c3x4=b3············aN-2xN-2+dN-1xN-1+cN-1xN=bN-1aN-1xN-1+dNxN=bN（i）根据方程组（1），式（2）和式（3），设计一个算法来求解上述方程组。算法必须有效地利用系数矩阵的稀疏性。a11x1+a12x2+···a1jxj+···+a1NxN=b1a21x1+a22x2+···a2jxj+···+a2NxN=b2···································方程组（1）aj1x1+aj2x2+···ajjxj+···+ajNxN=bj···············aN1x1+aN2x2+···aNjxj+···+aNNxN=bN(1)kjx=()()()()111111······kkkkjjjjjjjjjNNjjbaxaxaxaxa，j=1，2，···，N···式（2）(1)kjx=(1)(!)()()111111······kkkkjjjjjjjjjNNjjbaxaxaxaxa，j=1，2，···，N··式（3）（ii）根据（i）中设计的算法构建一个C++程序，并求解下列上三角线性方程组（a）4m1+m2=3（b）4m1+m2=1《数值方法》实验报告5m1+4m2+m3=3m1+4m2+m3=2m2+4m3+m4=3m2+4m3+m4=1m3+4m4+m5=3m3+4m4+m5=2························m48+4m49+m50=3m48+4m49+m50=1m49+m50=3m49+m50=26.P1294利用高斯—赛德尔迭代法求解下列带状方程。12x1-2x2+x3=5-2x1+12x2-2x3+x4=5x1-2x2+12x3-2x4+x5=5x2-2x3+12x4-2x5+x6=5···x46-2x47+12x48-2x49+x50=5x47-2x48+12x49-2x50=5x48-2x49+12x50=5四、实验结果及分析1.P1081实验描述：本次实验使用系数矩阵的第k行消去第k+1行的xk，消除方法为第k行减去第k-1行乘上系数ak-1/bk-1，待消至第N行时，求解出xN，并依次会带求出各xN-1至x1，为了检验结果的正确性使用上面的方程组组（1）及方程组（2）进行验证。方程组（1）的结果为12341322xxxx，方程组（2）的结果为12342321xxxx。算法流程图：《数值方法》实验报告6startInputA,B,P,delta,max1N=length(B);k=1kmax1j=1jNj==1j==NX(1)=(B(1)-A(1,2)*P(2))/A(1,1)X(N)=(B(N)-A(N,N-1)*(X(N-1))')/A(N,N)X(j)=(B(j)-A(j,j-1)*X(j-1)'-A(j,j+1)*P(j+1))/A(j,j)j=j+1YNYNNNYYk=k+1《数值方法》实验报告7实验结果：结果截图(1)errdelta)|(relerrdeltaoutputenderr=abs(norm(X'-P));relerr=err/(norm(X)+eps);P=X';YN《数值方法》实验报告8图1输入矩阵Ab=1272319423102412，输出结果为X=1322，与预期结果一致（2）图2《数值方法》实验报告9输入矩阵Ab=115215934219262，输出结果为X=2321，与预期结果一致实验结论：通过对系数矩阵的增广矩阵进行高斯消元和回带容易得到线性方程组的解，同时，利用这种方法可以求得矩阵的逆。2.P1201实验描述：本次实验的解法为使用LU矩阵求解X，该解法的内容为将系数矩阵A分解为一上三角矩阵及一下三角矩阵，且有A=LU，之后由LY=B，UX=Y分别求解出Y，X。实验结果：图3输入矩阵A=1357213500262631，B=1234，输出结果为X=1.342860.68571431.8实验结论：将X代入AX后结果与矩阵B一致，运行结果正确无误，该程序正确，且有《数值方法》实验报告10X=1.342860.68571431.83.P1202实验描述：（1）本次实验仍使用三角矩阵求解矩阵X的值，求解方法与实验二大致相同；（2）矩阵A编写函数buildA完成，buildA的输入为矩阵A的阶数N，输出结果为矩阵A的地址，对于矩阵A有A=[aij]N×N，,aij=ij-1；（3）矩阵B编写函数buildB完成，buildB的输入为矩阵B的阶数N，输出结果为矩阵B的地址，对于矩阵B有B=[bij]N×1，b11=N，当i≥时，bij=(iN-1)/(i-1)。实验结果：（1）N=3时图4（2）N=7时《数值方法》实验报告11图5（3）N=11时图6实验结论：由图4，图5，图6与结果对比可知，程序运行结果与预期结果相一致，程序正确无误。4.P1203实验描述：（1）本次实验的目的为求解A的逆矩阵A-1，求解方法为利用A[C1C2…CN]=[E1E2…EN]，A-1=[C1C2…CN]求解，故可将之分解为对于ACk=Ek，k=1，2，3，···，N中对于Ck的求解，之后另A-1=[C1C2…CN]；(2)由于A-1=[C1C2…CN]，A[C1C2…CN]=[E1E2…EN]，固有[E1E2…EN]=I，故Ek=[a1j]，j=1，2，···，N，其中a1k=1，other=0；(3)对于ACk=Ek，k=1，2，3，···，N的求解方法使用LU三角分解法求《数值方法》实验报告12解，用此方法求解出各个Ek对应的Ck，最后以此构成A-1；(4)求出LU后，应判断矩阵对角线上是否存在为0的元素，若存在，则A不存在逆矩阵；若不存在，则可求解逆矩阵A-1；(5)上述方法中的LU分解只需要进行一次；(6)对于程序的正确性使用矩阵132264251及矩阵123213231进行验证，其中，矩阵132264251不存在逆矩阵，矩阵123213231的逆矩阵为-0.66670.58330.250.3333-0.41670.250.33330.0833-0.25实验结果：（1）输入为矩阵132264251图