c語言解線性方程組
❶ c語言 解線性方程組
給，下面的代碼已經編譯運行確認，肯定好用了，試試吧:)
#include<conio.h>
#include<stdio.h>
#include<math.h>
#include<stdlib.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 k,i,j;
char celect;
system("cls");
printf("\n用Gauss列主元消元法解線性方程組");
printf("\n1.解方程組請按Enter.");
printf("\n2.退出程式請按Esc.");
celect=getch();
if(celect==Esc)
exit(0);
printf("\n 輸入方程組的維數：n=");
scanf("%d",&n);
printf("\n現在輸入系數矩陣A和向量b:");
for(i=1;i<=n;i++)
{
printf("\n請輸入a%d1--a%d%d系數和向量b%d: \n",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此方程組不合法!");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("\nx%d=%f",i,x[i]);
}
message();
getch();
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 繼續運算按 Enter ,退出程式按 Esc!");
switch(getch())
{
case Enter: main();
case Esc: exit(0);
default:{printf("\n不合法的輸入!");message();}
}
}
❷ 高斯先列主消元法求解線性方程組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語言求解N階線性矩陣方程Ax=b的簡單解法
基礎解系的求法
（1）對A作行初等變換，化為最簡階梯形
（2）寫出原方程組的同解方程組
（3）取定自由未知量，得基礎解系
a.每個非零行中第一個非零系數所代表的變數是主元，共R(A)個，剩餘的變數就是自由變數，共n-R(A)個；
b.在最簡階梯形矩陣中找出秩為R(A)的行列式，那麼其他各列的變數就是自由變數
3.其次線性方程組的解的判定
(1)AX=0隻有零解：R(A)=n
(2)AX=0有無窮多個非零解：R(A)=r<n
A的列向量線性相關
特別的：n階矩陣AX=0有無窮多個非零解，|A|=0
注意：若AB=0，則B的每一列都是AX=0的解
當B≠0時，意味著AX=0有非零解，進而R(B)≤n-R(A),R(A)+R(B)≤n
❹ C語言實現doolittle演算法解線性方程組
我以前剛好寫過，在我的博客里，你可以去看看
http://hi..com/ycdoit/blog/item/832586b066955d5d082302ef.html
c++解答齊次方程組的幾種方法——Cramer，Gauss列主元，Gauss全主元,Doolittle演算法/
/解線性方程組
//By JJ,2008
#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;
}
❺ 求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語言程序
void Solve ( double dCoef [] , double dY [] , unsigned int iOrder , double dErr)
{//用Jacobi迭代法解方程組, dCoef[]系數陣, Y[]向量, iOrder給出方程階數, dErr給出精度
double res [Max]; //方程解
double res2[Max]; //保存上一階方程解
if ( Max < iOrder )
{
printf ("最多支持%d階方程組.", Max);
return;
}
for ( unsigned int i = 0 ; i < iOrder ; res2 [i++] = 0.0 ); //初始解向量 (0,0...)
while ( true )
{
bool bStopIterative = true;
for (unsigned int i = 0 ; i < iOrder ; ++i)
{
double dSum2 = 0;
for (unsigned int j = 0 ; j < iOrder ; j++)
{//求第二項
if ( j == i ) continue;
dSum2 += dCoef [i * iOrder + j] * res2 [j];
}
res[i] = 1/dCoef[i * iOrder + i] * ( dY[i] - dSum2 );
if ( abs ( res2[i] - res [i] ) > dErr )
bStopIterative = false;
}
if ( bStopIterative )
break;
for (unsigned int i = 0 ; i < iOrder ; i++ )
res2[ i ] = res[ i ];
}
//輸出結果
for (unsigned int i = 0 ; i < iOrder ; i++)
printf ("x%d = %lf\n", i+1 , res[i]);
}
int main(int argc, char* argv[])
{
double a[] =
{
8,-3,2,
4,11,-1,
2,1,4
};
double b[3] =
{
20,33,12
};
Solve ( a , b , 15 , 1e-10);
getchar();
return 0;
}
❼ 用C語言編寫線性方程組數值解法程序
題目就不懂……高斯消元倒是簡單
❽ 線性方程組求解——C語言程序設計代碼，高手們，拜託啦！急急急啊！！！
這是矩陣運算的兩個文件，調用相應的函數即可：
// Matrix.h: interface for matrix calculation functions.
//
//////////////////////////////////////////////////////////////////////
#if !defined(AFX_MATRIXCALCULATE_H__CCBC1D7D_4466_4E8B_87DD_0A98B462C18D__INCLUDED_)
#define AFX_MATRIXCALCULATE_H__CCBC1D7D_4466_4E8B_87DD_0A98B462C18D__INCLUDED_
#if _MSC_VER > 1000
#pragma once
#endif // _MSC_VER > 1000
////求矩陣matrix的行列式值，n為維數
float CalculateLiner(float *matrix,int n);
////求矩陣matrix的第i行，j列的代數餘子式,n為維數
float CalculateCofacter(float *matrix,int i,int j,int n);
////matrixAT=(matrixA)T，m,n為matrixA的行、列數
void CalculateAT(float *matrixA,float *matrixAT,int m,int n);
////matrixAB=matrixA*matrixB，i,j為matrixA的行、列數，j,k為為matrixB的行、列數
void CalculateAB(float *matrixA,float *matrixB,float *matrixAB,int i,int j,int k);
////matrixATA=(matrixA)T*matrixA，m,n分別為matrixA的行、列數
void CalculateATA(float *matrixA,float *matrixATA,int m,int n);
////matrixA_為matrixA的逆，m為matrixA的行、列數
void CalculateA_(float *matrixA,float *matrixA_,int m);
///*矩陣求逆子程序(Good)*/
int Invers_matrix(float *m1,int n);
////求正定矩陣a的逆矩，n為階數
int MatrixInversion(float *a, int n);
void MatInversion(float *a,int n);
////解矩陣方程matrixA*matrixX=matrixL,m，n分別為matrixA矩陣的行，列數
void EquationResolution(float *matrixA,float *matrixL,float *matrixX,int m,int n);
#endif // !defined(AFX_MATRIXCALCULATE_H__CCBC1D7D_4466_4E8B_87DD_0A98B462C18D__INCLUDED_)
// Matrix.cpp: implementation of the matrix calculation functions.
//
//////////////////////////////////////////////////////////////////////
#include "StdAfx.h"
#include "Matrix.h"
#include "math.h"
#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif
////求矩陣matrix的行列式值，n為維數
float CalculateLiner(float *matrix,int n)
{
float liner=0;
int i=0,j=0,k=0;
int p=0,q=0;
if(n==1)
return matrix[0];
else
{
float *tr=(float *)calloc((n-1)*(n-1),sizeof(float));
for(k=0;k<n;k++)
{
p=0;
for(i=0;i<n;i++)
{
if(i!=0)
{
q=0;
for(j=0;j<n;j++)
{
if(j!=k)
{
tr[p*(n-1)+q]=matrix[i*n+j];
q++;
}
}
p++;
}
}
liner+=matrix[k]*pow(-1,k)*CalculateLiner(tr,n-1);
}
free(tr);
return liner;
}
}
////求矩陣matrix的第i行，j列的代數餘子式,n為維數
float CalculateCofacter(float *matrix,int i,int j,int n)
{
int x=0,y=0;
int p=0,q=0;
if(n==1)
return matrix[0];
else
{
float *tr=(float *)calloc((n-1)*(n-1),sizeof(float));
p=0;
for(x=0;x<n;x++)
{
if(x!=i)
{
q=0;
for(y=0;y<n;y++)
{
if(y!=j)
{
tr[p*(n-1)+q]=matrix[x*n+y];
q++;
}
}
p++;
}
}
float cc=pow(-1,i+j)*CalculateLiner(tr,n-1);
free(tr);
return cc;
}
}
////matrixAT=(matrixA)T，m,n為matrixA的行、列數
void CalculateAT(float *matrixA,float *matrixAT,int m,int n)
{
for (int i=0;i<m;i++)
{
for (int j=0;j<n;j++)
{
matrixAT[j*m+i]=matrixA[i*n+j];
}
}
}
////matrixAB=matrixA*matrixB，i,j為matrixA的行、列數，j,k為為matrixB的行、列數
void CalculateAB(float *matrixA,float *matrixB,float *matrixAB,int i,int j,int k)
{
for (int m=0;m<i;m++)
{
for (int n=0;n<k;n++)
{
matrixAB[m*k+n]=0;
for (int s=0;s<j;s++)
{
matrixAB[m*k+n]+=matrixA[m*j+s]*matrixB[s*k+n];
}
}
}
}
////matrixATA=(matrixA)T*matrixA，m,n為分別為matrixA的行、列數
void CalculateATA(float *matrixA,float *matrixATA,int m,int n)
{
float *at=(float *)calloc((m)*(n),sizeof(float));
CalculateAT(matrixA,at,m,n);
CalculateAB(at,matrixA,matrixATA,n,m,n);
free(at);
}
////matrixA_為matrixA的逆，m為matrixA的行、列數
void CalculateA_(float *matrixA,float *matrixA_,int m)
{
float liner=CalculateLiner(matrixA,m);
for(int i=0;i<m;i++)
{
for(int j=0;j<m;j++)
matrixA_[j*m+i]=CalculateCofacter(matrixA,i,j,m)/liner;
}
}
////////////////////////////////////////////////////////////////////
////求正定矩陣a的逆矩，n為階數
int MatrixInversion(float *a, int n)
{
int i, j, k, m;
float w, g, *b;
b = new float [n];
for(k = 0; k <= n - 1; k++)
{
w = a[0];
w=a[0]+1.0e-15;
/*
if(fabs(w)+1.0 == 1.0)
{
delete b;
printf("fail\n");
return(-2);
}
*/
m = n - k - 1;
for(i = 1; i <= n - 1; i++)
{
g = a[i * n];
b[i] = g / w;
if(i <= m)
b[i] = -b[i];
for(j = 1; j <= i; j++)
a[(i - 1) * n + j - 1] = a[i * n + j] + g * b[j];
}
a[n * n - 1] = 1.0 / w;
for(i = 1; i <= n - 1; i++)
a[(n - 1) * n + i - 1] = b[i];
}
for(i = 0; i <= n - 2; i++)
for(j = i + 1; j <= n - 1; j++)
a[i * n + j] = a[j * n + i];
delete b;
return(2);
}
////求正定矩陣a的逆矩，n為階數
void MatInversion(float *a,int n)
{
int i,j,k;
for(k=0;k<n;k++)
{
for(i=0;i<n;i++)
{
if(i!=k) *(a+i*n+k) = -*(a+i*n+k)/(*(a+k*n+k));
}
*(a+k*n+k)=1/(*(a+k*n+k));
for(i=0;i<n;i++)
{
if(i!=k)
{
for(j=0;j<n;j++)
{
if(j!=k) *(a+i*n+j) += *(a+k*n+j)* *(a+i*n+k);
}
}
}
for(j=0;j<n;j++)
{
if(j!=k) *(a+k*n+j)*=*(a+k*n+k);
}
}
}
/*矩陣求逆子程序(Good)*/
int Invers_matrix(float *m1,int n)
{
int *is,*js;
int i,j,k,l,u,v;
float temp,max_v;
is=(int *)malloc(n*sizeof(int));
js=(int *)malloc(n*sizeof(int));
if(is==NULL||js==NULL)
{
printf("out of memory!\n");
return(0);
}
for(k=0;k<n;k++)
{
max_v=0.0;
for(i=k;i<n;i++)
{
for(j=k;j<n;j++)
{
temp=fabs(m1[i*n+j]);
if(temp>max_v)
{
max_v=temp; is[k]=i; js[k]=j;
}
}
}
if(max_v==0.0)
{
free(is); free(js);
printf("invers is not availble!\n");
return(0);
}
if(is[k]!=k)
{
for(j=0;j<n;j++)
{
u=k*n+j; v=is[k]*n+j;
temp=m1[u]; m1[u]=m1[v]; m1[v]=temp;
}
}
if(js[k]!=k)
for(i=0;i<n;i++)
{
u=i*n+k; v=i*n+js[k];
temp=m1[u]; m1[u]=m1[v]; m1[v]=temp;
}
l=k*n+k;
m1[l]=1.0/m1[l];
for(j=0;j<n;j++)
{
if(j!=k)
{
u=k*n+j;
m1[u]*=m1[l];
}
}
for(i=0;i<n;i++)
{
if(i!=k)
{
for(j=0;j<n;j++)
{
if(j!=k)
{
u=i*n+j;
m1[u]-=m1[i*n+k]*m1[k*n+j];
}
}
}
}
for(i=0;i<n;i++)
{
if(i!=k)
{
u=i*n+k;
m1[u]*=-m1[l];
}
}
}
for(k=n-1;k>=0;k--)
{
if(js[k]!=k)
for(j=0;j<n;j++)
{
u=k*n+j; v=js[k]*n+j;
temp=m1[u]; m1[u]=m1[v]; m1[v]=temp;
}
if(is[k]!=k)
for(i=0;i<n;i++)
{
u=i*n+k; v=i*n+is[k];
temp=m1[u]; m1[u]=m1[v]; m1[v]=temp;
}
}
free(is); free(js);
return(1);
}
void EquationResolution(float *matrixA,float *matrixL,float *matrixX,int m,int n)
{
if (m<n) return;
float *at=(float *)malloc((m)*(n)*sizeof(float));
float *ata=(float *)malloc((n)*(n)*sizeof(float));
float *atl=(float *)malloc((n)*sizeof(float));
CalculateATA(matrixA,ata,m,n);
MatrixInversion(ata,n);
CalculateAT(matrixA,at,m,n);
CalculateAB(at,matrixL,atl,n,m,1);
CalculateAB(ata,atl,matrixX,n,n,1);
free(at);
free(ata);
free(atl);
}
❾ c語言解線性方程組的編程題 【做的好會多給分】
以下演算法的適用條件：A的各階主子式不為零
另外還可以採用
直接法：
消元法：Gauss-Jordan消元法，
分解法：Dolittle分解 （我用的是Courant分解法），追趕法，對稱正定矩陣的LDL『分解
----------
迭代法：
Jacobi迭代
Gauss-Seidel迭代
鬆弛迭代
-----------------
你上網可以搜索一下，或者看看數值計算方面的書
OK, 你看看這個, 另外還加了注釋 ：
Courant分解演算法：
aX = b, 作 A=LU, L是下三角矩陣， U是上三角矩陣
即L =
| L11
| L21 L22
| L31 L32 L33
| ..............
| Ln1 Ln2 ........Lnn
U =
| 1 U12 ..... U1n
| 空格 1 ..... U2n
| 空格 空格 ........
| 空格 空格 空格 空格 空格1
---------------------------------------------------
aX = b -----> LUX = b
記 UX = y,
由Ly = b得到
因為無法輸出數學符號，以下採用[i, j]Ai 表示對Ai從i到j求和
yi = (bi - [j=1, i-1]Lij yj) / Lii i = 1, 2, ..., n
由UX = y得到
xi = yi - [j=i+1, n]uij xj j = n, n-1, ..., 2, 1
你在紙上驗證一下就明白了
--------------------------------------------------------------
以下採用Courant分解 解 aX = b, 經檢查，程序運行正確
這是運行結果：
--------------------------------------------------------------
Input n value(dim of Ax=b): 3
Now input the matrix a(i, j), i, j = 0, ..., 2:
1 2 1 -2 -1 -5 0 -1 6
Now input the matrix b(i), i = 0, ..., 2:
24 -63 50
Solve...x_i =
7.000000
4.000000
9.000000
--------------------------------------------------------------
#include "stdafx.h"
#include <stdio.h>
#include <math.h>
#define MAX_N 20
int main(int argc, char* argv[])
{
int n; // 未知數個數
int i, j, k;
static double a[MAX_N][MAX_N], b[MAX_N], x[MAX_N], y[MAX_N];
static double l[MAX_N][MAX_N], u[MAX_N][MAX_N];
printf("\nInput n value(dim of Ax=b): ");
scanf("%d", &n);
if(n >MAX_N)
{
printf("The input n is larger than MAX_N, please redefine the MAX_N.\n");
return 1;
}
if(n <= 0)
{
printf("Please input a number between 1 and %d.\n", MAX_N);
return 1;
}
// {{ 程序輸入
printf("Now input the matrix a(i, j), i, j = 0, ..., %d:\n", n-1);
for (i=0; i<n; i++)
for (j=0; j<n; j++)
scanf("%lf", &a[i][j]);
printf("Now input the matrix b(i), i = 0, ..., %d:\n", n-1);
for(i=0; i<n; i++)
scanf("%lf", &b[i]);
// }} 程序輸入
for(i=0; i<n; i++)
u[i][i] = 1; //
for(k=0; k<n; k++)
{
for(i=k; i<n; i++) // 計算L的第k列元素
{
l[i][k] = a[i][k];
for(j=0; j<=k-1; j++)
l[i][k] -= (l[i][j]*u[j][k]);
}
for(j=k+1; j<n; j++) //計算U的第k行元素
{
u[k][j] = a[k][j];
for(i=0; i<=k-1; i++)
u[k][j] -= (l[k][i]*u[i][j]);
u[k][j] /= l[k][k];
}
}
for(i=0; i<n; i++) // 解Ly = b
{
y[i] = b[i];
for(j=0; j<=i-1; j++)
y[i] -= (l[i][j]*y[j]);
y[i] /= l[i][i];
}
for(i=n-1; i>=0; i--) // 解UX = Y
{
x[i]=y[i];
for(j=i+1; j<n; j++)
x[i] -= (u[i][j]*x[j]);
}
printf("Solve...x_i = \n"); // 輸出結果
for(i=0; i<n; i++)
printf("%f\n", x[i]);
return 0;
}
❿ 用高斯消去法求線性方程組，怎樣用c語言編寫
#include<conio.h>
#include<math.h>
#include <stdio.h>
#include <assert.h>
#define A(j,k) (*(A+(n+1)*j+k))
#define B(j,k) (*(B+n*j+k)) int GS(double *A, double *result , int n) //編寫高斯消去法求線性法
{ double t;
int i,j,k,tj;
for(k=0;k<n-1;k++)
{t=A(k,k);<br> tj=k;<br> for(j=k+1;j<n;j++) //A(j,k)=*(*(A+j)+k)<br> if(fabs(A(j,k))>fabs(t)){t=A(j,k);tj=j;}
if(tj!=k)
{ for(i=0;i<=n;i++)
{ t=A(tj,i);A(tj,i)=A(k,i);A(k,i)=t;}
}
for(i=k+1;i<n;i++)
{ if(A(k,k)==0) return 0;
t=A(i,k)/A(k,k);//t=i/k±�0�6
for(j=k;j<n;j++) A(i,j)-=t*A(k,j);
A(i,n)-=t*A(k,n);
}
}
for(i=n-1;i>=0;i--)
{ t=A(i,n);
for(j=n-1;j>i;j--) t-=A(i,j)*A(j,n);//X[j]=A(j,n);
if(A(i,i)==0) return 0;
A(i,n)=t/A(i,i);// X[i]=A(i,n);
}
for(i=0;i<n;i++)
result[i]=A(i,n); return 1;
}
main()
{ double *A,*result;
int i,k,n,r;
double rd;
printf("Input N: "); //輸入有幾個方程
scanf("%d",&n);
A=new double[n*(n+1)];
result=new double[n];
assert(A&&result);
printf("Input Matrix[%d][%d] In Row Order\n",n,n+1);
for(k=0;k<n*(n+1);k++) //依次輸入每個方程的左邊的系數和右邊的值
scanf("%lf",A+k);
r=GS(A, result, n);
if(r==0) printf("無解或無窮解\n");
else for(i=0;i<n;i++) printf("X%d=%.15g\n",i,result[i]);
delete[] A;
delete[] result;
getch();
}