Gauss-Jordan解法方程式グループ
原文のリンク:http://blog.csdn.net/zhangxaochen/article/details/8020668
本のルーチン:
ここを見てくださいhttp://is.gd/VoBVUJ
自分でもう一度書きました。
1.行の一番大きい元を見つけて、主元として対角線上にあるかどうかを確認します。(1、3)でないなら、第一行と第三行を交換して、この主元を(3、3)に位置させます。
2.(1行目は3行目)位置要素を1にし、3行目の他の要素を同時に(3,3)の値で割ってください。マトリクスBの3行目の要素は、Aの(3,3)の値で同時に除算されます。
3.(3行目(3、3)はもう1つになりました。)他の行を等式3で減らして、2行目の3列を0にします。対角線にいないからです。
これはサイクル終了です。n個のdirtyまで全部trueをn回繰り返します。この時Aは単位行列です。マトリックスBは解マトリックスです。
個人の疑惑:なぜ一番大きな要素をメインとして選んだのか分かりません。個人的にはゼロ以外の要素であれば、メインを作ることができると思います。安定性を考えているという人がいますが、まだ分かりません。疑いを抱く
原文のリンク:http://blog.csdn.net/zhangxaochen/article/details/8020668
{OVER}
本のルーチン:
ここを見てくださいhttp://is.gd/VoBVUJ
void NR::gaussj(Mat_IO_DP &a, Mat_IO_DP &b){
int i,icol,irow,j,k,l,ll;
DP big,dum,pivinv;
int n=a.nrows();
int m=b.ncols();
Vec_INT indxc(n),indxr(n),ipiv(n);
for (j=0;j<n;j++) ipiv[j]=0;
for (i=0;i<n;i++) {
big=0.0;
for (j=0;j<n;j++){
if (ipiv[j] != 1){
for (k=0;k<n;k++) {
if (ipiv[k] == 0) {
if (fabs(a[j][k]) >= big) {
big=fabs(a[j][k]);
irow=j;
icol=k;
}
}
}
}
}
++(ipiv[icol]);
if (irow != icol) {
for (l=0;l<n;l++) SWAP(a[irow][l],a[icol][l]);
for (l=0;l<m;l++) SWAP(b[irow][l],b[icol][l]);
}
indxr[i]=irow;
indxc[i]=icol;
if (a[icol][icol] == 0.0) nrerror("gaussj: Singular Matrix");
pivinv=1.0/a[icol][icol];
//a[icol][icol]=1.0;
for (l=0;l<n;l++) a[icol][l] *= pivinv;
for (l=0;l<m;l++) b[icol][l] *= pivinv;
for (ll=0;ll<n;ll++){
if (ll != icol) {
dum=a[ll][icol];
//a[ll][icol]=0.0;
for (l=0;l<n;l++) a[ll][l] -= a[icol][l]*dum;
for (l=0;l<m;l++) b[ll][l] -= b[icol][l]*dum;
}
}
}
for (l=n-1;l>=0;l--) {
if (indxr[l] != indxc[l])
for (k=0;k<n;k++)
SWAP(a[k][indxr[l]],a[k][indxc[l]]);
}
}
上のコードは、とにかく分かりにくいです。二つのところを回って、わざとあなたをからかうようです。自分でもう一度書きました。
void myGaussj(double a[][3], double b[][1], int n, int m){
vector<bool> dirty(n);
double tmp;
int row, col;
for(int re=0; re<n; re++){
for(int i=0; i<n; i++){
if(true==dirty[i])
continue;
double big=0;
// :
for(int j=0; j<n; j++){
if( (tmp=fabs(a[i][j]))>big){
big=tmp;
row=i;
col=j;
}
}
assert(big!=0);
// big ,
dirty[col]=true; //
if(row!=col){
for(int c=0; c<n; c++){
tmp=a[row][c];
a[row][c]=a[col][c];
a[col][c]=tmp;
}
for(int c=0; c<m; c++){
tmp=b[row][c];
b[row][c]=b[col][c];
b[col][c]=tmp;
}
}// if row!=col
//
tmp=a[col][col];
for(int c=0; c<n; c++){
a[col][c]*=1.0/tmp;
}
for(int c=0; c<m; c++){
b[col][c]*=1.0/tmp;
}
//col ,
for(int r=0; r<n; r++){
if(r==col)
continue;
tmp=a[r][col];
for(int c=0; c<n; c++){
a[r][c]-=tmp*a[col][c];
}
for(int c=0; c<m; c++){
b[r][c]-=tmp*b[col][c];
}
}
break;
}// for i
}// for re
}
基本的な手順は、 行ごとに:1.行の一番大きい元を見つけて、主元として対角線上にあるかどうかを確認します。(1、3)でないなら、第一行と第三行を交換して、この主元を(3、3)に位置させます。
2.(1行目は3行目)位置要素を1にし、3行目の他の要素を同時に(3,3)の値で割ってください。マトリクスBの3行目の要素は、Aの(3,3)の値で同時に除算されます。
3.(3行目(3、3)はもう1つになりました。)他の行を等式3で減らして、2行目の3列を0にします。対角線にいないからです。
これはサイクル終了です。n個のdirtyまで全部trueをn回繰り返します。この時Aは単位行列です。マトリックスBは解マトリックスです。
個人の疑惑:なぜ一番大きな要素をメインとして選んだのか分かりません。個人的にはゼロ以外の要素であれば、メインを作ることができると思います。安定性を考えているという人がいますが、まだ分かりません。疑いを抱く
原文のリンク:http://blog.csdn.net/zhangxaochen/article/details/8020668
{OVER}