ガウス消去のロジックエラー

Question

ガウス消去コードのロジックエラーの問題...このコードは1990年の数値メソッドテキストのものでした。解決策は、必要がある

  SOLUTION OF SIMULTANEOUS LINEAR EQUATIONS 
     USING GAUSSIAN ELIMINATION 

This program uses Gaussian Elimination to solve the 
system Ax = B, where A is the matrix of known 
coefficients, B is the vector of known constants 
and x is the column matrix of the unknowns. 
Number of equations: 3 

Enter elements of matrix [A] 
A(1,1) = 0 
A(1,2) = -6 
A(1,3) = 9 
A(2,1) = 7 
A(2,2) = 0 
A(2,3) = -5 
A(3,1) = 5 
A(3,2) = -8 
A(3,3) = 6 

Enter elements of [b] vector 
B(1) = -3 
B(2) = 3 
B(3) = -4 

     SOLUTION OF SIMULTANEOUS LINEAR EQUATIONS 

     The solution is 
     x(1) = 0.000000 
     x(2) = -1.#IND00 
     x(3) = -1.#IND00 
     Determinant = -1.#IND00 
Press any key to continue . . . 

テキストからコピーしたように、コード...

//Modified Code from C Numerical Methods Text- June 2009 

#include <stdio.h> 
#include <math.h> 
#define MAXSIZE 20 

//function prototype 
int gauss (double a[][MAXSIZE], double b[], int n, double *det); 

int main(void) 
{ 
    double a[MAXSIZE][MAXSIZE], b[MAXSIZE], det; 
    int i, j, n, retval; 

    printf("\n \t SOLUTION OF SIMULTANEOUS LINEAR EQUATIONS"); 
    printf("\n \t USING GAUSSIAN ELIMINATION \n"); 
    printf("\n This program uses Gaussian Elimination to solve the"); 
    printf("\n system Ax = B, where A is the matrix of known"); 
    printf("\n coefficients, B is the vector of known constants"); 
    printf("\n and x is the column matrix of the unknowns."); 

    //get number of equations 
    n = 0; 
    while(n <= 0 || n > MAXSIZE) 
    { 
     printf("\n Number of equations: "); 
     scanf ("%d", &n); 
    } 

    //read matrix A 
    printf("\n Enter elements of matrix [A]\n"); 
    for (i = 0; i < n; i++) 
     for (j = 0; j < n; j++) 
     { 
      printf(" A(%d,%d) = ", i + 1, j + 1); 
      scanf("%lf", &a[i][j]); 
     } 
    //read {B} vector 
    printf("\n Enter elements of [b] vector\n"); 
    for (i = 0; i < n; i++) 
    { 
     printf(" B(%d) = ", i + 1); 
     scanf("%lf", &b[i]); 
    } 

    //call Gauss elimination function 
    retval = gauss(a, b, n, &det); 

    //print results 
    if (retval == 0) 
    { 
     printf("\n\t SOLUTION OF SIMULTANEOUS LINEAR EQUATIONS\n"); 
     printf("\n\t The solution is"); 
     for (i = 0; i < n; i++) 
      printf("\n \t x(%d) = %lf", i + 1, b[i]); 
     printf("\n \t Determinant = %lf \n", det); 
    } 
    else 
     printf("\n \t SINGULAR MATRIX \n"); 

    return 0; 
} 

/* Solves the system of equations [A]{x} = {B} using  */ 
/* the Gaussian elimination method with partial pivoting. */ 
/* Parameters:            */ 
/*  n   - number of equations     */ 
/*  a[n][n] - coefficient matrix      */ 
/*  b[n]  - right-hand side vector     */ 
/*  *det  - determinant of [A]      */ 

int gauss (double a[][MAXSIZE], double b[], int n, double *det) 
{ 
    double tol, temp, mult; 
    int npivot, i, j, l, k, flag; 

    //initialization 
    *det = 1.0; 
    tol = 1e-30;  //initial tolerance value 
    npivot = 0; 
     //mult = 0; 

    //forward elimination 
    for (k = 0; k < n; k++) 
    { 
     //search for max coefficient in pivot row- a[k][k] pivot element 
     for (i = k + 1; i < n; i++) 
     { 
      if (fabs(a[i][k]) > fabs(a[k][k])) 
      { 
       //interchange row with maxium element with pivot row 
       npivot++; 
       for (l = 0; l < n; l++) 
       { 
        temp = a[i][l]; 
        a[i][l] = a[k][l]; 
        a[k][l] = temp; 
       } 
       temp = b[i]; 
       b[i] = b[k]; 
       b[k] = temp; 
      } 
     } 
     //test for singularity 
     if (fabs(a[k][k]) < tol) 
     { 
      //matrix is singular- terminate 
      flag = 1; 
      return flag; 
     } 
     //compute determinant- the product of the pivot elements 
     *det = *det * a[k][k]; 

     //eliminate the coefficients of X(I) 
     for (i = k; i < n; i++) 
     { 
      mult = a[i][k]/a[k][k]; 
      b[i] = b[i] - b[k] * mult; //compute constants 
      for (j = k; j < n; j++)  //compute coefficients 
       a[i][j] = a[i][j] - a[k][j] * mult; 
     } 
    } 
    //adjust the sign of the determinant 
    if(npivot % 2 == 1) 
     *det = *det * (-1.0); 

    //backsubstitution 
    b[n] = b[n]/a[n][n]; 
    for(i = n - 1; i > 1; i--) 
    { 
     for(j = n; j > i + 1; j--) 
      b[i] = b[i] - a[i][j] * b[j]; 
     b[i] = b[i]/a[i - 1][i]; 
    } 
    flag = 0; 
    return flag; 
} 

：コードが正しい出力を生成しませbook-からで入力した...

サンプルを実行しています：1.058824、DETと1.823529、0.882353を-102.000000

任意の洞察力が高く評価されて...

Answer 1

//eliminate the coefficients of X(I) 
for (i = k; i < n; i++) 

これは多分それは今ある

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

方法でなければならず、私は、これはそれをゼロ、あなたはそれ自体で、ピボット行を分割するようになりますと信じています。

Answer 2

これはおそらくあなたが意図した方法であなたの質問に答えることはできませんが、あなた自身の数値的に安定したマトリックスアルゴリズムをプログラミングすることは、あなた自身の手術と同様に推奨されます。

C++の基本行列演算を行う信頼できるソース（NIST）のTNT/JAMAという非常に素晴らしいライブラリがあります。 Ax = B、最初のファクタAを解決するには（QR decompositionは良い一般的な方法ですが、LUを使用できますが数値的には安定しません）、solve(B)を呼び出します。これは、正確な（数値計算の問題の影響を受ける）正方行列と、最小自乗解を得る過決定システムの両方で機能します。