数值分析计算实习题列主元高斯消去法解线性方程组

数值分析计算实习题第5 章 解线性方程组的直接方法【选题】 列主元高斯消去法解线性方程组。 书上的计算实习题1、2、3都要求用列主元高斯消去法解线性方程组，所以考虑写一个普适的程序来实现。 对于线性方程组Ax=b，程序允许用户从文件读入矩阵数据或直接在屏幕输入数据。 文件输入格式要求：（1）第一行为一个整数n（2<= n <= 100）,表示矩阵阶数。（2）第2 n+1行为矩阵A各行列的值。（3）第n+2 n+n+2行为矩阵b各行的值。屏幕输入：按提示输入各个数据。 输出：A、b、det(A)、列主元高斯消去计算过程、解向量x。【算法说明】 设有线性方程组Ax=b，其中设A为非奇异矩阵。方程组的增广矩阵为 第1步（k=1）：首先在A的第一列中选取绝对值最大的元素，作为第一步的主元素: 然后交换（A，b）的第1行与第l行元素，再进行消元计算。 设列主元素消去法已经完成第1步到第k-1步的按列选主元，交换两行，消元计算得到与原方程组等价的方程组 A(k)x=b(k) 第k步计算如下： 对于k=1，2，n-1 （1）按列选主元：即确定t使 （2）如果tk，则交换A，b第t行与第k行元素。 （3）消元计算 消元乘数mik满足: （4）回代求解【程序】/* 【普适列主元消去法解线性方程组】 对于线性方程组：Ax=b 输入： 选择屏幕直接输入 1. A的行阶数n(1 <= n <= 100) 2. A的值 3. b的值 选择读取文件 文件名（和主程序同级文件夹下） 输出： 1. A 2. b 3. det(A) 4. 解向量x */#include <stdio.h>#include <stdlib.h>#include <math.h>double A105105,A_B105105,b105,x105;double detA;int n,mark = 1;/读入数据void input() int i,j; char ch20,name100; FILE *f; printf("n-n是否从文件读取数据（Y/N）:"); scanf("%s",&ch); if(ch0 = 'Y' | ch0 = 'y') printf("请输入文件名（包括扩展名）："); scanf("%s",name); f = fopen(name,"r"); fscanf(f,"%d",&n); for(i = 0;i < n;i +) for(j = 0;j < n;j +) fscanf(f,"%lf",&Aij); for(i = 0;i < n;i +) fscanf(f,"%lf",&bi); else printf("请输入A的阶数："); scanf("%d %d",&n); printf("请输入A的值："); for(i = 0;i < n;i +) for(j = 0;j < n;j +) scanf("%lf",&Aij); printf("请输入b的值："); for(i = 0;i < n;i +) scanf("%lf",&bi); /计算行列式的值double det(double s105105,int m) int z,j,k; double b105105,total = 0,r; /*bNN用于存放，在矩阵sNN中元素s0的余子式*/ if(m > 2) for(z = 0;z < m;z+) for(j = 0;j < m - 1;j +) for(k = 0;k < m-1;k +) if(k >= z) bjk = sj+1k+1; else bjk = sj+1k; if(z % 2=0) r = s0z * det(b,m - 1); /*递归调用*/ else r = (-1) * s0z * det(b,m - 1); total = total + r; else if(m = 2) total = s00 * s11 - s01 * s10; else if(m = 1) total = s00; return total; / 输出A和b和det(A) void output_1() int i,j; printf("A = n"); for(i = 0;i < n;i +) for(j = 0;j < n;j +) printf("%15.4f",Aij); printf("n"); printf("nb = n"); for(i = 0;i < n;i +) printf("%15.4fn",bi); printf("ndet(A) = %.4fn",detA); / 主计算函数void count_x() int i,j,k; int max; double tmp,mik; /构造增广矩阵 for(i = 0;i < n;i +) for(j = 0;j < n;j +) A_Bij = Aij; A_Bin = bi; printf("n展示消元过程(A_B为A的增广矩阵):n"); printf("A_B = n"); for(i = 0;i < n;i +) for(j = 0;j <= n;j +) printf("%15.4f",A_Bij); printf("n"); /按列选主元 for(k = 0;k < n;k +) max = k; for(i = k;i < n;i +) if(abs(A_Bik) > abs(A_Bmaxk) max = i; /若Amaxk = 0,即det(A) = 0,则无唯一解，退出计算 if(A_Bmaxk = 0) mark = 0; break; /换行 if(max != k) for(j = k;j <= n;j +) tmp = A_Bkj; A_Bkj = A_Bmaxj; A_Bmaxj = tmp; /消元计算 for(i = k + 1;i < n;i +) A_Bik = A_Bik / A_Bkk; for(j = k + 1;j < n;j +) A_Bij = A_Bij - A_Bik * A_Bkj; A_Bin = A_Bin - A_Bik * A_Bkn; A_Bik = 0; printf("A_B = n"); for(i = 0;i < n;i +) for(j = 0;j <= n;j +) printf("%15.4f",A_Bij); printf("n"); / 若Ann = 0,即det(A) = 0,则无唯一解，退出计算 if(A_Bn - 1n - 1 = 0) mark = 0; else xn - 1 = A_Bn - 1n = A_Bn - 1n / A_Bn - 1n - 1; for(i = n - 2;i >= 0;i -) xi = A_Bin = A_Bin; for(j = i + 1;j < n;j +) xi = A_Bin -= A_Bij * A_Bjn; xi = A_Bin /= A_Bii; / 输出解向量x void output_2() int i,j; if(!mark) printf("n该方程没有唯一解！n"); else printf("nb = n"); for(i = 0;i < n;i +) printf("%15.5fn",xi); int main() char ch20; printf("普适列主元消去法解线性方程组n"); while(1) input(); / 读取数据 detA = det(A,n); / 计算A的行列式 output_1(); / 输出A和b和det(A) count_x(); / 主计算函数 output_2(); / 输出解向量x printf("nn是否继续（Y/N）:"); scanf("%s",&ch); if(ch0 = 'Y' | ch0 = 'y') continue; else break; 【结果测试】一共采用了书上的例题和习题内的5组数据（d1d5），均采用文件读入方式读取，计算结果，与matlab运算结果进行比对均相符合，基本验证了程序的正确性。下面给出各组数据和运行结果：数据3311 1-3111-3-410 -7 0 1-3 2.099999 6 25 -1 5 -12 1 0 2851-410 7 8 77 5 6 58 6 10 97 5 9 1032233331-运行结果普适列主元消去法解线性方程组-是否从文件读取数据（Y/N）:yA =b = 1.0000展示消元过程(A_B为A的增广矩阵):A_B =A_B = 3.0100 A_B =A_B =b =是否继续（Y/N）:y-是否从文件读取数据（Y/N）:yA = 0.9900 b =展示消元过程(A_B为A的增广矩阵):A_B =A_B =A_B =A_B =b =是否继续（Y/N）:y-是否从文件读取数据（Y/N）:yA =b =展示消元过程(A_B为A的增广矩阵):A_B =A_B =A_B =A_B =b =是否继续（Y/N）:y-是否从文件读取数据（Y/N）:yA =b =展示消元过程(A_B为A的增广矩阵):A_B =A_B =A_B =A_B =A_B =b =是否继续（Y/N）:y-是否从文件读取数据（Y/N）:yA =b =展示消元过程(A_B为A的增广矩阵):A_B =A_B =A_B =A_B =A_B =b =是否继续（Y/N）: