① 用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;
}
② 用高斯消元法解三元一次方程組，C語言
參閱我的文章：http://wenku..com/view/d4ea2273650e52ea5418981d.html
#include "stdafx.h" //VS 預編譯頭文件，其他系統請刪除
#include<stdio.h>
#include<stdlib.h>
#include<memory.h>
#include<math.h>
#include<time.h>
//VS 2013 否決了 scanf 等函數，為了使用，加上下句。
//其他系統請刪除
#pragma warning(disable:4996)
int GaussJordanElimination(int n, const double *pCoef, double *pOut);
//VS 主函數簽名格式。其他系統請改變簽名，如：
//int main()
int _tmain(int argc, _TCHAR* argv[])
{
double cf[3][4] = { {-0.02, 2.0, 2.0, 0.4}, {1.0, 0.78125, 0.0, 1.3816}, {3.996, 5.526, 4.0, 7.4178} };
double rs[3];
int i;
i = GaussJordanElimination(3, (double*)cf, rs);
printf("x1 = %lf, x2 = %lf, x3 = %lf\n", rs[0], rs[1], rs[2]);
system("pause"); //避免窗口一閃而退
return 0;
}
//絕對值函數
__inline double _abs(double v)
{
return v < 0 ? -v : v;
}
//線性方程組列主元高斯消元法
//n 方程元數；pCoef 系數，必須以行主序方式存放的二維數組；
//pOut 長度為 n 的一維數組(調用者負責維護)，用於輸出數據
//返回值：0 成功，-1 無解，1 申請內存失敗, 2 不定解。
int GaussJordanElimination(int n, const double *pCoef, double *pOut)
{
double *pcf;
int rows = n, columns = n + 1;
//pcf = new double[rows * columns];
pcf = (double*)malloc(rows * columns * sizeof(double));
if (pcf == 0) return 1; //巧婦難為無米之炊，內存都申請不到，還能幹嘛！
memcpy(pcf, pCoef, (rows * columns) * sizeof(double)); //據說這個運行效率很高
int r, c, i; //循環變數
int a, b;
double x, y;
//開始消元，將 pcf 方陣區處理成到直角三角形（直角在右上角）矩陣
for (r = 0; r < rows - 1; r++)
{
//選取主元
a = r; x = _abs(pcf[r * columns + r]);
for (i = r + 1; i < rows; i++)
{ //查找主元在哪行
if (x < _abs(pcf[i * columns + r])) a = i;
}
if (a > r)
{ //主元不是當前行（r）,比較麻煩，需要將第 a 行與第 r 行兌換
//第 r 列前面的就不要對換了，因為這些項已經被消元，變成 0 了
for (c = r; c < columns; c++)
{
x = pcf[r * columns + c];
pcf[r * columns + c] = pcf[a * columns + c];
pcf[a * columns + c] = x;
}
}
//開始消元
a = r * columns; //記住將主元的行地址偏移量，以提高程序運行效率
x = -pcf[a + r]; //要多次使用，記下她，以提高程序運行效率
if (x == 0) //主元居然為 0，純粹是想坑爹，豈能上當！
continue; //繼續後面的消元，以便最終判斷是無解還是任意解
for (i = r + 1; i < rows; i++)
{ //正在消元
b = i * columns;//記住將要消元的行地址偏移量，以提高程序運行效率
y = pcf[b + r]; //要多次使用，記下她，以提高程序運行效率
if (y != 0)
{ //y == 0，本行不需要消元
y /= x; //要多次使用，記下她，以提高程序運行效率
pcf[b + r] = 0; //肯定為 0，不用計算。
for (c = r + 1; c < columns; c++)
pcf[b + c] += pcf[a + c] * y;
}
}
}//至此，pcf 方陣區已經處理成到直角三角形（直角在右上角）矩陣
//回代，將 pcf 方陣區處理成主對角線為 1，其他為 0 的矩陣
int columns_1 = c = columns - 1; //多次用到，提高效率
for (r = rows - 1; r >= 1; r--)
{
b = r * columns;
if (pcf[b + r] == 0)
{ //經過前面的消元，除主元外，其他元應該都為 0
if (pcf[b + columns_1] == 0)
{ //常數項為 0，方程有不定解
free(pcf);
return 2;
}
else
{ //常數項為 0，方程有無解
free(pcf); //釋放內存
return -1;
}
}
pcf[b + columns_1] /= pcf[b + r];
pcf[b + r] = 1; //肯定為 1，不用計算。
y = -pcf[b + columns_1];
//回代
for (i = r - 1; i >= 0; i--)
{
pcf[i * columns + columns_1] += pcf[i * columns + r] * y;
pcf[i * columns + r] = 0; //已經回代，此項已消，置為 0。
}
}
//處理第一行數據
pcf[columns_1] /= pcf[0];
pcf[0] = 1;
//至此，回代過程結束，pcf 矩陣的最後一列就是結果
//返回結果
for (r = 0; r < rows; r++)
pOut[r] = pcf[r * columns + columns_1];
free(pcf);
return 0;
}
③ 【編程求助】用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語言
其中用到了高斯先列主消元法 #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語言編寫一個程序用高斯消元法求解線性方程組的解程序流程圖
#include <iostream>
#include <cmath>
using namespace std;
int main()
{
int n,i,j,k;
double a[100][100],b[100],o;
cout<<"輸入未知數個數:"<<endl;
cin>>n;
cout<<"輸入"<<n<<"個方程組(未知數系數及等號右邊常數):"<<endl;
for (i=1;i<=n;i++)
for (j=1;j<=n+1;j++)
cin>>a[i][j];
for (i=1;i<=n;i++)
for (j=i+1;j<=n;j++)
if (fabs(a[j][i])>1e-7)
{
o=a[i][i]/a[j][i];
for (k=i;k<=n+1;k++)
a[j][k]=a[j][k]*o-a[i][k];
}
for (i=n;i>0;i--)
{
b[i]=a[i][n+1]/a[i][i];
for (j=i-1;j>0;j--)
a[j][n+1]=a[j][n+1]-b[i]*a[j][i];
}
cout<<"解得:"<<endl;
for (i=1;i<=n;i++) cout<<b[i]<<endl;
system("pause");
return 0;
}
⑥ 用C語言寫一個高斯消元法解方程組的程序
我們以方程組 2x1 + 6x2 - x3 = -12
5x1 - x2 +2x3 = 29
-3x1 - 4x2 + x3 = 5
為例 來說明樓主自己把方程組化為矩陣形式。以下為源代碼 。
#include <stdio.h>
#include <stdlib.h>
#include <malloc.h>
#include <math.h>
int GS(int,double**,double *,double);
double **TwoArrayAlloc(int,int);
void TwoArrayFree(double **);
int main(void)
{
int i,n;
double ep,**a,*b;
n = 3;
ep = 1e-4;
a = TwoArrayAlloc(n,n);
b = (double *)calloc(n,sizeof(double));
if(b == NULL)
{
printf("memory get error\n");
exit(1);
}
a[0][0]= 2; a[0][1]= 6; a[0][2]= -1;
a[1][0]= 5; a[1][1]=-1; a[1][2]= 2;
a[2][0]=-3; a[2][1]=-4; a[2][2]= 1;
b[0] = -12; b[1] = 29; b[2] = 5;
if(!GS(n,a,b,ep))
{
printf("can't solve it with GS elimination\n");
exit(0);
}
printf("The solution of equations is as follows：\n");
for(i=0;i<3;i++)
{
printf("x%d = %.2f\n",i,b[i]);
}
TwoArrayFree(a);
free(b);
return 0;
}
int GS(n,a,b,ep)
int n;
double **a;
double *b;
double ep;
{
int i,j,k,l;
double t;
for(k=1;k<=n;k++)
{
for(l=k;l<=n;l++)
if(fabs(a[l-1][k-1])>ep)
break;
else if(l==n)
return(0);
if(l!=k)
{
for(j=k;j<=n;j++)
{
t = a[k-1][j-1];
a[k-1][j-1] =a[l-1][j-1];
a[l-1][j-1] =t;
}
t=b[k-1];
b[k-1]=b[l-1];
b[l-1]=t;
}
t=1/a[k-1][k-1];
for(j=k+1;j<=n;j++)
a[k-1][j-1]=t*a[k-1][j-1];
b[k-1]*=t;
for(i=k+1;i<=n;i++)
{
for(j=k+1;j<=n;j++)
a[i-1][j-1]-=a[i-1][k-1]*a[k-1][j-1];
b[i-1]-=a[i-1][k-1]*b[k-1];
}
}
for(i=n-1;i>=1;i--)
for(j=i+1;j<=n;j++)
b[i-1]-=a[i-1][j-1]*b[j-1];
return(1);
}
double **TwoArrayAlloc(int r,int c)
{
double *x,**y;
int n;
x=(double *)calloc(r*c,sizeof(double));
y=(double **)calloc(r,sizeof(double*));
for(n=0;n<=r-1;++n)
{
y[n]=&x[c*n];
}
return y ;
}
void TwoArrayFree(double **x)
{
free(x[0]);
free(x);
}
⑦ 求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();
}
}
}
⑧ 用無回代的高斯消元法和簡單迭代法 解四元一次方程組
N元通用程序，用四元調用.四元方程組:
a[4][5]={4.,-1.,0.,2.,-1.0,-1,4,-1,0,-7.0,0,-1,4,-1,9.0,2,0,-1,4.0,0.0};
高斯消元法（通用程序）：
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
float *ColPivot(float *c,int n);
void main()
{
int i,j;
float *x;
int n=4;
float a[4][5]={4.,-1.,0.,2.,-1.0,-1,4,-1,0,-7.0,0,-1,4,-1,9.0,2,0,-1,4.0,0.0};
x=ColPivot(a[0],n);
for (i=0;i<n;i++) printf("x[%d]=%f\n",i,x[i]);
getch();
}
float *ColPivot(float *c,int n)
{
int i,j,t,k;
float *x,p;
x=(float *)malloc(n*sizeof(float));
for(i=0;i<=n-2;i++)
{
k=i;
for(j=i+1;j<=n-1;j++)
if(fabs(*(c+j*(n+1)+i))>(fabs(*(c+k*(n+1)+i)))) k=j;
if(k!=i)
for(j=i;j<=n;j++)
{
p=*(c+i*(n+1)+j);
*(c+i*(n+1)+j)=*(c+k*(n+1)+j);
*(c+k*(n+1)+j)=p;
}
for(j=i+1;j<=n-1;j++)
{
p=(*(c+j*(n+1)+i))/(*(c+i*(n+1)+i));
for(t=i;t<=n;t++)
*(c+j*(n+1)+t)-=p*(*(c+i*(n+1)+t));
}
}
for(i=n-1;i>=0;i--)
{
for(j=n-1;j>=i+1;j--)
(*(c+i*(n+1)+n))-=x[j]*(*(c+i*(n+1)+j));
x[i]=*(c+i*(n+1)+n)/(*(c+i*(n+1)+i));
}
return x;
}
-------------------------------------------
高斯賽德爾迭代法（通用程序）：
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define N 100
float *GauseSeidel(float *a,int n)
{
int i,j,nu=0;
float *x,dx;
x=(float *)malloc(n*sizeof(float));
for(i=0;i<=n-1;i++)
x[i]=0.0;
do {
for(i=0;i<=n-1;i++) {
float d=0.0;
for(j=0;j<=n-1;j++)
d+=*(a+i*(n+1)+j)*x[j];
dx=(*(a+i*(n+1)+n)-d)/(*(a+i*(n+1)+i));
x[i]+=dx;
}
if(nu>=N)
{
printf("iter divergence\n");
exit(1);
}
nu++;
}
while(fabs(dx)>1e-6);
return x;
}
void main()
{
int i;
float *x;
float c[20]={4.,-1.,0.,2.,-1.0,-1,4,-1,0,-7.0,0,-1,4,-1,9.0,2,0,-1,4.0,0.0};
float *GauseSeidel(float *,int);
x=GauseSeidel(c,4);
for(i=0;i<=3;i++)
printf("x[%d]=%f\n",i,x[i]);
getch();
}