今天講了線性代數,順帶復習了一下之前沒有認真學的高斯消元以及矩陣求逆。
高斯消元:
考慮一個滿秩的系數矩陣,它意味着有唯一解;而不存在唯一解的充要條件就是其行列式為 \(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;
}
關於應用:待補