高斯消元法c语言
‘壹’ 高斯先列主消元法求解线性方程组AX=b c语言
其中用到了高斯先列主消元法 #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++;
}
}
‘贰’ 【编程求助】用c语言或者c++编程,实现用高斯消元法求解线性方程组Ax=b。
void gaussj(double a[], int n, double b[])
{
int i,j,k,l,ll,irow,icol;
double big,pivinv,m;
int ipiv[50], indxr[50], indxc[50];
for (j=0;j<=n-1;j++)
{
ipiv[j]=0;
}
for (i=0;i<=n-1;i++)
{
big=0.0;
for (j=0;j<=n-1;j++)
{
if(ipiv[j]!=1)
{
for(k=0;k<=n-1;k++)
{
if(ipiv[k]==0)
{
if(fabs(a[j*n+k])>=big)
{
big=fabs(a[j*n+k]);
irow=j;
icol=k;
}
else if(ipiv[k]>1)
{
cout<<"singular matrix";
}
}
}
}
}
ipiv[icol]=ipiv[icol]+1;
if(irow!=icol)
{
for(l=0;l<=n-1;l++)
{
m=(a[irow*n+l]);
a[irow*n+l]=a[icol*n+l];
a[icol*n+l]=m;
}
m=b[irow];
b[irow]=b[icol];
b[icol]=m;
}
indxr[i]=irow;
indxc[i]=icol;
if(a[icol*n+icol]==0.0)
{
cout<< "singular matrix.";
}
pivinv=1.0/(a[icol*n+icol]);
a[icol*n+icol]=1.0;
for(l=0;l<=n-1;l++)
{
a[icol*n+l]=a[icol*n+l]*pivinv;
}
b[icol]=b[icol]*pivinv;
for(ll=0;ll<=n-1;ll++)
{
if(ll!=icol)
{
m=a[ll*n+icol];
a[ll*n+icol]=0.0;
for(l=0;l<=n-1;l++)
{
a[ll*n+l]=a[ll*n+l]-a[icol*n+l]*m;
}
b[ll]=b[ll]-b[icol]*m;
}
}
}
for(l=n-1;l<=0;l--)
{
if(indxr[l]!=indxc[l])
{
for(k=0;k<=n-1;k++)
{
m=a[k*n+indxr[l]];
a[k*n+indxr[l]]=a[k*n+indxc[l]];
a[k*n+indxr[l]]=m;
}
}
}
}
‘叁’ 采用高斯先列主元消元法求解线性方程组AX=b,编写一个程序C语言,急需
这个程序我做过的。LZ检验下: // 高斯消元求矩阵逆。
#include<stdio.h>
#include<math.h>#define N 100//定义矩阵的最大行int n;//表示矩阵的行,列。
double matix[N][N];//矩阵的最大行,最大列不
double unit[N][N];bool findmax(int s)//从s到n行选择最大的,作为主元。
{
int i,j,k;
double mas,temp;
mas=fabs(matix[s][s]);k=s;
for(i=s+1;i<n;i++)
{
if(mas<fabs(matix[i][s])) {mas=fabs(matix[i][s]);k=i;}
}
if(mas==0) return false;
//交换两行
for(j=0;j<n;j++)
{
temp=matix[s][j];
matix[s][j]=matix[k][j];
matix[k][j]=temp;
temp=unit[s][j];
unit[s][j]=unit[k][j];
unit[k][j]=temp;
}
return true;
}void chuli(int s)
{
int i,j;
double mas=matix[s][s],r;
for(i=s+1;i<n;i++)
{
r=matix[i][s]/mas;
for(j=0;j<n;j++)
{
matix[i][j]=matix[i][j]-matix[s][j]*r;
unit[i][j]=unit[i][j]-unit[s][j]*r;
}
}
}void solve(int s)
{
int i,j;
double mas;
mas=matix[s][s];
for(i=s;i<n;i++)
matix[s][i]=matix[s][i]/mas;
for(i=0;i<n;i++)
unit[s][i]=unit[s][i]/mas;
for(i=s-1;i>=0;i--)
{
mas=matix[i][s];
matix[i][s]=0;
for(j=0;j<n;j++)
unit[i][j]=unit[i][j]-mas*unit[s][j];
}
}
int main()
{
int i,j;
//输入
scanf("%d",&n);
for(i=0;i<n;i++)
{
for(j=0;j<n;j++)
{
if(i==j) unit[i][j]=1;
else unit[i][j]=0;
scanf("%lf",&matix[i][j]);
}
}
//输出
//转成上三角
for(i=0;i<n-1;i++)
{
//每一次选一个绝对值最大的值作为aii
if(findmax(i))
{
chuli(i);
}
else
{
printf("该矩阵不可逆");
return 0;
}
}
//转成单位阵
for(i=n-1;i>=0;i--)
{
solve(i);
}
//其逆矩阵是
for(i=0;i<n;i++)
{
for(j=0;j<n;j++)
printf("%.2lf ",unit[i][j]);
printf("\n");
}
return 0;
}
‘肆’ C语言用高斯消元法解n元线性方程
#include<math.h>
#include<stdio.h>
#include<stdlib.h>
#include<conio.h>
#define NUMBER 20
#define Esc 0x1b
#define Enter 0x0d
float A[NUMBER][NUMBER+1] ,ark;
int flag,n;
void exchange(int r,int k);
float max(int k);
void message();
int main()
{
float x[NUMBER]; /*此数组用于存放方程解*/
int r,k,i,j;
char celect;
system("cls");
printf("\n\n用Gauss列主元消元法解线性方程组");
printf("\n\n1.解方程组请按Enter.");
printf("\n\n2.退出程式请按Esc.");
celect=getch();
if(celect==Esc)
exit(0);
printf("\n\n 输入方程组的维数:n=");
scanf("%d",&n);
printf(" \n\n现在输入系数矩阵A和向量b:");
for(i=1;i<=n;i++)
{
printf("\n\n请输入a%d1--a%d%d系数和向量b%d:",i,i,n,i);
/*实现将每一行中的系数和向量一次性输入,数之间用空格格开,输完后回车确定*/
for(j=1;j<=n+1;j++) /*将刚才输入的数存入数组*/
scanf("%f",&A[i][j]);
}
for(k=1;k<=n-1;k++)
{
ark=max(k);
if(ark==0) /*判断方程是否为线性方程,即是否合法*/
{
printf("\n\n此方程组不合法!");message();
}
else if(flag!=k)
exchange(flag,k);
for(i=k+1;i<=n;i++)
for(j=k+1;j<=n+1;j++)
A[i][j]=A[i][j]-A[k][j]*A[i][k]/A[k][k];
}
x[n]=A[n][n+1]/A[n][n];
for( k=n-1;k>=1;k--)
{
float me=0;
for(j=k+1;j<=n;j++)
{
me=me+A[k][j]*x[j];
}
x[k]=(A[k][n+1]-me)/A[k][k];
}
for(i=1;i<=n;i++)
{
printf(" \n\nx%d=%f",i,x[i]);
}
message();
return 1;
}
void exchange(int r,int k) /*交换行的矩函数*/
{
int i;
for(i=1;i<=n+1;i++)
A[0][i]=A[r][i];
for(i=1;i<=n+1;i++)
A[r][i]=A[k][i];
for(i=1;i<=n+1;i++)
A[k][i]=A[0][i];
}
float max(int k) /*比校系数大小的函数*/
{
int i;
float temp=0;
for(i=k;i<=n;i++)
if(fabs(A[i][k])>temp)
{
temp=fabs(A[i][k]);
flag=i;
}
return temp;
}
void message() /*实现菜单选择的函数*/
{
printf("\n\n 继续运算按 Enter ,退出程式按 Esc!");
switch(getch())
{
case Enter: main();
case Esc: exit(0);
default:{printf("\n\n不合法的输入!");message();}
}
}
‘伍’ 用C语言编程高斯全主元消元法
//TurboC 2.0太落后了,建议使用VC++6.0。
#include"stdio.h"
#include"math.h"
//最大49阶
#define N 50
void Gauss(float U[N][N],int n);
void main()
{
int n,i,j;
float U[N][N];
printf("------------特殊说明---------------\n");
printf("当输出的数据含有<-1.#IND>时,表示在计算过程中数据已经出现溢出!\n");
printf("-----------------------------------\n");
printf("输入对应方程的阶数:");
scanf("%d",&n);
for(i=0;i<N;i++)
for(j=0;j<N;j++)
U[i][j]=0;
printf("输入方程组的增广矩阵:\n");
for(i=0;i<n;i++)
for(j=0;j<=n;j++)
scanf("%f",&U[i][j]);
Gauss(U,n);
}
//高斯选列主元消去法
void Gauss(float U[N][N],int n)
{
int i,j,m,row;
float max,t,sum;
float result[50];
for(m=0;m<n-1;m++)
{
//选取主元
max=U[m][m];
for(i=m;i<n;i++)
{
if(fabs(max)<fabs(U[i][m]))
{
max=U[i][m];
row=i;
}
}
if(fabs(max)<0.01)
{
printf("主元接近于零,方法失效!\n");
return;
}
else
{
if(max!=U[m][m])
{
for(j=m;j<=n;j++)
{
t=U[m][j];
U[m][j]=U[row][j];
U[row][j]=t;
}
}
}
//消元
for(i=m+1;i<n;i++)
{
float t1,t2;
t1=U[i][m];
t2=U[m][m];
U[i][m]=0;
for(j=m+1;j<=n;j++)
U[i][j]=U[i][j]*t2-U[m][j]*t1;
}
}
//回代求解
for(i=n-1;i>=0;i--)
{
if(i==n-1) result[i]=U[i][i+1]/U[i][i];
else
{
sum=0;
for(j=i+1;j<n;j++)
sum=U[i][j]*result[j]+sum;
result[i]=(U[i][n]-sum)/U[i][i];
}
}
//输出根
printf("高斯选列主元消去法求得的解为:\n");
for(i=0;i<n;i++)
printf("%3.3f ",result[i]);
printf("\n");
}
‘陆’ 用c语言实现高斯消去法,解三元一次方程组。求具体程序!!
#include<iostream>
#include<cmath>
usingnamespacestd;
#defineMAX50
voidinput(doublea[MAX][MAX+1],intn)
{
cout<<"输入原方程组的增广矩阵"<<endl;
for(inti=0;i<n;i++)
for(intj=0;j<n+1;j++)
cin>>a[i][j];
}
voidoutput(doublex[],intn)
{
cout<<"Gauss消去法得到的原方程组的解为"<<endl;
for(intk=0;k<n;k++)
cout<<x[k]<<"";
}
intmain()
{
doublea[MAX][MAX+1],x[MAX],sum,max,t;
intn,i,j,k,max_i;
cout<<"输入原方程组的阶"<<endl;cin>>n;
input(a,n);
for(k=0;k<n-1;k++)//选主元素
{max=a[k][k];
max_i=k;
for(i=k+1;i<n;i++)
if(fabs(a[i][k])>fabs(max))
{
max=a[i][k];
max_i=i;
}
if(max==0)
break;
if(max_i!=k)//交换两行
for(j=k;j<n+1;j++)
{
t=a[k][j];
a[k][j]=a[max_i][j];
a[max_i][j]=t;
}
for(i=k+1;i<n;i++)
{
a[i][k]=a[i][k]/-a[k][k];
for(j=k+1;j<n+1;j++)
a[i][j]=a[i][j]+a[i][k]*a[k][j];
}//消元
}
if(max==0)cout<<"原方程组无解"<<endl;
else
{
for(k=n-1;k>=0;k--)
{
sum=0;
for(j=k+1;j<n;j++)
sum=sum+a[k][j]*x[j];
x[k]=(a[k][n]-sum)/a[k][k];
}//回代
output(x,n);
cout<<endl;
}
return0;
}