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

32数值方法实验报告线性方程组AX=B的数值计算方法实验【摘要】在自然科学与工程技术中很多问题的解决常常归结为解线性代数方程组。例如电学中的网络问题，船体数学放样中建立三次样条函数问题，用最小二乘法求实验数据的曲线拟合问题，解非线性方程组的问题，用差分法或者有限元法解常微分方程，偏微分方程边值问题等都导致求解线性方程组。线性代数方面的计算方法就是研究求解线性方程组的一些数值解法与研究计算矩阵的特征值及特征向量的数值方法。关于线性方程组的数值解法一般有两类：直接法和迭代法。关 键 字 高斯消元法、三角分解法、高斯-赛德尔迭代、稀疏矩阵1、 实验目的1.掌握高斯消元法、三角分解法、高斯赛德尔迭代发的编程技巧。2.掌握线性方程组AX=B的数值计算方法。3.掌握矩阵的基本编程技巧。2、 实验原理1. 高斯消元法 数学上，高斯消元法是线性代数规划中的一个算法，可用来为线性方程组求解。高斯（Gauss）夏鸥按法其实是将一般的线性方程组变换为三角形（上三角）方程组求解问题（消元法），只是步骤规范，便于编写计算机程序。 一般高斯消元法包括两过程：先把方程组化为同解的上三角形方程组，再按相反顺序求解上三角方程组。前者称为消去或消元过程，后者称回代过程。消去过程实际上是对增广矩阵作行初等变换。 对一般的n阶方程组，消去过程分n-1步：第一步消去下方元素。第二步消去下方元素，.，第n-1步消去下方元素。即第k步将第k行的适当倍数加于其后各行，或可说是从k+1n行减去第k行的适当倍数，使它们第k列元素变为零，而其余列元素减去第k行对应列元素的倍数。 2. 三角分解法三角分解法是将原正方 (square)矩阵分解成一个上三角形矩阵或是排列(permuted) 的上三角形矩阵和一个 下三角形矩阵，这样的分解法又称为LU分解法。它的用途主要在简化一个大矩阵的行列式值的计算过程，求 反矩阵，和求解联立方程组。不过要注意这种分解法所得到的上下三角形矩阵并非唯一，还可找到数个不同 的一对上下三角形矩阵，此两三角形矩阵相乘也会得到原矩阵。3. 高斯赛德尔迭代 高斯赛德尔迭代（GaussSeidelmethod）是数值线性代数中的一个迭代法，可用来求出线性方程组解的近似值。研究雅可比迭代法，我们发现在逐个求的分量时，当计算到时,分量，.，都已经求得，而仍用旧分量，.，计算。由于新计算出的分量比旧分量准确些，因此设想一旦新分量，.，求出，马上就用新分量，.，代替雅可比迭代法中，.，来求这就是高斯-赛德尔(Gauss-Seidel)迭代法。 把矩阵A分解成 (6) 其中,分别为的主对角元除外的下三角和上三角部分,于是,方程组(1)便可以写成 即其中 (7)以为迭代矩阵构成的迭代法(公式) (8)称为高斯塞德尔迭代法(公式),用变量表示的形式为 (9)4. 稀疏矩阵矩阵中非零元素的个数远远小于矩阵元素的总数，并且非零元素的分布没有规律，则称该矩阵为稀疏矩阵(sparse matrix)；与之相区别的是，如果非零元素的分布存在规律（如上三角矩阵、下三角矩阵、对角矩阵），则称该矩阵为特殊矩阵。常见于进行大量数据计算。3、 实验内容1. P108 1许多科学应用包含的矩阵带有很多的零。在实际情况中很重要的三角形线性方程组有如下形式： d1x1+c1x2 =b1 a1x1+d2x2+c2x3 =b2 a2x2+d3x3+c3x4 =b3 aN-2xN-2+dN-1xN-1+cN-1xN =bN-1 aN-1xN-1+dNxN =bN构造一个程序求解三角形线性方程组。可假定不需要行变换，而且可用第k行消去第k+1行的xk。2. P120 1求解线性方程组AX=B，其中 A= B=使用三角分解法求解X。3. P120 2求解线性方程组AX=B，其中A=aijNN，,aij=ij-1；而且B=bijN1，b11=N，当i时，bij=(iN-1)/(i-1)。对N=3，7，11的情况分别求解。精确解为X=1 1 1 1。对得到的结果与精确解的差异进行解释。4. P120 3通过重复求解N各线性方程组 ACJ=EJ，其中J=1，2，N来得到A-1,则 AC1 C2 CN=E1 E2 EN而且 A-1=C1 C2 CN保证对LU分解只计算一次!5.P129 3设有如下三角线性方程组，而且系数矩阵具有严格对角优势：d1x1+c1x2 =b1 a1x1+d2x2+c2x3 =b2 a2x2+d3x3+c3x4 =b3 aN-2xN-2+dN-1xN-1+cN-1xN =bN-1 aN-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=，j=1，2，N 式（2）=，j=1，2，N 式（3）（ii）根据（i）中设计的算法构建一个C+程序，并求解下列上三角线性方程组（a）4m1+ m2 =3 （b）4m1+ m2 =1 m1+4m2+m3=3 m1+4m2+m3=2 m2+4m3+m4=3 m2+4m3+m4=1m3+4m4+m5=3 m3+4m4+m5=2 m48+4m49+m50=3 m48+4m49+m50=1m49+m50=3 m49+m50=26. P129 4利用高斯赛德尔迭代法求解下列带状方程。 12x1 -2x2+x3 =5 -2x1+12x2-2x3+x4 =5 x1-2x2+12x3-2x4+x5=5 x2-2x3+12x4-2x5+x6=5 x46-2x47+12x48-2x49+x50=5 x47-2x48+12x49-2x50=5 x48-2x49+12x50=54、 实验结果及分析1. P108 1 实验描述：本次实验使用系数矩阵的第k行消去第k+1行的xk，消除方法为第k行减去第k-1行乘上系数ak-1/bk-1，待消至第N行时，求解出xN，并依次会带求出各xN-1至x1，为了检验结果的正确性使用上面的方程组组（1）及方程组（2）进行验证。方程组（1）的结果为，方程组（2）的结果为。算法流程图：start InputA,B,P,delta,max1N=length(B); k=1Yk=k+1kmax1Nj=1Yj=j+1jNYNX(1)=(B(1)-A(1,2)*P(2)/A(1,1)j=1YNX(N)=(B(N)-A(N,N-1)*(X(N-1)/A(N,N)j=NNX(j)=(B(j)-A(j,j-1)*X(j-1)-A(j,j+1)*P(j+1)/A(j,j)err=abs(norm(X-P); relerr=err/(norm(X)+eps);P=X;Nerrdelta)|(relerrdeltaYoutputend 实验结果：结果截图(1) 图1输入矩阵=，输出结果为X=，与预期结果一致（2） 图2输入矩阵=，输出结果为X=，与预期结果一致实验结论： 通过对系数矩阵的增广矩阵进行高斯消元和回带容易得到线性方程组的解，同时，利用这种方法可以求得矩阵的逆。2.P120 1 实验描述：本次实验的解法为使用LU矩阵求解X，该解法的内容为将系数矩阵A分解为一上三角矩阵及一下三角矩阵，且有A=LU，之后由LY=B，UX=Y分别求解出Y，X。 实验结果： 图3输入矩阵A=，B=，输出结果为X=实验结论：将X代入AX后结果与矩阵B一致，运行结果正确无误，该程序正确，且有X=3. P120 2 实验描述：（1）本次实验仍使用三角矩阵求解矩阵X的值，求解方法与实验二大致相同；（2）矩阵A编写函数buildA完成，buildA的输入为矩阵A的阶数N，输出结果为矩阵A的地址，对于矩阵A有A=aijNN，,aij=ij-1；（3）矩阵B编写函数buildB完成，buildB的输入为矩阵B的阶数N，输出结果为矩阵B的地址，对于矩阵B有B=bijN1，b11=N，当i时，bij=(iN-1)/(i-1)。 实验结果：（1） N=3时 图4（2） N=7时 图5（3） N=11时 图6 实验结论：由图4，图5，图6与结果对比可知，程序运行结果与预期结果相一致，程序正确无误。4. P120 3 实验描述：（1）本次实验的目的为求解A的逆矩阵A-1，求解方法为利用AC1 C2 CN=E1 E2 EN ，A-1=C1 C2 CN求解，故可将之分解为对于ACk=Ek，k=1，2，3，N中对于Ck的求解，之后另A-1=C1 C2 CN；(2)由于A-1=C1 C2 CN，AC1 C2 CN=E1 E2 EN，固有E1 E2 EN=I，故Ek=a1j，j=1，2，N，其中a1k=1，other=0；(3)对于ACk=Ek，k=1，2，3，N的求解方法使用LU三角分解法求解，用此方法求解出各个Ek对应的Ck，最后以此构成A-1；(4)求出LU后，应判断矩阵对角线上是否存在为0的元素，若存在，则A不存在逆矩阵；若不存在，则可求解逆矩阵A-1；(5)上述方法中的LU分解只需要进行一次；(6)对于程序的正确性使用矩阵及矩阵进行验证，其中，矩阵不存在逆矩阵，矩阵的逆矩阵为实验结果：（1）输入为矩阵 图7（2） 输入矩阵为 图8 实验结论： 如果系数矩阵能分解为LU的形式，其中L为下三角矩阵，U为上三角矩阵，通过对系数矩阵的分解再求解可应用简单的迭代进行求解x。5.P129 3 实验描述：（1）本实验的目的为使用高斯赛德尔迭代求解带状方程组；（2）由于实验中的带状方程有极强的稀疏性和相似性，故编程时应考虑该矩阵的以上特点以减少运算量及运行时占用的内存；（3）为了减少程序的运行次数，故选择式（3）作为运行的程序；（4）应无明确的结束标志，故选择delta=1e-5作为结束迭代的标志，此时再进行一轮迭代后输出结果。实验结果：a）4m1+ m2 =3 m1+4m2+m3=3 m2+4m3+m4=3m3+4m4+m5=3 m48+4m49+m50=3 m49+m50=3 图9（b）4m1+ m2 =1 m1+4m2+m3=2 m2+4m3+m4=1m3+4m4+m5=2 m48+4m49+m50=1m49+m50=2 图10 实验结论：求解带状线性方程组的解可使用高斯-赛德尔迭代法。6. P129 4 实验描述：（1）本实验的目的为使用高斯赛德尔迭代求解带状方程组；（2）由于实验中的带状方程有极强的稀疏性和相似性，故编程时应考虑该矩阵的以上特点以减少运算量及运行时占用的内存；（3）为了减少程序的运行次数，故选择式（3）作为运行的程序；（4）应无明确的结束标志，故选择delta=max1Nj=1j=j+1YjNNYj=1X(1)=(B(1)-A(1,2)*P(2)/A(1,1)NYj=2X(2)=(B(2)-A(2,1)*X(1)-A(2,3:4)*P(3:4)/A(2,2)NYj=N-1X(N-1)=(B(N-1)-A(N-1,N-3:N-2)*X(N-3:N-2) -A(N-1,N)*P(N)/A(N-1,N-1)NYj=NX(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)err=abs(norm(X-P); relerr=err/(norm(X)+eps);P=X;Nerrdelta)|(relerrdeltaYoutputend 实验结果：表1.通过高斯赛德尔迭代得到的x的解1100.463795523816550.53728460519996560.50902292460133060.49822163443617410.49894186023976190.49998535124813080.50008872389013570.5000153188460520.49999479326697530.499997856913467511200.50000010842519920.50000020157668730.5000000226109450.49999998623857220.4999999958739790.50000000053029380.50000000043903880.50000000002393920.4999999999660160.499999999992557321300.50000000000174290.50000000000091690.50.49999999999992050.50.50.4999999999999212 0.50.50000000000091760.500000000001744531400.4999999999925550.4999999999660170.5000000000239390.50000000043903910.50000000053029360.4999999958739790.49999998623857230.50000002261094510.50000020157668730.50000010842519941500.49999785691346750.49999479326697530.5000153188460520.50008872389013570.49998535124813080.49894186023976190.49822163443617410.50902292460133060.53728460519996570.4637955238165501 图11实验结论：X的值如上图所示。附件（代码）：1. P108 1#include#includeusing namespace std;/此函数用于计算矩阵的解float *uptrbk(float *A,int N)float c;float *x=new floatN-1; /生成一维数组xNint n;for(n=1;n=0;n-)AnN=AnN-xn+1*Ann+1;xn=AnN/Ann;return x; /带回计算结果数组xN的地址x/实验的main函数int main()int N;coutN;int i,k;float *uptrbk(float*,int);float *x;float *A=new float*N-1; /生成动态增广矩阵for(i=0;iN;i+)Ai=new floatN;cout请输入增广矩阵的值endl; /输入增广矩阵的值for(i=0;iN;i+)for(k=0;kAik;x=uptrbk(A,N); /计算矩阵的解列向量Xcoutx的值为：endl;for(i=0;iN;i+) /输出矩阵的解的列向量Xcoutxi+1=xiendl;system(pause);return 0;2. P120 1#include#includeusing namespace std;/实验的main函数int main()int N,i,k;float *x;coutN;N=N-1;float *B=new float N;float *lufact(float*,float*,int);float *A=new float*N; /生成动态系数矩阵Afor(i=0;i=N;i+)Ai=new floatN;cout请输入矩阵A的值：endl; /输入系数矩?阵的值for(i=0;i=N;i+)for(k=0;kAik;cout请输入矩阵B的值：endl; /生成动态矩阵Bfor(i=0;iBi; /输入矩阵B的值x=lufact(A,B,N);coutx的值为a：endl;for(i=0;i=N;i+) /输出AX=B解的列向量Xcoutxi+1=xiendl;system(pause);return 0;/使用LU法求解X的函数，A为系数矩阵，B为计算结果，N为矩阵阶数float *lufact(float *A,float *B,int N)int i,j,k;float c;float *U; /生成二维数组Ufloat *x=new float N; /生产保存结果的列矩阵Xfloat *y=new float N; /生产用于保存中间值的列矩阵YU=A;for(i=0;i=N;i+) /将A转换为LUfor(j=i+1;j=N;j+)c=Uji/Uii;for(k=i;k=N;k+)Ujk=Ujk-c*Uik;Uji=c;for(i=0;i=N;i+) /计算中间矩阵Y的值for(j=0;j=0;i-) /计算解X的值for(j=N;ji;j-)yi=yi-xj*Uij;xi=yi/Uii;return x;3. P120 2#include#include#includeusing namespace std;/实验的main函数int main()int N,i;double *x;double *A;double *B;double *buildA(int);double *buildB(int);double *lufact(double*,double*,int);coutN;N=N-1; A=buildA(N);B=buildB(N);x=lufact(A,B,N);coutx的值为：endl;for(i=0;i=N;i+) /输出AX=B解的列矩阵Xcoutxi+1=xiendl;system(pause);return 0;/使用LU法求解X的函数，A为系数矩阵，B为计算结果，N为矩阵阶数double *lufact(double *A,double *B,int N)int i,j,k;double c;double *U; /生成二维数组U，用于储存L&U矩阵double *x=new double N; /生成保存结果的列矩阵Xdouble *y=new double N; /生成用于保存中间值的列矩阵YU=A;for(i=0;i=N;i+) /将A转换为LU矩阵for(j=i+1;j=N;j+)c=Uji/Uii;for(k=i;k=N;k+)Ujk=Ujk-c*Uik;Uji=c;for(i=0;i=N;i+) /计算中间矩阵Y的值for(j=0;j=0;i-) /计算解X的值for(j=N;ji;j-)yi=yi-xj*Uij;xi=yi/Uii;return x;/用于产生阶数为N的矩阵A的函数double *buildA(int N)int i,j;double *A=new double*N; /生成二维动态数组AN-1N-1for(i=0;i=N;i+)Ai=new double N;for(i=0;i=N;i+) /产生矩阵A的元素并存储到AN-1N-1中for(j=0;j=N;j+)Aij=pow(double(i+1),j);return A;/用于产生1*N的列矩阵Bdouble *buildB(int N)int i;double *B=new double N; /生成元素个数为N的动态矩阵BN=N+1;B0=N; /产生矩阵的元素并储存到BN-1中for(i=1;i=N;i+)Bi=(pow(double(i+1),N)-1)/i;return B;4. P120 3#include#include#includeusing namespace std;/实验的main函数int main()int N,i,j;float *LU;float *luchange(float*,int*,int);float *LUsave(float*,int,int,int);coutN;N=N-1;float *A=new float *N;int *B=new int N;float *invA=new float *N; /生成用于存放矩阵A的二维数组for(i=0;i=N;i+)Ai=new float N;cout请输入需要求逆矩阵的矩阵A：endl; /输入矩阵A的值for(i=0;i=N;i+)for(j=0;jAij;for(i=0;i=N;i+) Bi=i;LU=luchange(A,B,N); /将A转换为LU矩阵for(i=0;i=N;i+) /判断A是否存在逆矩阵if(LUii=0) cout矩阵A不存在逆矩阵！endl;system(pause);return 0;for(i=0;i=N;i+) /若存在求解逆矩阵invAi=LUsave(LU,Bi,i,N);coutA的逆矩阵的值为a：endl; /输出矩阵A的逆矩阵for(i=0;i=N;i+)for(j=0;j=N;j+)coutstd:leftsetw(15)invAji;coutendl;system(pause);return 0;/求解LU矩阵的函数，A为系数矩阵，N为矩阵阶数float *luchange(float *A,int *B,int N)int i,j,k,a;float c;float *b=new float N;float *LU; /生成二维数组LU，用于存放L矩阵及U矩阵LU=A;for(i=0;i=N;i+) /将A转换为LUk=i+1;while(LUii=0&k=N)b=LUi;LUi=LUk;LUk=b;a=Bi;Bi=Bk;Bk=a;k+;for(j=i+1;j=N;j+)c=LUji/LUii;for(k=i;k=N;k+)LUjk=LUjk-c*LUik;LUji=c;return LU; /返回LU矩阵的地址/此函数用于求解对应B的解float *LUsave(float *LU,int n,int k,int N)int i,j;float *y=new floatN;float *x=new floatN;float *B=new floatN;for(i=0;i=N;i+)if(i=n)Bi=1;elseBi=0;for(i=0;i=N;i+) /计算中间矩阵Y的值for(j=0;j=0;i-) /计算解X的值for(j=N;ji;j-)yi=yi-xj*LUij;xi=yi/LUii;return x;5.P129 3#include#includeusing namespace std;/定义用于储存数组的类struct typeAdouble *A;int b;/实验的主函数int main()int N,n1,i;double *save(typeA,typeA,int,int);coutN;double *X=new double N-1; /生成用于存放解的数组XtypeA A,B;coutn1; A.A=new double n1;coutA.b;cout请输入系数元素的序列：;for(i=0;iA.Ai;coutB.b;B.A=new doubleB.b;for(i=0;iB.Ai;X=save(A,B,N,n1); /求解方程组cout方程组的解为：endl;for(i=0;iN;i+) /输出方程的解coutxi+1=Xiendl;system(pause);return 0;/用于计算方程组的解double *save(typeA A,typeA B,int N,int n1)int i,j,k,m,deltan;double b;deltan=int(n1-A.b);double *X=new double N-1; /生成用于存储方程组解的数组Xfor(i=0;iN;i+) /初始化数组XXi=0;double x;double delta=1e-6; /定义计算精度double deltax=1;bool run=true; /定于用于判定是否结束计算的变量 while(ru