线性方程组的各种解法

题目：

用三种不同的方法计算线性方程组Ax=b，输入方程组的阶数n,矩阵A的元素和常向量b的元素,输出方程组的解。 一． 高斯列主元法 1．方法原理：.

高斯消去法包括两过程：先把方程组化为同解的上三角方程组，再按相反顺序求解上三角方程组。前者称消去或消元过程，后者称回代过程。

消去过程实际上是对增广矩阵作初等变换。对一般的n阶方程组，消去过程分n-1步：第一步消去a[1][1]下方元素，第二部消去a[2][2]下方元素，······，第n-1步消去a[n-1][n-1]下方元素。第k步将第k行的适当倍数加于其后各行，也可说是从k+1~n行减去第行的适当倍数，使它们的第k列元素变为0，而k+1~n+1列元素减去第k行对应列元素的倍数。

为了避免出现小主元，可在第k步的第k列的元素a[k][k],a[k+1][k],···,a[n][k]中选主元，即在其中找出绝对值最大的元素a[p][k]，然后交换第k和第p行，继续进行消去过程。交换行相当于改变方程顺序，不会影响原方程组的解。这种方法称为列主元消去法。 2. 源代码：

void Gaosi(float a[15][15],int n,float b[15],float x[15]) /*高斯列主元法的函数*/ {

int i,j,k;

for(k=1;k<=n-1;k++) {

float max=fabs(a[k][k]); int p_max=k;

for(int p=k;p<=n;p++) /*列主元的选主元*/ if(max

for(int m=1;m<=n+1;m++) {

float temp=0; temp=a[k][m]; a[k][m]=a[p][m]; a[p][m]=temp; } }

for(i=k+1;i<=n;i++) {

float c=a[i][k]/a[k][k]; for(j=k;j<=n+1;j++)

a[i][j]=a[i][j]-c*a[k][j]; } }

for(int p=1;p<=n;p++)

{

b[p]=a[p][n+1]; }

x[n]=b[n]/a[n][n]; for(k=n-1;k>=1;k--) {

float sum=0;

for(j=k+1;j<=n;j++) {

sum+=(a[k][j]*x[j]); }

x[k]=(b[k]-sum)/a[k][k]; }

cout<<\; cout<<\用高斯消去法解方程组：\<

for(int q=1;q<=n+1;q++) {

printf(\,a[m][q]); }

cout<

Output(x,n); /*依次输出方程的解*/ } 二． Jacobi迭代法 1. 方法原理：

由方程组的形式可构成迭代式：

x[i](k+1)=(-a[i][1]x[1](k)-···-a[i][i-1]x[i-1](k)-a[i][i+1]x[i+1](k)-···-a[1]x[n](k)+b[i])/a[i][i], i=1~n

给定初值x(0)=x(x[1](0),x[2](0),···x[n](0))T,令k=0,1,···，由此可得向量序列{x(k)}.如果此序列收敛于x,那么每个分量序列{x[i](k)}必收敛于想x[i],x[i](i=1~n)就必然是原方程组的解。此方法称为雅可比(Jacobi)迭代法。 2. 源代码：

void Jacobi(float a[15][15],float x[15],int n,float b[15]) /*雅可比迭代法函数*/ {

float y[15]={0}; double D;

for(int i=1;i

D=0;

for(int i=1;i

y[i]=b[i];

for(int j=1;j

y[i]=y[i]-a[i][j]*x[j]; y[i]=y[i]/a[i][i];

if(fabs(x[i]-y[i])>D) /*退出循环的条件*/ D=fabs(x[i]-y[i]); }

for(int i=1;i=Level) continue; else

break; }

cout<<\; cout<<\用Jacobi法解方程组：\<

Seidel迭代法是在Jacobi方法的基础上进行了改进。它在计算想x(k+1)时，将已经算出的分量立即代替x(k)对应分量，又称作逐个迭代法。 2. 源代码：

void Seidel(float a[15][15],float x[15],int n,float b[15]) /*Seidel迭代法函数*/ {

float y=0; double D;

for(int i=1;i

D=0;

for(int i=1;i

y=b[i];

for(int j=1;j

y=y-a[i][j]*x[j]; y=y/a[i][i];

if(fabs(x[i]-y)>D) D=fabs(x[i]-y); x[i]=y; }

if(D>=Level) continue; else

break; }

cout<<\; cout<<\用Seidel法解方程组：\<

cout<<\方程的解为:\<

cout<<\<

void Input(float a[15][15],int &n,float b[]) /*输入方程参数的函数*/ {

float m_in;

cout<<\请输入方程阶数n:\<>n;

cout<<\请依次输入方程对应的增广矩阵中的数据：\<

for(int q=1;q<=n+1;q++) {

cin>>m_in; a[m][q]=m_in; }

for(int p=1;p<=n;p++) {

b[p]=a[p][n+1]; } }

int _tmain(int argc, _TCHAR* argv[]) /*主函数*/ {

}

运行结果：

float a[15][15]={0},x[15]={0},b[15]={0}; int n=0; while(1) {

char in;

Input(a,n,b); Gaosi(a,n,b,x); Jacobi(a,x,n,b);

Seidel(a,x,n,b);

/*用三种方法解方程组*/

cout<<\再来解一个？y/n\<>in; if(in=='y') continue; else break; }

return 0;

五． 结果分析：

上方的方程组用三种方法得出的结果基本相同，但在计算其它的方程组时，会发现高斯

列主元法基本不会出大的差错，而两种迭代法会出现一些错误。这是由于高斯列主元法是一种比较通用的方法，而迭代法是要求方程组的解必须是收敛的。如果不收敛的话，算法会无限进行下去，程序也就运行不出结果。所以，在实际操作中，应该先判断方程组的解是否收敛，然后选择合适的算法进行计算。