高斯消元 - dblank

高斯消元

方法也很见简单,从第1行开始,我们利用当前行第i列不为0,就可以通过变换将i+1..M行第一列全部变换为0,接着对于第2行,我们用同样的方法将第3..M行第2列也变换为0...不断重复直到第n行为止。

假如计算到第i行时,第i列已经为0,则我们需要在第i+1..M行中找到一行第i列不为0的行k,并交换第i行和第k行,来保证ai != 0。但这时候还有可能出现一个情况,就是第i..M行中的i列均为0,此时可以判定,该方程组有多解。

当得到上三角矩阵后,就可以从第n行开始逆推,一步一步将a[i]j也变换为0.

因为第n行为an * x[n] = y'[n],则x[n] = y'[n] / an。

第n-1行为an-1 x[n - 1] + an x[n] = y'[n - 1]。我们将得到的x[n]代入,即可计算出x[n-1]。

同样的依次类推就可以得到所有的x[1]..x[n]。

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<map>
#include<cstring>
#include<map>
#include<set>
#include<queue>
#include<stack>
#include<cmath>
using namespace std;
const double eps = 1e-6;
const int inf = 0x3f3f3f3f, N = 1000 + 10, Mod = 1000000000;
typedef long long ll;
int n, m, ans[N];
double aN, b[N];
inline void swap_l(const int x, const int y)
{
    if(x!=y)
    {
        for(int i = 1; i<=n; i++)
            swap(ax, ay);
        swap(b[x], b[y]);
    }
}
inline bool zero(double x)
{
    return fabs(x) < eps;
}
int guass()
{
    for(int i = 1; i<=n; i++)
    {
        int r = 0;
        for(int j = i; j<=m; j++)// 找到第j列不为0的那一行
        {
            if(fabs(aj) > fabs(ar))
            {
                r = j;
                break;
            }
        }
        if(!r) return 2;
        swap_l(r, i);// 如果第i行的第i列为0,则与之前找到的那一行交换
        for(int j = i+1; j<=m; j++)//将 i+1 ~ m 行的 第i列全部消为0
        {
            double f = aj / ai;
            for(int k = 1; k<=n; k++)
                aj -= ai*f;
            b[j] -= f*b[i];
        }
    }
    for(int i = 1; i<=m; i++) // 检查是否多解或者无解
    {
        int flag = 1;
        for(int j = 1; j<=n; j++)
            if(!zero(ai))
                flag = 0;
        if(flag && zero(b[i])) flag = 0; // 如果某一行全为0,但是等式另一边为0,有唯一解
        if(flag) return 0; // 左边系数全为0,右边不为0,没有解
    }
    for(int i = n; i>=1; i--) // 对方程求解
    {
        for(int j = i+1; j<=n; j++)
        {
            b[i] -=  ai * ans[j];
            ai = 0;
        }
        ans[i] = b[i] / ai + 0.5;
    }
    return 1;
}
int main()
{
    scanf("%d%d", &n, &m);
    for(int i = 1; i<=m; i++)
    {
        for(int j = 1; j<=n; j++)
        {
            scanf("%lf", &ai);
        }
        scanf("%lf", &b[i]);
    }
    int flag = guass();
    if(flag == 2)
        printf("Many solutionsn");
    else if(flag == 0)
        printf("No solutionsn");
    else
    {
        for(int i = 1; i<=n; i++)
        {
            printf("%dn", ans[i]);
        }
    }
    return 0;
}

利用高斯消元还可以求解一类开关问题,例如http://POJ 1681,

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<map>
#include<cstring>
#include<map>
#include<set>
#include<queue>
#include<stack>
#include<cmath>
using namespace std;
const double eps = 1e-6;
const int inf = 0x3f3f3f3f, N = 300 + 10, Mod = 1000000000;
int equ, var;//有equ个方程,有var个待解
int aN, x[N], free_x[N], free_num;//分别是增广矩阵,解集,自由元,自由元的个数
int guass()
{
    int max_r, col, k;
    free_num = 0;
    for(k = 0, col = 0; k<equ && col <  var; k++, col ++)
    {
        max_r = k;
        for(int i = k+1; i<equ; i++)
        {
            if(abs(ai) > abs(amax_r))//找到一行第col列不为0的max_r
                max_r = i;
        }
        if(amax_r == 0)//存在自由元
        {
            k--;
            free_x[free_num++] = col;
            continue;
        }
        if(max_r != k)//交换这个两行
        {
            for(int j = col; j<var+1; j++)
                swap(ak, amax_r);
        }
        for(int i = k+1; i<equ; i++)//进行行变换,将k下面第col列不为0的那一行变换
        {
            if(ai)
            {
                for(int j = col; j<var + 1; j++)
                    ai ^= ak;
            }
        }
    }
    for(int i = k; i<equ; i++)//无解的情况
        if(ai)
            return -1;
    if(k<var) return var - k;
    for(int i = var - 1; i>=0; i--)//回代,并解出解集
    {
        x[i] = ai;
        for(int j = i+1; j<var; j++)
            x[i]^= (ai && x[j]);
    }
    return 0;
}
int n;
void init()//初始化,构造方程组
{
    memset(a, 0, sizeof(a));
    memset(x, 0, sizeof(x));
    equ = n*n;
    var = n*n;
    for(int i = 0; i<n; i++)
        for(int j = 0; j<n; j++)
        {
            int t = i * n + j;
            at = 1;
            if(j>0) at-1 = 1;
            if(i>0) at-n = 1;
            if(i<n-1) at+n = 1;
            if(j<n-1) at+1 = 1;
        }
}
int solve()
{
    int t = guass();
    if(t==-1)
    {
        printf("infn");
    }
    else if(t==0)
    {
        int ans = 0;
        for(int i = 0; i<var; i++)
            ans += x[i];
        printf("%dn", ans);
    }
    else//枚举自由元,找到最小的解
    {
        int ans = inf;
        int tot = (1<<t);
        for(int i = 0; i<tot; i++)
        {
            int cnt = 0 ;
            for(int j = 0; j<t; j++)
            {
                if(i&(1<<j))
                {
                    x[free_x[j]] = 1;
                    cnt++;
                }
                else x[free_x[j]] = 0;
            }
            for(int j = var - t - 1; j>=0; j--)
            {
                int idx;
                for(idx = j; idx < var ; idx++)
                    if(aj)
                        break;
                x[idx] = aj;
                for(int l = idx + 1; l<var; l++)
                    if(aj)
                        x[idx] ^= x[l];
                cnt += x[idx];
            }
            ans = min(cnt, ans);
        }
        printf("%dn", ans);
    }
}
int main()
{
    char str50;
    int T;
    scanf("%d", &T);
    for(int cas = 1; cas<=T ; cas++)
    {
        scanf("%d", &n);
        init();
        for(int i = 0; i<n; i++)
        {
            scanf("%s", str[i]);
            for(int j = 0; j<n; j++)
            {
                if(stri == 'y')
                    ai*n+j = 0;
                else ai*n+j = 1;
            }
        }
        solve();
    }
    return 0;
}

 

发表新评论