初次做线性方程组的题。
可用高斯-约当消元法将原增广矩阵化为对角阵。需要注意无穷解、不定解等特殊情况。。
用一组数据解释一些特殊情况:
5
1 2
2 3
3 2
4 5
5 4
0 0
5
1 2 3 4 5
0
化简后的对角阵为:
1.000000 0.000000 0.000000 0.000000 0.000000 1.000000
0.000000 1.000000 -1.000000 0.000000 0.000000 1.000000
0.000000 0.000000 0.000000 0.000000 0.000000 1.000000
0.000000 0.000000 0.000000 1.000000 -1.000000 0.000000
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
明显答案为: f(1) = 1 , f(2) = f(3 ) = inf , f(4) = f(5 ) = 0;
判断方法为:
1. 无穷解:" m[ i ][ i ] = 0 && m[ i ][ n ] != 0 " || "m[ i ][ j ] != 0 && f( j ) = inf "
2. 不定解:"m[ i ][ i ] = m[ i ][ n ] = 0" ( 注意这题应输出0.000 ,看上面数据中的5)
3.其他情况有唯一解: f( i ) = m[ i ][ n ] / m[ i ][ i ];
参考代码:
#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
using namespace std;
#define MEM(a) memset(a,0,sizeof(a))
const int maxn = 110;
const double eps = 1e-8;
typedef double Mat[maxn][maxn];
int head[maxn] , next[maxn*maxn] , pre[maxn*maxn] ,en;
int inf[maxn] , d[maxn] , N;
Mat m;
void init(){
MEM(head) , MEM(m) ,MEM(d) ,MEM(inf);
en = 1;
}
void add_fa(int u,int v){
next[en] = head[u] , pre[en] = v , head[u] = en++;
}
void Debug(){
printf("\nMat:\n");
for(int i=0;i<N;i++){
for(int j=0;j<=N;j++){
printf("%f ",m[i][j]);
}
printf("\n");
}
}
void gause_jordan(Mat &m ,int n){
int i , j , k ,r;
double K;
for(i=0;i<n;i++){
r = i;
for(j=i+1 ;j<n;j++){
if(abs(m[j][i]) > abs(m[r][i])) r = j;
}
if(abs(m[r][i]) < eps) continue;
if(r!=i) for(int j=0;j<=n;j++) swap(m[i][j] , m[r][j]);
for(k=0;k<n;k++) {
if(k==i) continue;
K = m[k][i] / m[i][i];
for(j=n;j>=i;j--)
m[k][j] -= K*m[i][j];
}
}
}
int main()
{
int cas = 1;
while(scanf("%d",&N)!=EOF && N){
init();
int u , v;
do{
scanf("%d%d",&u,&v);
add_fa(v-1,u-1);
d[u-1]++;
}while(u+v);
m[0][N] = 1;
for(int i=0;i<N;i++){
m[i][i] = 1.0;
for(int p=head[i] ; p ; p=next[p]){
int f = pre[p];
m[i][f] -= 1.0 / d[f];
}
}
gause_jordan(m , N);
for(int i=N-1 ;i>=0;i--){
if(abs(m[i][i] < eps && abs(m[i][N]) > eps)) {
inf[i] = 1; continue;
}
for(int j=i+1;j<N;j++) {
if(abs(m[i][j]) > eps && inf[j]) {
inf[i] = 1; break;
}
}
}
//Debug();
int Q,q;
printf("Case #%d:\n",cas++);
scanf("%d",&Q);
while(Q--){
scanf("%d",&q); q--;
if(inf[q]) printf("infinity\n");
else printf("%.3f\n",abs(m[q][q]) < eps ? 0.0 : m[q][N] / m[q][q]);
}
}
return 0;
}