【線性代數】高斯消元與矩陣求逆


今天講了線性代數,順帶復習了一下之前沒有認真學的高斯消元以及矩陣求逆。

高斯消元:

考慮一個滿秩的系數矩陣,它意味着有唯一解;而不存在唯一解的充要條件就是其行列式為 \(0.\)

那么考慮如何求解方程組:用初等行變換的形式將矩陣消成上三角矩陣,從而我們得到了最后一個未知數的解,再進行回代即可。

也可以直接消成一個對角矩陣。這個地方是和求逆一一樣的。

\(O(n^3)\) 可以完成高斯消元的過程,其中注意要確定用每次數絕對值最大的來消元,這樣可以最大限度避免小數類計算的誤差。

下面考慮矩陣求逆:

設已知矩陣為 \(A,\) 目前要求 \(P=A^{-1}.\)

\[A\cdot A^{-1}=AP=I \]

那么如果我們通過 初等行變換 的形式,將 \(A\) 消為 \(I,\) 那么對左側的 \(I\) 進行同樣的初等變換,對應地,它運算后的矩陣就應該是 \(P\) 了。

於是同樣在 \(O(n^3)\) 的復雜度完成。

消元:

#include<bits/stdc++.h>
using namespace std;
const int N=105;
int n;
const double eps=1e-9;
double ans[N];
struct Matrix{
    double a[N][N];
    Matrix(){memset(a,0xfe,sizeof a);}
    void Guess(){
        for(int i=0;i<n;++i){
            int pos=i;
            for(int j=i+1;j<n;++j)
                if(fabs(a[pos][i])<fabs(a[j][i]))
                    pos=j;
            if(fabs(a[pos][i])<eps){
                puts("No Solution");
                exit(0);
            }
            if(i!=pos)swap(a[i],a[pos]);
            double div=a[i][i];
            for(int j=i;j<=n;++j)a[i][j]/=div;
            for(int j=i+1;j<n;++j){
                div=a[j][i];
                for(int k=i;k<=n;++k)a[j][k]-=div*a[i][k];
            }
        }
    }
}A;

int main(){
    scanf("%d",&n);
    for(int i=0;i<n;++i)
        for(int j=0;j<=n;++j)
            scanf("%lf",&A.a[i][j]);
    A.Guess();
    ans[n-1]=A.a[n-1][n];
    for(int i=n-2;~i;--i){
        ans[i]=A.a[i][n];
        for(int j=i+1;j<n;++j)ans[i]-=ans[j]*A.a[i][j];
    }
    for(int i=0;i<n;++i)printf("%.2lf\n",ans[i]);
    return 0;
}

求逆:

#include<bits/stdc++.h>
using namespace std;
const int N=500;
const int mod=1e9+7;
int n;
inline int Add(int x,int y){return (x+y)%mod;}
inline int Mul(int x,int y){return 1ll*x*y%mod;}
inline int Dec(int x,int y){return (x-y+mod)%mod;}
inline int Min(int x,int y){return x<y?x:y;}
inline int Max(int x,int y){return x>y?x:y;}
inline int Abs(int x){if(x<0)x=-x;return x;}
inline int qpow(int x,int y){
    int res=1;
    while(y){
        if(y&1)res=Mul(res,x);
        x=Mul(x,x);y>>=1;
    }
    return res;
}
inline int getinv(int x){return qpow(x,mod-2);}
struct Matrix{
    int a[N][N];
    Matrix(){memset(a,0,sizeof a);}
    void I(){for(int i=0;i<N;++i)a[i][i]=1;}
    void print(){
        for(int i=0;i<n;++i)
            for(int j=0;j<n;++j)
                printf("%d%c",a[i][j],j==(n-1)?'\n':' ');
    }
    void Init(){
        for(int i=0;i<n;++i)
            for(int j=0;j<n;++j)
                scanf("%d",&a[i][j]);
    }
    Matrix Getinv(){
        Matrix res;
        res.I();
        for(int i=0;i<n;++i){
            int pos=i;
            for(int j=i+1;j<n;++j)
                if(Abs(a[pos][i])<Abs(a[j][i]))
                    pos=j;
            if(!a[pos][i]){
                puts("No Solution");
                exit(0);
            }
            if(i!=pos)swap(a[i],a[pos]),swap(res.a[i],res.a[pos]);
            int div=a[i][i];
            div=getinv(div);
            for(int j=0;j<n;++j){
                res.a[i][j]=Mul(res.a[i][j],div);
                a[i][j]=Mul(a[i][j],div);
            }
            for(int j=0;j<n;++j){
                if(i==j)continue;
                div=a[j][i];
                for(int k=0;k<n;++k){
                    res.a[j][k]=Dec(res.a[j][k],Mul(div,res.a[i][k]));
                    a[j][k]=Dec(a[j][k],Mul(div,a[i][k]));
                }
            }
        }
        return res;
    }
}A;
int main(){
    scanf("%d",&n);
    A.Init();
    Matrix Ans=A.Getinv();
    Ans.print();
    // A.print();
    return 0;
}

關於應用:待補


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM