第3章:线性方程组求解代码汇编问题:求Ax=b的解,A是M阶可逆方阵;
约定:算法中用到的是M×N增广矩阵,N=M+1;
变量:i,j,k等为整型变量,x,y,z为实型变量;
把任意M阶线性方程组化为上三角形方程组的C语言代码:
void ToTriangle()
{ int i,j,k;
float x;
for(k=0;k<M;k++)
{ x = A[k][k];
for(j=k;j<N;j++) A[k][j] /= x;
for(i=k+1;i<M;i++)
{ x=A[i][k];
for(j=k;j<N;j++)A[i][j] -= x*A[k][j];
}
}
return;
}
求解上三角形线性方程组c语言代码:
for(k=M-1;k>0;k--)
for(i=0;i<k;i++)
{ A[i][N-1] -= A[i][k]*A[k][N-1];
A[i][k]=0.0;
}
3.约当消去法C语言代码:
for(K=0;K<M;K++)
{ x=A[K][K];
for(j=K;j<N;j++) A[K][j] /= x;
for(i=0;i<M;i++)
{ if(i==K) continue;
x=A[i][K];
for(j=K;j<N;j++)A[i][j] -=x*A[K][j];
}
}
4.利用约当消去法求矩阵的逆对于M×M可逆矩阵A,我们可以构造矩阵T=(A I),那么用约当消去法把T的左边M列化为单位矩阵后,其效果相当于用A-1左乘T,从而得到A-1(A I)=(I A-1),所以T的后面的M列就是我们所要求的A-1
4实现选主元的源代码
j=K;
x=Math.abs(A[K][K]);
for(i=K+1;i<M;i++)
{ if(Math.abs(A[i][K])<x)
continue;
x=Math.abs(A[i][K]);
j=i;
}
//say("j="+j+",x="+A[j][K]);
for(i=K;i<N;i++)
{ x=A[K][i];
A[K][i]=A[j][i];
A[j][i]=x;
}
提示:把Math.abs()改写为fabs()即成为C语言代码。say()是用于调试的一个小脚本程序。
追过程c语言代码:
U[0]=D[0]/B[0];
V[0]=C[0]/B[0];
for(i=1;i<N;i++)
{ x=B[i]-A[i]*V[i-1];
U[i]=(D[i]- U[i-1]*A[i])/x;
V[i]=C[i]/x;
}
//赶过程c语言代码:
X[N-1]=U[N-1];
for(i=N-2;i>=0;i--)X[i]=U[i]-V[i]*X[i+1];