P3389 【模板】高斯消元法

題目背景

Gauss消元

題目描述

給定一個線性方程組，對其求解

輸入輸出格式

輸入格式：

第一行，一個正整數 n

第二至 n+1行，每行 n+1 個整數，為a1​,a2​⋯an​ 和 b，代表一組方程。

輸出格式：

共n行，每行一個數，第 i行為 xi​ （保留2位小數）

如果不存在唯一解，在第一行輸出"No Solution".

輸入輸出樣例

輸入樣例#1：

3

1 3 4 5

1 4 7 3

9 3 2 2
輸出樣例#1：

-0.97

5.18

-2.39

說明

1≤n≤100,∣ai​∣≤104,∣b∣≤104

sol：究竟是什麼東西（感覺像小學奧數）

我們首先確定一個方程組作為例子

x-2y+3z=6

4x-5y+6z=12

7x-8y+10z=21

先將它轉化為矩陣

1 -2 3 6

4 -5 6 12

7 -8 10 21

解決這個方程組

我們會希望它變成如下形式

1 0 0 a

0 1 0 b

0 0 1 c

這樣就可以表示為x=a,y=b,z=cx=a,y=b,z=c

我們使用高斯消元，就要一步一步將每個未知數約去。

這種方法是以列為單位消去的

首先我們將第一列轉化為1 0 0的形式

在這裡要注意一下，我們往往是將這個係數絕對值最大的方程轉移到被減的這一行，這樣就可以減小誤差

所以我們先將矩陣變成這樣

7 -8 10 21

4 -5 6 12

1 -2 3 6

然後將正在處理的方程式化簡，讓正被處理的係數化1

1 -8/7 10/7 3

4 -5 6 12

1 -2 3 6

然後使用加減法將第二個與第三個方程組的第一個係數化0

1 -8/7 10/7 3

0 -3/7 2/7 0

0 -6/7 11/7 3

然後這時候第一列就被化簡完成

同理我們化去第二行與第三行，步驟如下：

1.化簡第二行

1 -8/7 10/7 3

0 1 -2/3 0

0 -6/7 11/7 3

2.用第一行減第二行×(-8/7),第三行減第二行×(-6/7)

1 0 2/3 3

0 1 -2/3 0

0 0 1 3

3.不需要化簡第三行，所以直接用第一行減去第三行×2/3，第二行減去第三行×(-2/3)

1 0 0 1

0 1 0 2

0 0 1 3

最後我們就得到了一組解x=1,y=2,y=3x=1,y=2,y=3。所以高斯消元其實是運用了小學解方程組的加減法的呢。

注意是一列列消去的

#include <bits/stdc++.h>

using namespace std;

typedef int ll;

inline ll read()

{

    ll S=;

    bool f=;

    char ch=' ';

    while(!isdigit(ch))

    {

        f|=(ch=='-'); ch=getchar();

    }

    while(isdigit(ch))

    {

        S=S*+(ch-''); ch=getchar();

    }

    return (f)?(-S):(S);

}

#define R(x) x=read()

inline void write(ll x)

{

    if(x<)

    {

        putchar('-'); x=-x;

    }

    if(x<)

    {

        putchar(x+''); return;

    }

    write(x/);

    putchar(x%+'');

    return;

}

#define W(x) write(x),putchar(' ')

#define Wl(x) write(x),putchar('\n')

const int N=;

const double eps=1e-;

int n;

double a[N][N];

inline bool Gauss()

{

    int i,j,k;

    for(i=;i<=n;i++)

    {

        k=i;

        for(j=i+;j<=n;j++) if(fabs(a[j][i])>fabs(a[k][i])) k=j;

        //找到最大的數

        if(fabs(a[k][i])<eps) return false;

        if(i!=k) for(j=i;j<=n+;j++) swap(a[k][j],a[i][j]);

        //對換一行或一列,屬於找最大當前係數的其中一步.(這樣就可以只處理當前行的係數啦!)

        double Div=a[i][i];

        for(j=i;j<=n+;j++) a[i][j]/=Div;

        for(j=;j<=n;j++) if(j!=i)

        {

            Div=a[j][i];

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

            {

                a[j][k]-=a[i][k]*Div;

            }

        }

        //把這列除這個數外都消成0

    }

    return true;

}

int main()

{

    int i,j;

    R(n);

    for(i=;i<=n;i++)

    {

        for(j=;j<=n+;j++) scanf("%lf",&a[i][j]);

    }

    if(!(Gauss()))

    {

        puts("No Solution");

    }

    else

    {

        for(i=;i<=n;i++) printf("%.2lf\n",a[i][n+]);

    }

    return ;

}

/*

input

3

1 3 4 5

1 4 7 3

9 3 2 2

output

-0.97

5.18

-2.39

*/

附上更加符合上面教程的程式碼（完全按照上面寫的）

#include <bits/stdc++.h>

using namespace std;

typedef int ll;

inline ll read()

{

    ll s=;

    bool f=;

    char ch=' ';

    while(!isdigit(ch))

    {

        f|=(ch=='-'); ch=getchar();

    }

    while(isdigit(ch))

    {

        s=(s<<)+(s<<)+(ch^); ch=getchar();

    }

    return (f)?(-s):(s);

}

#define R(x) x=read()

inline void write(ll x)

{

    if(x<)

    {

        putchar('-'); x=-x;

    }

    if(x<)

    {

        putchar(x+''); return;

    }

    write(x/);

    putchar((x%)+'');

    return;

}

#define W(x) write(x),putchar(' ')

#define Wl(x) write(x),putchar('\n')

const double eps=1e-;

const int N=;

int n;

double a[N][N],b[N];

inline void Debug()

{

    int i,j;

    for(i=;i<=n;i++)

    {

        for(j=;j<=n;j++) printf("%.2lf ",a[i][j]);

        printf("%.2lf",b[i]);

        puts("");

    }

    puts("");

}

inline void Gauss(int n)

{

    int i,j,k;

    double Div;

    for(i=;i<=n;i++)

    {

        int Pos=i;

        for(j=i+;j<=n;j++) if(fabs(a[Pos][i])<fabs(a[j][i])) Pos=j;

        if(fabs(a[Pos][i])<eps)

        {

            puts("No Solution");

            exit();

        }

        if(Pos!=i)

        {

            swap(a[i],a[Pos]); swap(b[i],b[Pos]);

        }

        Div=a[i][i];

        for(j=i;j<=n;j++) a[i][j]/=Div; b[i]/=Div;

        for(j=;j<=n;j++) if(j!=i)

        {

            Div=a[j][i];

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

            {

                a[j][k]-=Div*a[i][k];

            }

            b[j]-=Div*b[i];

        }

//        Debug();

    }

}

int main()

{

    int i,j;

    R(n);

    for(i=;i<=n;i++)

    {

        for(j=;j<=n;j++) scanf("%lf",&a[i][j]); scanf("%lf",&b[i]);

    }

    Gauss(n);

    for(i=;i<=n;i++) printf("%.2lf\n",b[i]);

    return ;

}

/*

input

3

1 3 4 5

1 4 7 3

9 3 2 2

output

-0.97

5.18

-2.39

*/