列主元高斯消去法c語言
Ⅰ 如何用C++用列主元高斯消去法求解線性方程組的解
大二的時候自己寫得，包你滿意！
四種方法：
Cramer演算法解方程組
Gauss列主元解方程組
Gauss全主元解方程組
用Doolittle演算法解方程組
//解線性方程組
#include<iostream.h>
#include<iomanip.h>
#include<stdlib.h>
//----------------------------------------------全局變數定義區
const int Number=15; //方程最大個數
double a[Number][Number],b[Number],_a[Number][Number],_b[Number]; //系數行列式
int A_y[Number]; //a[][]中隨著橫坐標增加列坐標的排列順序,如a[0][0],a[1][2],a[2][1]...則A_y[]={0,2,1...};
int lenth,_lenth; //方程的個數
double a_sum; //計算行列式的值
char * x; //未知量a,b,c的載體
//----------------------------------------------函數聲明區
void input(); //輸入方程組
void print_menu(); //列印主菜單
int choose (); //輸入選擇
void cramer(); //Cramer演算法解方程組
void gauss_row(); //Gauss列主元解方程組
void guass_all(); //Gauss全主元解方程組
void Doolittle(); //用Doolittle演算法解方程組
int Doolittle_check(double a[][Number],double b[Number]); //判斷是否行列式>0,若是,調整為順序主子式全>0
void xiaoqu_u_l(); //將行列式Doolittle分解
void calculate_u_l(); //計算Doolittle結果
double & calculate_A(int n,int m); //計算行列式
double quanpailie_A(); //根據列坐標的排列計算的值,如A_y[]={0,2,1},得sum=a[0][ A_y[0] ] * a[1][ A_y[1] ] * a[2][ A_y[2] ]=a[0][0]*a[1][2]*a[2][1];
void exchange(int m,int i); //交換A_y[m],A_y[i]
void exchange_lie(int j); //交換a[][j]與b[];
void exchange_hang(int m,int n); //分別交換a[][]和b[]中的m與n兩行
void gauss_row_xiaoqu(); //Gauss列主元消去法
void gauss_all_xiaoqu(); //Gauss全主元消去法
void gauss_calculate(); //根據Gauss消去法結果計算未知量的值
void exchange_a_lie(int m,int n); //交換a[][]中的m和n列
void exchange_x(int m,int n); //交換x[]中的x[m]和x[n]
void recovery(); //恢復數據
//主函數
void main()
{
int flag=1;
input(); //輸入方程
while(flag)
{
print_menu(); //列印主菜單
flag=choose(); //選擇解答方式
}
}
//函數定義區
void print_menu()
{
system("cls");
cout<<"------------方程系數和常數矩陣表示如下:\n";
for(int j=0;j<lenth;j++)
cout<<"系數"<<j+1<<" ";
cout<<"\t常數";
cout<<endl;
for(int i=0;i<lenth;i++)
{
for(j=0;j<lenth;j++)
cout<<setw(8)<<setiosflags(ios::left)<<a[i][j];
cout<<"\t"<<b[i]<<endl;
}
cout<<"-----------請選擇方程解答的方案----------";
cout<<"\n 1. 克拉默(Cramer)法則";
cout<<"\n 2. Gauss列主元消去法";
cout<<"\n 3. Gauss全主元消去法";
cout<<"\n 4. Doolittle分解法";
cout<<"\n 5. 退出";
cout<<"\n 輸入你的選擇:";
}
void input()
{ int i,j;
cout<<"方程的個數:";
cin>>lenth;
if(lenth>Number)
{
cout<<"It is too big.\n";
return;
}
x=new char[lenth];
for(i=0;i<lenth;i++)
x[i]='a'+i;
//輸入方程矩陣
//提示如何輸入
cout<<"====================================================\n";
cout<<"請在每個方程里輸入"<<lenth<<"系數和一個常數:\n";
cout<<"例：\n方程:a";
for(i=1;i<lenth;i++)
{
cout<<"+"<<i+1<<x[i];
}
cout<<"=10\n";
cout<<"應輸入:";
for(i=0;i<lenth;i++)
cout<<i+1<<" ";
cout<<"10\n";
cout<<"==============================\n";
//輸入每個方程
for(i=0;i<lenth;i++)
{
cout<<"輸入方程"<<i+1<<":";
for(j=0;j<lenth;j++)
cin>>a[i][j];
cin>>b[i];
}
//備份數據
for(i=0;i<lenth;i++)
for(j=0;j<lenth;j++)
_a[i][j]=a[i][j];
for(i=0;i<lenth;i++)
_b[i]=b[i];
_lenth=lenth;
}
//輸入選擇
int choose()
{
int choice;char ch;
cin>>choice;
switch(choice)
{
case 1:cramer();break;
case 2:gauss_row();break;
case 3:guass_all();break;
case 4:Doolittle();break;
case 5:return 0;
default:cout<<"輸入錯誤,請重新輸入:";
choose();
break;
}
cout<<"\n是否換種方法求解(Y/N):";
cin>>ch;
if(ch=='n'||ch=='N') return 0;
recovery();
cout<<"\n\n\n";
return 1;
}
//用克拉默法則求解方程.
void cramer()
{
int i,j;double sum,sum_x;char ch;
//令第i行的列坐標為i
cout<<"用克拉默(Cramer)法則結果如下:\n";
for(i=0;i<lenth;i++)
A_y[i]=i;
sum=calculate_A(lenth,0);
if(sum!=0)
{
cout<<"系數行列式不為零,方程有唯一的解:";
for(i=0;i<lenth;i++)
{ ch='a'+i;
a_sum=0;
for(j=0;j<lenth;j++)
A_y[j]=j;
exchange_lie(i);
sum_x=calculate_A(lenth,0);
cout<<endl<<ch<<"="<<sum_x/sum;
exchange_lie(i);
}
}
else
{
cout<<"系數行列式等於零,方程沒有唯一的解.";
}
cout<<"\n";
}
double & calculate_A(int n,int m) //計算行列式
{ int i;
if(n==1) {
a_sum+= quanpailie_A();
}
else{for(i=0;i<n;i++)
{ exchange(m,m+i);
calculate_A(n-1,m+1);
exchange(m,m+i);
}
}
return a_sum;
}
double quanpailie_A() //計算行列式中一種全排列的值
{
int i,j,l;
double sum=0,p;
for(i=0,l=0;i<lenth;i++)
for(j=0;A_y[j]!=i&&j<lenth;j++)
if(A_y[j]>i) l++;
for(p=1,i=0;i<lenth;i++)
p*=a[i][A_y[i]];
sum+=p*((l%2==0)?(1):(-1));
return sum;
}
//高斯列主元排列求解方程
void gauss_row()
{
int i,j;
gauss_row_xiaoqu(); //用高斯列主元消區法將系數矩陣變成一個上三角矩陣
for(i=0;i<lenth;i++)
{
for(j=0;j<lenth;j++)
cout<<setw(10)<<setprecision(5)<<a[i][j];
cout<<setw(10)<<b[i]<<endl;
}
if(a[lenth-1][lenth-1]!=0)
{
cout<<"系數行列式不為零,方程有唯一的解：\n";
gauss_calculate();
for(i=0;i<lenth;i++) //輸出結果
{
cout<<x[i]<<"="<<b[i]<<"\n";
}
}
else
cout<<"系數行列式等於零,方程沒有唯一的解.\n";
}
void gauss_row_xiaoqu() //高斯列主元消去法
{
int i,j,k,maxi;double lik;
cout<<"用Gauss列主元消去法結果如下:\n";
for(k=0;k<lenth-1;k++)
{
j=k;
for(maxi=i=k;i<lenth;i++)
if(a[i][j]>a[maxi][j]) maxi=i;
if(maxi!=k)
exchange_hang(k,maxi);//
for(i=k+1;i<lenth;i++)
{
lik=a[i][k]/a[k][k];
for(j=k;j<lenth;j++)
a[i][j]=a[i][j]-a[k][j]*lik;
b[i]=b[i]-b[k]*lik;
}
}
}
//高斯全主元排列求解方程
void guass_all()
{
int i,j;
gauss_all_xiaoqu();
for(i=0;i<lenth;i++)
{
for(j=0;j<lenth;j++)
cout<<setw(10)<<setprecision(5)<<a[i][j];
cout<<setw(10)<<b[i]<<endl;
}
if(a[lenth-1][lenth-1]!=0)
{
cout<<"系數行列式不為零,方程有唯一的解：\n";
gauss_calculate();
for(i=0;i<lenth;i++) //輸出結果
{
for(j=0;x[j]!='a'+i&&j<lenth;j++);
cout<<x[j]<<"="<<b[j]<<endl;
}
}
else
cout<<"系數行列式等於零,方程沒有唯一的解.\n";
}
void gauss_all_xiaoqu() //Gauss全主元消去法
{
int i,j,k,maxi,maxj;double lik;
cout<<"用Gauss全主元消去法結果如下:\n";
for(k=0;k<lenth-1;k++)
{
for(maxi=maxj=i=k;i<lenth;i++)
{
for(j=k;j<lenth;j++)
if(a[i][j]>a[maxi][ maxj])
{ maxi=i;
maxj=j;
}
}
if(maxi!=k)
exchange_hang(k,maxi);
if(maxj!=k)
{
exchange_a_lie(maxj,k); //交換兩列
exchange_x(maxj,k);
}
for(i=k+1;i<lenth;i++)
{
lik=a[i][k]/a[k][k];
for(j=k;j<lenth;j++)
a[i][j]=a[i][j]-a[k][j]*lik;
b[i]=b[i]-b[k]*lik;
}
}
}
void gauss_calculate() //高斯消去法以後計算未知量的結果
{
int i,j;double sum_ax;
b[lenth-1]=b[lenth-1]/a[lenth-1][lenth-1];
for(i=lenth-2;i>=0;i--)
{
for(j=i+1,sum_ax=0;j<lenth;j++)
sum_ax+=a[i][j]*b[j];
b[i]=(b[i]-sum_ax)/a[i][i];
}
}
void Doolittle() //Doolittle消去法計算方程組
{
double temp_a[Number][Number],temp_b[Number];int i,j,flag;
for(i=0;i<lenth;i++)
for(j=0;j<lenth;j++)
temp_a[i][j]=a[i][j];
flag=Doolittle_check(temp_a,temp_b);
if(flag==0) cout<<"\n行列式為零.無法用Doolittle求解.";
xiaoqu_u_l();
calculate_u_l();
cout<<"用Doolittle方法求得結果如下:\n";
for(i=0;i<lenth;i++) //輸出結果
{
for(j=0;x[j]!='a'+i&&j<lenth;j++);
cout<<x[j]<<"="<<b[j]<<endl;
}
}
void calculate_u_l() //計算Doolittle結果
{ int i,j;double sum_ax=0;
for(i=0;i<lenth;i++)
{
for(j=0,sum_ax=0;j<i;j++)
sum_ax+=a[i][j]*b[j];
b[i]=b[i]-sum_ax;
}
for(i=lenth-1;i>=0;i--)
{
for(j=i+1,sum_ax=0;j<lenth;j++)
sum_ax+=a[i][j]*b[j];
b[i]=(b[i]-sum_ax)/a[i][i];
}
}
void xiaoqu_u_l() //將行列式按Doolittle分解
{ int i,j,n,k;double temp;
for(i=1,j=0;i<lenth;i++)
a[i][j]=a[i][j]/a[0][0];
for(n=1;n<lenth;n++)
{ //求第n+1層的上三角矩陣部分即U
for(j=n;j<lenth;j++)
{ for(k=0,temp=0;k<n;k++)
temp+=a[n][k]*a[k][j];
a[n][j]-=temp;
}
for(i=n+1;i<lenth;i++) //求第n+1層的下三角矩陣部分即L
{ for(k=0,temp=0;k<n;k++)
temp+=a[i][k]*a[k][n];
a[i][n]=(a[i][n]-temp)/a[n][n];
}
}
}
int Doolittle_check(double temp_a[][Number],double temp_b[Number]) //若行列式不為零,將系數矩陣調整為順序主子式大於零
{
int i,j,k,maxi;double lik,temp;
for(k=0;k<lenth-1;k++)
{
j=k;
for(maxi=i=k;i<lenth;i++)
if(temp_a[i][j]>temp_a[maxi][j]) maxi=i;
if(maxi!=k)
{ exchange_hang(k,maxi);
for(j=0;j<lenth;j++)
{ temp=temp_a[k][j];
temp_a[k][j]=temp_a[maxi][j];
temp_a[maxi][j]=temp;
}
}
for(i=k+1;i<lenth;i++)
{
lik=temp_a[i][k]/temp_a[k][k];
for(j=k;j<lenth;j++)
temp_a[i][j]=temp_a[i][j]-temp_a[k][j]*lik;
temp_b[i]=temp_b[i]-temp_b[k]*lik;
}
}
if(temp_a[lenth-1][lenth-1]==0) return 0;
return 1;
}
void exchange_hang(int m,int n) //交換a[][]中和b[]兩行
{
int j; double temp;
for(j=0;j<lenth;j++)
{ temp=a[m][j];
a[m][j]=a[n][j];
a[n][j]=temp;
}
temp=b[m];
b[m]=b[n];
b[n]=temp;
}
void exchange(int m,int i) //交換A_y[m],A_y[i]
{ int temp;
temp=A_y[m];
A_y[m]=A_y[i];
A_y[i]=temp;
}
void exchange_lie(int j) //交換未知量b[]和第i列
{ double temp;int i;
for(i=0;i<lenth;i++)
{ temp=a[i][j];
a[i][j]=b[i];
b[i]=temp;
}
}
void exchange_a_lie(int m,int n) //交換a[]中的兩列
{ double temp;int i;
for(i=0;i<lenth;i++)
{ temp=a[i][m];
a[i][m]=a[i][n];
a[i][n]=temp;
}
}
void exchange_x(int m,int n) //交換未知量x[m]與x[n]
{ char temp;
temp=x[m];
x[m]=x[n];
x[n]=temp;
}
void recovery() //用其中一種方法求解後恢復數據以便用其他方法求解
{
for(int i=0;i<lenth;i++)
for(int j=0;j<lenth;j++)
a[i][j]=_a[i][j];
for(i=0;i<lenth;i++)
b[i]=_b[i];
for(i=0;i<lenth;i++)
x[i]='a'+i;
a_sum=0;
lenth=_lenth;
}
#include <iostream>
#include <iomanip.h>
using namespace std;
#define N 20
double a[N][N];
double x[N+1];
double b[N+1];
int n;//n方程個數，n未知數個數
int set( )
{
cout<<"請輸入方程的個數和未知數個數： "<<endl;
cin>>n;
int i,j;
cout<<"請輸入方程組（逐個輸入方程 i）"<<endl;
for(i = 1;i <= n;i++)
{
for(j = 1;j<=n;j++)
{
cin>>a[i][j];//系數
}
cin>>b[i];//結果
}
return 0;
}
int find(int k)//尋找第k列主元
{
int i,tag = k;
double maxv=0;
for(i = k;i <= n;i++)
{
if(a[i][k] > maxv)
{
maxv = a[i][k];
tag = i;
}
}
return tag;
}
void exchange(int i1,int i2)//將第 i1 i2行互換
{
int j;
for(j = 1;j <= n;j++)
{
swap(a[i1][j],a[i2][j]);
}
swap(b[i1],b[i2]);
}
void div(int k)//讓第k個方程的首項系數為1
{
double temp = a[k][k];
for(int j = k; j <= n;j++)
{
a[k][j]/=temp;
}
b[k]/=temp;
}
void disME(int k)
{
int i,j;
for(i =1 ;i<=n;i++)
{
for(j = i;j<= n;j++)
{
if(a[i][j])
{
if(a[i][j]==1)
{ if(j==n)
cout<<"x"<<j;
else
cout<<"x"<<j<<" + ";
}
else
{
if(j==n)
cout<<a[i][j]<<"x"<<j;
else
cout<<a[i][j]<<"x"<<j<<" + ";
}
}
}
cout<<" = "<<b[i]<<endl;
}
system("pause");
}
void eliminate(int k)//消元
{
int i,j;
double t;
for(i = k+1;i<= n;i++)
{
t = a[i][k]/a[k][k];
for(j = k;j <= n;j++)
{
a[i][j]-=a[k][j] * t;
}
b[i] -= b[k] * t;
}
}
void Gauss()//高斯消元法
{
int i,j,k;
for(k = 1;k < n;k++)//共進行n - 1次消元
{
int l = find(k);//尋找第k次消元法的主元方程
if(l!=k)
{
exchange(l,k);
}
//消元
div(k);
eliminate(k);
cout<<"第"<<k<<"次消元結果："<<endl;
disME(k);
}
div(k);
x[k] = b[k];
//disM();
for(i = n - 1;i>=1;i--)
{
for(j = i+1;j <=n ;j++)
{
b[i] -= a[i][j] * b [j];
}
x[i] = b[i];
}
}
void dis()
{
int i;
cout<<"解方程得："<<endl;
for(i = 1;i<=n;i++)
{
cout<<"x"<<i<<" = ";
printf("%.5f\n",x[i]);
}
}
int main()
{
set();
Gauss();
dis();
system("pause");
return 0;
}
Ⅲ 高斯列主元消去法,求解齊次線性方程組的C語言實現
C/C++ code #include<stdio.h> #include <math.h> #define N 20 int main() { int n,i,j,k; int mi,tmp,mx; float a[N][N],b[N],x[N]; printf("\nInput n:"); scanf("%d",&n); if(n>N) { printf("The input n should in(0,N)!\n"); getch(); return 1; } if(n<=0) { printf("The input n should in(0,N)!\n"); getch(); return 1; } printf("Now input a(i,j),i,j=0...%d:\n",n-1); for(i=0;i<n;i++) { for(j=0;j<n;j++) scanf("%f",&a[i][j]);} printf("Now input b(i),i,j=0...%d:\n",n-1); for(i=0;i<n;i++) scanf("%f",&b[i]); for(i=0;i<n-2;i++) { for(j=i+1,mi=i,mx=fabs(a[i][j]);j<n-1;j++) if(fabs(a[j][i])>mx) { mi=j; mx=fabs(a[j][i]); } if(i<mi) { tmp=b[i];b[i]=b[mi];b[mi]=tmp; for(j=i;j<n;j++) { tmp=a[i][j]; a[i][j]=a[mi][j]; a[mi][j]=tmp; } } for(j=i+1;j<n;j++) { tmp=-a[j][i]/a[i][i]; b[j]+=b[i]*tmp; for(k=i;k<n;k++) a[j][k]+=a[i][k]*tmp; } } x[n-1]=b[n-1]/a[n-1][n-1]; for(i=n-2;i>=0;i--) { x[i]=b[i]; for(j=i+1;j<n;j++) x[i]-=a[i][j]*x[j]; x[i]/=a[i][i]; } for(i=0;i<n;i++) printf("Answer:\n x[%d]=%f\n",i,x[i]); getch(); return 0; } #include<math.h> #include<stdio.h> #define NUMBER 20 #define Esc 0x1b #define Enter 0x0d float A[NUMBER][NUMBER+1] ,ark; int flag,n; exchange(int r,int k); float max(int k); message(); main() { float x[NUMBER]; int r,k,i,j; char celect; clrscr(); printf("\n\nUse Gauss."); printf("\n\n1.Jie please press Enter."); printf("\n\n2.Exit press Esc."); celect=getch(); if(celect==Esc) exit(0); printf("\n\n input n="); scanf("%d",&n); printf(" \n\nInput matrix A and B:"); for(i=1;i<=n;i++) { printf("\n\nInput a%d1--a%d%d and 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\nIt's wrong!");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(); } 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; } message() { printf("\n\n Go on Enter ,Exit press Esc!"); switch(getch()) { case Enter: main(); case Esc: exit(0); default:{printf("\n\nInput error!");message();} } }
Ⅳ 如何用編程實現高斯列主元消去法
用C語言編制以下程序：
#include <stdio.h>
#include <math.h>
void main()
{
int i,j,k,ik,n;
float max,t,a[10][10],x[10],sum;
printf("請輸入方程組的階數:");
scanf("%d",&n);
printf("請輸入增廣矩陣\n");
for (i=1;i<=n;i++)
for (j=1;j<=n+1;j++)
scanf("%f",&a[i][j]);
for (k=1;k<=n-1;k++)
{
max=0.0;
for (i=k;i<=n;i++)
if (max<fabs(a[i][k]))
{
max=fabs(a[i][k]);
ik=i;
}
if (max==0)
{
printf("A is singular");
break;
}
if (ik!=k)
for (j=k;j<=n+1;j++)
{
t=a[k][j];
a[k][j]=a[ik][j];
a[ik][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 (a[n][n]==0)
printf("A is singular");
else
x[n]=a[n][n+1]/a[n][n];
for (k=n-1;k>=1;k--)
{
sum=0.0;
for (j=k+1;j<=n;j++)
sum=sum+a[k][j]*x[j];
x[k]=(a[k][n+1]-sum)/a[k][k];
}
for (i=1;i<=n;i++)
printf("x[%d]=%f\n",i,x[i]);
}
Ⅳ 用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語言編譯 採用高斯先列主元消元法求解線性方程組AX=b
545
Ⅶ 求C語言課程設計：用高斯列主元消元法解線性方程組
這里向你推薦一下克魯特演算法（其實就是對高斯列主元消元法進行優化，使之更適合於計算機編程），首先將矩陣A進行LU分解（將系數矩陣分解成一個上三角矩陣和一個下三角矩陣），分解的過程中用到了隱式的主元尋找法，同時利用克魯特演算法可以將兩個n*n矩陣壓縮到一個n*n矩陣中，大大節省了存儲空間提高了計算速度。
方程可化為L*U*x=B，令U*x=y --->L*y=B
然後利用回代先求y，再利用y求x
因為該方法在求解過程中不涉及增廣矩陣所以矩陣B幾乎不參與什麼運算，所以它的計算速度應該能夠達到高斯列主元消元法的三倍，但原理與其基本一致。
而且我在程序中使用了動態數組方便你今後進行擴展。
以下程序按照《矩陣論第二版》和《C語言數值計演算法方法大全》編寫，LU分解部分程序主要參考了《C語言數值計演算法方法大全》第二章的程序
如果你需要詳細的理論講解我可以將這兩本書和源程序發給你，上面的論述相當詳細足夠你答辯用的了，我的郵箱[email protected]
計算結果：
A矩陣:
2 2 5
3 4 7
1 3 3
B矩陣:
5
6
5
解矩陣:
x 1=-7
x 2=0.333333
x 3=3.66667
Press any key to continue
#include <cmath>
#include <iostream>
#include <iomanip>
#include <cstdlib>
#include <functional>
#include <vector>
#include <algorithm>
using namespace std;
#define TINY 1.0e-20 //A small number.
#define N 3
void ludcmp(vector<vector<float> > &a, int n, vector<int> &indx, float &d);//對矩陣進行LU分解
void lubksb(vector<vector<float> > &a, int n, vector<int> &indx, vector<float> &b);//對矩陣進行前向和後向回代
void root(vector<vector<float> > &x,vector<float> &col);//解方程結果保存在y中
void iniv(vector<vector<float> > &x,vector<float> line,int n);//對二維動態數組進行初始化
void main()
{
int i,j,n=N;//輸入矩陣的維數
float A[N][N]={{2,2,5},{3,4,7},{1,3,3}};//左邊A矩陣
float B[N]={5,6,5};//右邊B矩陣
vector<vector<float> > x;//建立動態二維數組存放A，保證你的程序進行擴展時只改A，B，N
vector<float> line;
vector<float> y(n);//建立動態數組存放B
iniv(x,line,n);
y.clear();
for(i=0;i<n;i++)//將A賦給x，B賦給y
{
y.push_back(B[i]);
for(j=0;j<n;j++)
{
x[i].push_back(A[i][j]);
}
}
cout<<"A矩陣:"<<endl;
for(i=0;i<n;i++)
{
for(j=0;j<n;j++)
{
cout<<setw(2)<<setiosflags(ios::left)<<setw(2)<<x[i][j]<<" ";
}
cout<<endl;
}
cout<<"B矩陣:"<<endl;
for(i=0;i<n;i++)
{
cout<<setw(2)<<setiosflags(ios::left)<<setw(2)<<y[i]<<endl;
}
root(x,y);//求根
cout<<"解矩陣:"<<endl;
for(i=0;i<n;i++)
{
cout<<setw(2)<<setiosflags(ios::left)<<"x"<<i+1<<"="<<y[i]<<endl;
}
cout<<endl;
}
void root(vector<vector<float> > &x,vector<float> &col)
{
int n=x.size(),i=0,j=0;
vector<int> index(n);//用於記錄尋找主元素過程中對矩陣的初等變換
index.clear();
float m=1.0;//記錄變換方式，此程序中無用
ludcmp(x,n,index,m);//進行LU分解
lubksb(x,n,index,col);//根據分解結果進行回帶
}
//以下程序按照《矩陣論第二版》和《C語言數值計演算法方法大全》編寫，LU分解部分程序主要參考了《C語言數值計演算法方法大全》第二章的程序
//如果你需要詳細的理論講解我可以將這兩本書和源程序發給你，我的郵箱[email protected]
void ludcmp(vector<vector<float> > &a, int n, vector<int> &indx, float &d)
{
int i,imax,j,k;
float big=0,m=0,sum=0,temp=0;
vector<float> vv(n);
vv.clear();
d=1.0;
for (i=0;i<n;i++)
{
big=0.0;
for (j=0;j<n;j++)
if ((temp=fabs(a[i][j])) > big)
big=temp;
vv[i]=1.0/big;
}
for (j=0;j<n;j++)
{
for (i=0;i<j;i++)
{
sum=a[i][j];
for (k=0;k<i;k++)
sum -= a[i][k]*a[k][j];
a[i][j]=sum;
}
big=0.0;
for (i=j;i<n;i++)
{
sum=a[i][j];
for (k=0;k<j;k++)
sum -= a[i][k]*a[k][j];
a[i][j]=sum;
if ( (m=vv[i]*fabs(sum)) >= big)
{
big=m;
imax=i;
}
}
if (j != imax)
{
for (k=0;k<n;k++)
{
m=a[imax][k];
a[imax][k]=a[j][k];
a[j][k]=m;
}
d = -(d);
vv[imax]=vv[j];
}
indx[j]=imax;
if (a[j][j] == 0.0)
a[j][j]=TINY;
if (j != n)
{
m=1.0/(a[j][j]);
for (i=j+1;i<n;i++)
a[i][j] *= m;
}
}
}
void lubksb(vector<vector<float> > &a, int n, vector<int> &indx, vector<float> &b)
{
int i,ii=0,ip,j;
float sum;
for(i=0;i<n;i++)//按LU分解時尋找主元所進行的初等變換進行反邊變換。
{
ip=indx[i];
sum=b[ip];
b[ip]=b[i];
b[i]=sum;
}
sum=0;
for (i=1;i<n;i++)
{
sum=0;
for(j=0;j<i;j++)
{
sum+=a[i][j]*b[j];
}
b[i]=b[i]-sum;
}
b[n-1]=b[n-1]/a[n-1][n-1];
for (i=n-2;i>=0;i--)
{
sum=0;
for(j=i+1;j<n;j++)
{
sum+=a[i][j]*b[j];
}
b[i]=(b[i]-sum)/a[i][i];
}
}
void iniv(vector<vector<float> > &x,vector<float> line,int n)
{
int i,j;
for(i=0;i<n;i++)
{
x.push_back(line);
for(j=0;j<n;j++)
{
x[i].clear();
}
}
}
Ⅷ 用c語言編寫高斯-約當全主元消去法
#include <iostream>
#include <iomanip.h>
using namespace std;
#define N 20
double a[N][N];
double x[N+1];
double b[N+1];
int n;//n方程個數，n未知數個數
int set( )
{
cout<<"請輸入方程的個數和未知數個數： "<<endl;
cin>>n;
int i,j;
cout<<"請輸入方程組（逐個輸入方程 i）"<<endl;
for(i = 1;i <= n;i++)
{
for(j = 1;j<=n;j++)
{
cin>>a[i][j];//系數
}
cin>>b[i];//結果
}
return 0;
}
int find(int k)//尋找第k列主元
{
int i,tag = k;
double maxv=0;
for(i = k;i <= n;i++)
{
if(a[i][k] > maxv)
{
maxv = a[i][k];
tag = i;
}
}
return tag;
}
void exchange(int i1,int i2)//將第 i1 i2行互換
{
int j;
for(j = 1;j <= n;j++)
{
swap(a[i1][j],a[i2][j]);
}
swap(b[i1],b[i2]);
}
void div(int k)//讓第k個方程的首項系數為1
{
double temp = a[k][k];
for(int j = k; j <= n;j++)
{
a[k][j]/=temp;
}
b[k]/=temp;
}
void disME(int k)
{
int i,j;
for(i =1 ;i<=n;i++)
{
for(j = i;j<= n;j++)
{
if(a[i][j])
{
if(a[i][j]==1)
{ if(j==n)
cout<<"x"<<j;
else
cout<<"x"<<j<<" + ";
}
else
{
if(j==n)
cout<<a[i][j]<<"x"<<j;
else
cout<<a[i][j]<<"x"<<j<<" + ";
}
}
}
cout<<" = "<<b[i]<<endl;
}
system("pause");
}
void eliminate(int k)//消元
{
int i,j;
double t;
for(i = k+1;i<= n;i++)
{
t = a[i][k]/a[k][k];
for(j = k;j <= n;j++)
{
a[i][j]-=a[k][j] * t;
}
b[i] -= b[k] * t;
}
}
void Gauss()//高斯消元法
{
int i,j,k;
for(k = 1;k < n;k++)//共進行n - 1次消元
{
int l = find(k);//尋找第k次消元法的主元方程
if(l!=k)
{
exchange(l,k);
}
//消元
div(k);
eliminate(k);
cout<<"第"<<k<<"次消元結果："<<endl;
disME(k);
}
div(k);
x[k] = b[k];
//disM();
for(i = n - 1;i>=1;i--)
{
for(j = i+1;j <=n ;j++)
{
b[i] -= a[i][j] * b [j];
}
x[i] = b[i];
}
}
void dis()
{
int i;
cout<<"解方程得："<<endl;
for(i = 1;i<=n;i++)
{
cout<<"x"<<i<<" = ";
printf("%.5f\n",x[i]);
}
}
int main()
{
set();
Gauss();
dis();
system("pause");
return 0;
}
———————————————————————————————
希望答案能幫到你，要是你有疑問可以追問
當然，採納也是您的美德，謝謝
Ⅸ 用C語言進行列主元素高斯消元法，遇到問題
for(i=k+1;i<N;i++)
{
m[i-1][k]=a[i][k]/a[k][k];
//這里m的坐標應該是[i-1][k]，如果是[i][k]會造成越界
for(j=0;j<N;j++)
{
temp=a[i][j];
a[i][j]=temp-m[i-1][k]*a[k][j];
//這里也一樣
}
}
m是
2X2
的數組，而a是
3X3
的數據，即a[1][0]與a[0][0]的比值應存在m[0][0]中！
希望可以幫到你！