1、数值方法实验报告1线性方程组 AX=B的数值计算方法实验【摘要】在自然科学与工程技术中很多问题的解决常常归结为解线性代数方程组。例如电学中的网络问题，船体数学放样中建立三次样条函数问题，用最小二乘法求实验数据的曲线拟合问题，解非线性方程组的问题，用差分法或者有限元法解常微分方程，偏微分方程边值问题等都导致求解线性方程组。线性代数方面的计算方法就是研究求解线性方程组的一些数值解法与研究计算矩阵的特征值及特征向量的数值方法。关于线性方程组的数值解法一般有两类：直接法和迭代法。关 键 字 高斯消元法、三角分解法、高斯-赛德尔迭代、稀疏矩阵1、实验目的1.掌握高斯消元法、三角分解法、高斯赛德尔迭代发的
2、编程技巧。2.掌握线性方程组 AX=B 的数值计算方法。3.掌握矩阵的基本编程技巧。2、实验原理1.高斯消元法数学上，高斯消元法是线性代数规划中的一个算法，可用来为线性方程组求解。高斯（Gauss）夏鸥按法其实是将一般的线性方程组变换为三角形（上三角）方程组求解问题（消元法），只是步骤规范，便于编写计算机程序。一般高斯消元法包括两过程：先把方程组化为同解的上三角形方程组，再按相反顺序求解上三角方程组。前者称为消去或消元过程，后者称回代过程。消去过程实际上是对增广矩阵作行初等变换。对一般的 n 阶方程组，消去过程分 n-1 步：第一步消去 下方元素。第二1a步消去 下方元素， ，第 n-1 步消
3、去 下方元素。即第 k 步将第 k 行2a 1-na，的适当倍数加于其后各行，或可说是从 k+1n 行减去第 k 行的适当倍数，使它们第 k 列元素变为零，而其余列元素减去第 k 行对应列元素的倍数。2.三角分解法数值方法实验报告2三角分解法是将原正方 (square)矩阵分解成一个上三角形矩阵或是排列(permuted) 的上三角形矩阵和一个 下三角形矩阵，这样的分解法又称为 LU 分解法。它的用途主要在简化一个大矩阵的行列式值的计算过程，求 反矩阵，和求解联立方程组。不过要注意这种分解法所得到的上下三角形矩阵并非唯一，还可找到数个不同 的一对上下三角形矩阵，此两三角形矩阵相乘也会得到原矩阵
4、。3.高斯 赛德尔迭代高斯赛德尔迭代（Gauss Seidelmethod）是数值线性代数中的一个迭代法，可用来求出线性方程组解的近似值。研究雅可比迭代法，我们发现在逐个求 的分量时，当计算到 时,）（ 1kX ）（ 1kiX分量 ， ， 都已经求得，而仍用旧分量 ， 计算 。）（ 1kX）（ 1k-i k1）（ 1-i）（ i由于新计算出的分量比旧分量准确些，因此设想一旦新分量 ，）（ k求出，马上就用新分量 ， 代替雅可比迭代法中 ，）（ 1k-i ）（ 1kX）（ 1k-i k1X， 来求 这就是高斯-赛德尔(Gauss-Seidel) 迭代法。）（ 1-iX）（ 1ki把矩阵 A 分解
5、成 ULDA(6) 其中 na,.diagD21, 分别为 的主对角元除外的下三角和上三角部分,于是,方程组(1) 便可以写成bx即 2fB其中 bLD,ULD112 (7)以 2B为迭代矩阵构成的迭代法(公式) 221fxkk(8)称为高斯塞德尔迭代法(公式),用变量表示的形式为 ,.k,ni xabaxij nij)k(j)k(ji)k(i 10111(9)4.稀疏矩阵数值方法实验报告3矩阵中非零元素的个数远远小于矩阵元素的总数，并且非零元素的分布没有规律，则称该矩阵为稀疏矩阵(sparse matrix)；与之相区别的是，如果非零元素的分布存在规律（如上三角矩阵、下三角矩阵、对角矩阵），
6、则称该矩阵为特殊矩阵。常见于进行大量数据计算。3、实验内容1.P108 1许多科学应用包含的矩阵带有很多的零。在实际情况中很重要的三角形线性方程组有如下形式：d1x1+c1x2 =b1a1x1+d2x2+c2x3 =b2a2x2+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=35721061234使用三角分解法求解 X。3.P120 2求解线性方程组 AX=B，其
7、中 A=aijNN，,a ij=ij-1；而且 B=bijN1，b 11=N，当 i时，b ij=(iN-1)/(i-1)。对 N=3，7， 11 的情况分别求解。精确解为 X=1 1 1 1。对得到的结果与精确解的差异进行解释。4.P120 3数值方法实验报告4通过重复求解 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 =b1a1x1+d2x2+c2x3 =b2a2x2+d3x3+c3x4
8、 =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 (1)kj()()()()111kkkkjjjj jjb ax式（2）数值方法实验报告5= ，j=1 ，2，N (1)kjx(1)(!
9、)()()11kkkkjj jj jNjbaxaxax式（3）（ii）根据（i）中设计的算法构建一个 C+程序，并求解下列上三角线性方程组（a）4m 1+ m2 =3 （b）4m 1+ m2 =1m1+4m2+m3=3 m1+4m2+m3=2m2+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 =5x1-2x2+12x3-2x4+x5=5x2-
10、2x3+12x4-2x5+x6=5x46-2x47+12x48-2x49+x50=5x47-2x48+12x49-2x50=5x48-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）的结果为 。134x1234x数值方法实验报告6算法流程图：startInputA,B,P,de
11、lta,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)j=j+1YNYNNNYYk=k+1数值方法实验报告7实验结果：结果截图(1) errmax1jNj=1YNYk=k+1j=j+1数值方法实验报告17实验结果：表 1.通过高斯赛德尔迭代得到的 x 的解x1100.463795523816550.53728460519996560.50902292460133060.49822163443617410.49894186023976190.49
12、998535124813080.50008872389013570.5000153188460520.49999479326697530.499997856913467511 0.50 0.50 0.50 0.49 0.49 0.50 0.50 0.50 0.49 0.49j=1X(1)=(B(1)-A(1,2)*P(2)/A(1,1)j=NX(N)=(B(N)-A(N,N-1)*(X(N-1)/A(N,N)err=abs(norm(X-P); relerr=err/(norm(X)+eps);P=X;err#includeusing namespace std;/此函数用于计算矩阵的解flo
13、at *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;数值方法实验报告21float *uptrbk(float*,int);float *x;float *A=new float*N-1; /生成动态增广矩阵for(i=0;iAik;x=uptrbk(A,N); /计算矩阵的解列向量Xc
14、out#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;iAik;coutBi; /输入矩阵B的值x=lufact(A,B,N);cout=0;i-) /计算解X的值for(j=N;ji;j-)yi=yi-xj*Uij;xi=yi/Uii;return x;数值方法实验报告233.P120 2#incl
15、ude#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);cout=0;i-) /计算解X的值for(j=N;ji;j-)yi=yi-xj*Uij;xi=yi/Uii;return x;/用于产生阶数为N的矩
16、阵A的函数double *buildA(int N)int i,j;double *A=new double*N; /生成二维动态数组AN-1N-1for(i=0;i#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 *
17、N; /生成用于存放矩阵A的二维数组for(i=0;iAij;for(i=0;i=0;i-) /计算解X的值for(j=N;ji;j-)yi=yi-xj*LUij;xi=yi/LUii;return x;5.P129 3#include#includeusing namespace std;数值方法实验报告28/定义用于储存数组的类struct typeAdouble *A;int b;/实验的主函数int main()int N,n1,i;double *save(typeA,typeA,int,int);coutN;double *X=new double N-1; /生成用于存放解的数组
18、XtypeA A,B;coutn1;A.A=new double n1;coutA.b;coutA.Ai;coutB.b;B.A=new doubleB.b;for(i=0;iB.Ai;X=save(A,B,N,n1); /求解方程组cout=0)b=b-Xj*A.Am;m+;m+;for(j=i+1;j=B.b)k=0;deltax=abs(XN-1-x); /计算两次计算的相对误差return X;数值方法实验报告306.P129 4#include#includeusing namespace std;/定义用于储存数组的类struct typeAdouble *A;int b;/实验的
19、主函数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;coutA.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+) /输出方程的解cout“x“i+1“=“Xiendl;system(“pause“);return 0;/用于计算方程组的解
