Question

高斯先列主消元法求解线性方程组AX=b C语言

采用高斯先列主元消元法求解线性方程组AX=b方法说明（以4阶为例）：（1）第1步消元——在增广矩阵（A，b）第一列中找到绝对值最大的元素，将其所在行与第一行交换，再对（A，b）做初等行变换使原方程组转化为如下形式：，注：“*”代表非0。（2）第2步消元——在增广矩阵（A，b）中的第二列中（从第二行开始）找到绝对值最大的元素，将其所在行与第二行交换，再对（A，b）做初等行变换使原方程组转化为：（3）第3步消元——在增广矩阵（A，b）中的第三列中（从第三行开始）找到绝对值最大的元素，将其所在行与第二行交换，再对（A，b）做初等行变换使原方程组转化为：（4）按x4 à x3à x2à x1 的顺序回代求解出方程组的解。要求方程的系数矩阵由用户输入，显示部分中间结果，提供如下菜单：（1）录入系数矩阵；（2）输出结果（显示部分中间结果）；（3）退出。

Answer 1

其中用到了高斯先列主消元法 #include <iostream.h>
#include <stdlib.h>
#include <math.h>

/*楼竞网站www.LouJing.com
拥有该程序的版权,转载请保留该版权.
谢谢合作!*/
double* allocMem(int ); //分配内存空间函数
void GaussLineMain(double*,double*,double*,int );//采用高斯列主元素消去法求解x的初始向量值
void Jacobi(double*,double*,double*,double*,int,int);//利用雅可比迭代公式求解x的值

void main()
{
short matrixNum; //矩阵的行数(列数)
double *matrixA; //矩阵A，初始系数矩阵
double *matrixD; //矩阵D为A中的主对角阵
double *matrixL; //矩阵L为A中的下三角阵
double *matrixU; //矩阵U为A中的上三角阵
double *B; //矩阵B为雅可比方法迭代矩阵
double *f; //矩阵f为中间的过渡的矩阵
double *x; //x为一维数组，存放结果
double *xk; //xk为一维数组，用来在迭代中使用
double *b; //b为一维数组，存放方程组右边系数

int i,j,k;

cout<<"<<请输入矩阵的行数（列数与行数一致）>>：";
cin>>matrixNum;

//分别为A、D、L、U、B、f、x、b分配内存空间
matrixA=allocMem(matrixNum*matrixNum);
matrixD=allocMem(matrixNum*matrixNum);
matrixL=allocMem(matrixNum*matrixNum);
matrixU=allocMem(matrixNum*matrixNum);
B=allocMem(matrixNum*matrixNum);
f=allocMem(matrixNum);
x=allocMem(matrixNum);
xk=allocMem(matrixNum);
b=allocMem(matrixNum);

//输入系数矩阵各元素值
cout<<endl<<endl<<endl<<"<<请输入矩阵中各元素值（为 "<<matrixNum<<"*"<<matrixNum<<"，共计 "<<matrixNum*matrixNum<<" 个元素）"<<">>："<<endl<<endl;
for(i=0;i<matrixNum;i++)
{
cout<<"请输入矩阵中第 "<<i+1<<" 行的 "<<matrixNum<<" 个元素：";
for(j=0;j<matrixNum;j++)
cin>>*(matrixA+i*matrixNum+j);
}

//输入方程组右边系数b的各元素值
cout<<endl<<endl<<endl<<"<<请输入方程组右边系数各元素值，共计 "<<matrixNum<<" 个"<<">>："<<endl<<endl;
for(i=0;i<matrixNum;i++)
cin>>*(b+i);

/* 下面将A分裂为A=D-L-U */
//首先将D、L、U做初始化工作
for(i=0;i<matrixNum;i++)
for(j=0;j<matrixNum;j++)
*(matrixD+i*matrixNum+j)=*(matrixL+i*matrixNum+j)=*(matrixU+i*matrixNum+j)=0;
//D、L、U分别得到A的主对角线、下三角和上三角；其中D取逆矩阵、L和U各元素取相反数
for(i=0;i<matrixNum;i++)
for(j=0;j<matrixNum;j++)
if(i==j&&*(matrixA+i*matrixNum+j)) *(matrixD+i*matrixNum+j)=1/(*(matrixA+i*matrixNum+j));
else if(i>j) *(matrixL+i*matrixNum+j)=-*(matrixA+i*matrixNum+j);
else *(matrixU+i*matrixNum+j)=-*(matrixA+i*matrixNum+j);
//求B矩阵中的元素
for(i=0;i<matrixNum;i++)
for(j=0;j<matrixNum;j++)
{
double temp=0;
for(k=0;k<matrixNum;k++)
temp+=*(matrixD+i*matrixNum+k)*(*(matrixL+k*matrixNum+j)+*(matrixU+k*matrixNum+j));
*(B+i*matrixNum+j)=temp;
}
//求f中的元素
for(i=0;i<matrixNum;i++)
{
double temp=0;
for(j=0;j<matrixNum;j++)
temp+=*(matrixD+i*matrixNum+j)*(*(b+j));
*(f+i)=temp;
}

/* 计算x的初始向量值 */
GaussLineMain(matrixA,x,b,matrixNum);

/* 利用雅可比迭代公式求解xk的值 */
int JacobiTime;
cout<<endl<<endl<<endl<<"<<雅可比迭代开始，请输入希望迭代的次数>>：";
cin>>JacobiTime;
while(JacobiTime<=0)
{
cout<<"迭代次数必须大于0，请重新输入：";
cin>>JacobiTime;
}
Jacobi(x,xk,B,f,matrixNum,JacobiTime);

//输出线性方程组的解 */
cout<<endl<<endl<<endl<<"<<方程组运算结果如下>>"<<endl;
cout.precision(20); //设置输出精度，以此比较不同迭代次数的结果
for(i=0;i<matrixNum;i++)
cout<<"x"<<i+1<<" = "<<*(xk+i)<<endl;

cout<<endl<<endl<<"求解过程结束..."<<endl<<endl;

//释放掉所有动态分配的内存
delete [] matrixA;
delete [] matrixD;
delete [] matrixL;
delete [] matrixU;
delete [] B;
delete [] f;
delete [] x;
delete [] xk;
delete [] b;
}

/*--------------------
分配内存空间函数
www.LouJing.com
--------------------*/
double* allocMem(int num)
{
double *head;
if((head=new double[num])==NULL)
{
cout<<"内存空间分配失败，程序终止运行！"<<endl;
exit(0);
}
return head;
}

/*---------------------------------------------
计算Ax=b中x的初始向量值，采用高斯列主元素消去法
www.LouJing.com
---------------------------------------------*/
void GaussLineMain(double* A,double* x,double* b,int num)
{
int i,j,k;

//共处理num-1行
for(i=0;i<num-1;i++)
{
//首先每列选主元，即最大的一个
double lineMax=fabs(*(A+i*num+i));
int lineI=i;
for(j=i;j<num;j++)
if(fabs(*(A+j*num+i))>fabs(lineMax)) lineI=j;

//主元所在行和当前处理行做行交换，右系数b也随之交换
for(j=i;j<num;j++)
{
//A做交换
lineMax=*(A+i*num+j);
*(A+i*num+j)=*(A+lineI*num+j);
*(A+lineI*num+j)=lineMax;
//b中对应元素做交换
lineMax=*(b+i);
*(b+i)=*(b+lineI);
*(b+lineI)=lineMax;
}

if(*(A+i*num+i)==0) continue; //如果当前主元为0，本次循环结束

//将A变为上三角矩阵，同样b也随之变换
for(j=i+1;j<num;j++)
{
double temp=-*(A+j*num+i)/(*(A+i*num+i));
for(k=i;k<num;k++)
{
*(A+j*num+k)+=temp*(*(A+i*num+k));
}
*(b+j)+=temp*(*(b+i));
}
}

/* 验证Ax=b是否有唯一解，就是验证A的行列式是否为0；
如果|A|!=0，说明有唯一解*/
double determinantA=1;
for(i=0;i<num;i++)
determinantA*=*(A+i*num+i);
if(determinantA==0)
{
cout<<endl<<endl<<"通过计算,矩阵A的行列式为|A|=0，即A没有唯一解。\n程序退出..."<<endl<<endl;
exit(0);
}

/* 从最后一行开始，回代求解x的初始向量值 */
for(i=num-1;i>=0;i--)
{
for(j=num-1;j>i;j--)
*(b+i)-=*(A+i*num+j)*(*(x+j));
*(x+i)=*(b+i)/(*(A+i*num+i));
}
}

/*------------------------------------
利用雅可比迭代公式求解x的精确值
www.LouJing.com
迭代公式为：xk=Bx+f
------------------------------------*/
void Jacobi(double* x,double* xk,double* B,double* f,int num,int time)
{
int t=1,i,j;
while(t<=time)
{
for(i=0;i<num;i++)
{
double temp=0;
for(j=0;j<num;j++)
temp+=*(B+i*num+j)*(*(x+j));
*(xk+i)=temp+*(f+i);
}

//将xk赋值给x，准备下一次迭代
for(i=0;i<num;i++)
*(x+i)=*(xk+i);
t++;
}
}

Answer 2

/*高斯消元法及矩阵乘法计求解n*n阶行列式不等于零的一次线性方程*林木*C语言程序经典实现 */
/* 转载自： http://user.qzone.qq.com/464372077/infocenter?ptlang=2052&ADUIN=464372077&ADSESSION=1271069269&ADTAG=CLIENT.QQ.2719_Mysrv.0 */
#include <stdio.h>
#include <stdlib.h>
void fun1(float *,int);
void fun2(float *,int);
void fun3(float *,int);
float fun4(float *,int);
void fun5(float *,float *,int);
void fun6(float *,float *,float *,int);
int jiemian();
void main()
{
int i,j,n,e;
float *a,*b,*c,*d,answer;
do
{
e=jiemian();
if(e==1)
{
for(i=0;i<3;i++)
printf("<----------------------------------------------------------------------------->\n");
printf(" 欢迎进入*雒森*n*n阶一次线性方程求解程序,请准确输入以便我为您正确计算! \n");
for(i=0;i<3;i++)
printf("<----------------------------------------------------------------------------->\n");
printf("\n\n请输入方程的阶数n(n>1):");
do
{
scanf("%d",&n);
if(n<2)
printf("\n\n请重新输入正确的方程的阶数n(n>1):");
}while(n<2); a=(float *)malloc(n*n*2*sizeof(float));
c=(float *)malloc(n*n*sizeof(float));
b=(float *)malloc(n*n*sizeof(float));
d=(float *)malloc(n*sizeof(float));
for(i=0;i<n;i++)
{
printf("请输入您要计算的方程的系数矩阵的:\n");
for(j=0;j<2*n;j++)
{
if(j<n)
{
printf("第%d行%d列:",i+1,j+1);
scanf("%f",&a[i*2*n+j]);
c[i*n+j]=a[i*2*n+j];
}
else if(j==n+i)
a[i*2*n+j]=1;
else if(j>n&&j!=n+i)
a[i*2*n+j]=0;
}
}
for(i=0;i<n;i++)
{
printf("\n请输入线性方程常数,列矩阵的第%d行第1列:",i+1);
scanf("%f",&d[i]);
}
fun3(c,n);
answer=fun4(c,n);
free(c);
printf("您输入的系数矩阵的行列式的值为:%f!\n",answer);
if(fabs(answer)+0.1!=0.1)
{
fun1(a,n);
for(i=0;i<n;i++)
for(j=0;j<n;j++)
b[i*n+j]=a[(i*2+1)*n+j];
free(a);
fun2(b,n);
fun5(b,d,n);
}
else
printf("\n如果您输入的方程的系数矩阵的行列式值等于0,则是是奇异矩阵,是不可逆的!\n方程的解不唯一该程序无法计算!抱歉！\n\n\n");
}
}while(e!=3);
}
void fun1(float a[],int n)
{
int i,j,d,f;
float m,w;
for(j=0;j<n;j++)
{
for(i=0;i<n;i++)
{
if(a[j*2*n+j]==0)
{
for(d=j+1;d<n;d++)
{
if(a[d*2*n+j]!=0)
break;
}
for(f=j;f<(2*n);f++)
{
m=a[j*2*n+f];
a[j*2*n+f]=a[d*2*n+f];
a[d*2*n+f]=m;
}
}
if(i==0&&a[j*2*n+j]!=1)
{
w=a[j*2*n+j];
for(f=j;f<2*n;f++)
a[j*2*n+f]/=w;
}
if(i!=j)
{
w=(-a[i*2*n+j]);
for(f=j;f<(2*n);f++)
a[i*2*n+f]+=a[j*2*n+f]*w;
}
}
}
}
void fun2(float b[],int n)
{
int i,j;
printf("您所要求的方程的系数矩阵的逆矩阵为:\n");
for(i=0;i<n;i++)
{
for(j=0;j<n;j++)
printf("b[%d行%d列]=%f",i+1,j+1,b[i*n+j]);
printf("\n\n\n");
}
}
void fun3(float c[],int n)
{
int i,j,d,f;
float m=1,w;
for(j=0;j<n;j++)
{
for(i=0;i<n;i++)
{
if(c[j*n+j]==0)
{
for(d=j+1;d<n;d++)
{
if(c[d*n+j]!=0)
break;
}
for(f=j;f<n;f++)
c[j*n+f]+=c[d*n+f]*m;
}
if(i!=j)
{
w=-(c[i*n+j]/c[j*n+j]);
for(f=j;f<n;f++)
c[i*n+f]+=c[j*n+f]*w;
}
}
}
}
float fun4(float c[],int n)
{
int i;
float answer=1;
for(i=0;i<n;i++)
answer*=c[i*n+i];
return answer;
}
void fun5(float b[],float d[],int n)
{
int i;
static float *x;
x=(float *)malloc(n*sizeof(float));
fun6(b,d,x,n);
i=jiemian();
if(i==2)
{
printf("你要解得线性方程的解为:\n");
for(i=0;i<n;i++)
printf("x[%d]=%f\n",i+1,x[i]);
free(x);free(b);free(d);
for(i=0;i<3;i++)
printf("<----------------------------------------------------------------------------->\n");
printf(" *雒森*n*n阶线性方程求解程序计算完毕 \n");
for(i=0;i<3;i++)
printf("<----------------------------------------------------------------------------->\n");
}
}
void fun6(float b[],float d[],float x[],int n)
{
int i,j,k;
for(i=0;i<n;i++)
{
x[i]=0;
for(k=0;k<n;k++)
x[i]+=b[i*n+k]*d[k];
}
}
int jiemian()
{
int a;
do
{
printf("请选择您要的功能：\n");
printf("1代表计算线性方程录入数据，2代表输出结果，3代表退出");
printf("请输入：");
scanf("%d",&a);
if(a==3)
exit(1);
else if(a==1)
return(1);
else if(a==2)
return(2);
else
printf("您输入错误请重新输入：\n");
}while(a<1||a>3);
}

Answer 3