数值分析计算实习题列主元高斯消去法解线性方程组
数值分析计算实习题
第5章解线性方程组的直接方法
【选题
列主元高斯消去法解线性方程组。
书上的计算实习题1、2、3都要求用列主元高斯消去法解线性方程组，所以考虑写一个普适的程序来实现。
对于线性方程组Ax二b,程序允许用户从文件读入矩阵数据或直接在屏幕输入数据。
文件输入格式要求:
(1)第一行为一个整数n (2<=n<=100)，表示矩阵阶数。
(2)第2~n+l行为矩阵A各行列的值。
(3)第n+2~n+n+2行为矩阵b各行的值。
屏幕输入：按提示输入各个数据。
输出：A. b、det(A).列主元高斯消去计算过程、解向量X。
【算法说明】
设有线性方程组Ax=b,其中设A为非奇异矩阵。方程组的增广矩阵为
?12
?21
[Nb] =
第1步（k=l ）：首先在A的第一列中选取绝对值最大的元素?I,作为第一步的主元素:
?|| H0
然后交换（A, b）的第1行与第I行元素，再进行消元计算。
设列主元素消去法已经完成第1步到第k?l步的按列选主元，交换两行，消元计算得到与原方程组等价的方程组 A（k）x=b（k）
4? …4;)
…唸)
?
忒
?
?
輕
■
[A.b]T[A ⑹,b")] = ??■
咲■
■
■
■
■
* *■
〃伏)
?? - %
■
第k步计算如下: 对于 k=l, 2, ?…，0-1
（1）按列选主元：即确定t使
（2）如果tHk,则交换［A, b］第t行与第k行元素。（3）消元计算
5
4* J 叫=一鱼(=^ + 1,…,H)
% 吗 <-?y + 〃如伽 (fJ = R + l,…/)
b- <-勺+加汝仇， (i = /c + l,…,《)
消元乘数mik 满足:
n (%-D 内)
X1 < ------ -- ---- 9(j = ? 一 1,?一2■…J)tk M 1,(,=斤 +1, ???,?)
fet e
(4)回代求解
【程序】
/*
【普适列主元消去法解线性方程组】对于线性方程组：Ax=b 输入：
［选择屏幕直接输入］
1.A的行阶数n(l <= 11<= 100)
2.A的值
3.b的值
［选择读取文件1
文件名(和主程序同级文件夹下)
输出：
I.A
2?b
3? det(A)
4解向疑X */
#inciude
#include
#include
double det A:
int n.mark = 1;
〃读入数据
void input(){
int ij:
char ch[20],name[ 100];
HLE *f;
printf("\iv-An是否从文件读取数据(Y/N):”); scanf(”％st&ch);
if(ch[0] = Y II ch[0] = y)( prinif(“请输入文件名(包括扩展需)：")； scanf("%s".name):
f = fop cn(namc「T');
fscanfCf/'%d'\&n);
ford = 0:i < n;i 卄) for(j = 0;j < nJ ++) fscanf(f,'*%ir\&A[i]|j))：for(i = 0:i < n;i
卄)
fscanf(「％F?&b[i]);
else{
prin氓”请输入A的阶数：”)；scanf( '%d %d\&n)： prinifC请输入A的值：”)；
for(i = 0:i V n:i ++)
for(j = 0;j < n:j ++) scanf("%lf\&A[i]U]);
phnifC请输入b的值：”)；
for(i = 0:i < n;i 卄)
scanf(''%lf\&b[i])：
〃讣算行列式的值
double det{double s[105](105] jni m){
int z.jk
double b[105][105Klotal = 0.r; /*b[Nl[N]用于存放，在矩阵s[Nl[N冲元素s[0]的余子式灯if(m>2){
for(z = 0：z V m:z++){
for(j = 0;j V m ? 1 j ++)
for(k = 0:k vn卜l;k ++)
if(k >= z)
bUKk] = sU+l](k+l];
else
bLi][k] = s[j+l][k];
if(z % 2==0)
r=s[0)⑵ * dcKb.m - 1)：/*递归调用制
else
r= (-1) * s[0](z] * det(b.m -1);
total = total + r;
else if(m == 2)
total = s[0][0] *s(l][l]-s[0](l] *s[l][01：
else if(m == 1)
total =
s[0][0];
return total;
// 输出A^llb和dcl(A)
void ouipui_l(){ int i j;
primlTAW);
for(i = 0;i < n:i ++){ for(j = 0:j < n:j ++)
p rintf('-%15.4f\A[i]lj]);
prinif(W);
prinlf(5b = \rf);
for{i = 0;i < n:i ++)
prinif("%15.4l\iV\b[i]):
printf("\ndet(A) = %?4f\n”?dclA);
//主il?算函数void couni_x(){
int ij,k:
int max; double tmpjnik;
〃构造增广矩阵
for{i = 0;i for。= 0:j < n;j ++) A_B[i]Ul=A[i][jl； A_B[i](n] = b[ij: prinlHAn展示消元过程的增广矩阵):\j门; prinlf("A_B=\n"'); for(i = 0;i < n:i ++){ for(j = 0:j <= n:j ++) p rimf(”％15?4代 prinif(W); 〃按列选主元 for(k = 0:k max = k: for(i = k:i V n;i ++) if(abs(A-B[i][k]) > abs(A_B(max][k])) max = i; 〃若A|maxl[k]=O.即dcl(A) = O?则无唯一解.退出il?算if(A_B(maxl[k] = 0){ mark = 0; break: 〃换行ifCmax != k)( for(j = k:j <= n;j ++)( tnip = A_B[k]ljl： A_B[k]U]=A_B[max]U]; A_B[max][j] = tmp; 〃消元il?算for(i = k + l;i vn:i ++){ A_BEi][kl = A_B[i](kl/A_B[k]rk]; for(j = k + 1 j V n;j ++) A_B(iIU] = A_B[i][j]? A_B[il(k] * A_B[k][j]： A_B[i][n] = A_B[i](nl - A_B[i][kl * A_B(k][nJ; prinif{"A_B = \n")： for(i = 0:i < n;i ++){ for。= 0:j <= n;j ++) p rinlfW); //若A[n]In] = OJiPdct(A) = 0,则无唯一解，退出计算if(A_B[n - ll[n-1] = 0) mark = 0; else{ x[n-1] = A_B[n - l][n] = A_B[n- ll[n]/A_B[n- l][n-1]: for(i = n- 2:i >= 0;i -){ x(il=A_B[i](n]=A_B[i](n]; for(j = i + l;j x[i]= -= A_B[i]U] * A_BU][n]；