使用LU分解求解Ax=b給出inf值

我正在嘗試使用 A = LU 演算法求解形式為 Ax = b 的線性系統，其中 A 是實數的 (nxn) 矩陣和實數的 ba (1xn) 向量。這是我的實作：

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

int LUPDecompose(double A[N][N], double Tol, int P[N])
{

    int i, j, k, imax;
    double maxA, ptr[N], absA;

    for (i = 0; i <= N; i  )
        P[i] = i; //Unit permutation matrix, P[N] initialized with N

    for (i = 0; i < N; i  ) {
        maxA = 0.0;
        imax = i;

        for (k = i; k < N; k  )
            if ((absA = abs(A[k][i])) > maxA) {
                maxA = absA;
                imax = k;
            }

        if (maxA < Tol) return 0; //failure, matrix is degenerate

        if (imax != i) {
            //pivoting P
            j = P[i];
            P[i] = P[imax];
            P[imax] = j;

            //pivoting rows of A
            for (int ii = 0; ii < N; ii  )
            {
                ptr[ii] = A[i][ii];
                A[i][ii] = A[imax][ii];
                A[imax][ii] = ptr[ii];
            }

            //counting pivots starting from N (for determinant)
            P[N]  ;
        }

        for (j = i   1; j < N; j  ) {
            A[j][i] /= A[i][i];

            for (k = i   1; k < N; k  )
                A[j][k] -= A[j][i] * A[i][k];
        }
    }

    return 1;  //decomposition done
}

/* INPUT: A,P filled in LUPDecompose; b - rhs vector; N - dimension
 * OUTPUT: x - solution vector of A*x=b
 */
void LUPSolve(double A[N][N], int P[N], double b[N], double x[N])
{

    for (int i = 0; i < N; i  ) {
        x[i] = b[P[i]];

        for (int k = 0; k < i; k  )
            x[i] -= A[i][k] * x[k];
    }

    for (int i = N - 1; i >= 0; i--) {
        for (int k = i   1; k < N; k  )
            x[i] -= A[i][k] * x[k];

        x[i] /= A[i][i];
    }
}

int main()
{
    double Am[N][N] = {{0.6289, 0, 0.0128, 0.3184, 0.7151},
                        {0,   1,   0,   0,   0},
                        {0.0128, 0, 0.0021, 0.0045, 0.0380},
                        {0.3184, 0, 0.0045, 0.6618, 0.3371},
                        {0.7151, 0, 0.0380, 0.3371, 1.1381}};
    double bm[N] = {1.6752, 0, 0.0574, 1.3217, 2.2283};
    int Pm[N] = {0};
    double X[N] = {0};

    LUPDecompose( Am, 0.0001, Pm);
    LUPSolve(Am, Pm,  bm, X);

    printf("%f %f %f %f %f",X[0],X[1],X[2],X[3],X[4]);
}

但是，我正在獲得 inf 值。

-1.#IND00 -1.#IND00 3.166387 0.849298 0.670689

不知道是代碼問題還是演算法問題。任何幫助解決這個問題？

uj5u.com熱心網友回復：

“不知道是代碼問題還是演算法問題。有什么可以幫助解決這個問題的嗎？”

我相信存在代碼和演算法問題。以下是您的代碼，并進行了更正以僅解決編譯錯誤和警告（請參閱內嵌注釋）。它沒有在 C 語法之外進行除錯，以實作干凈的編譯，并且沒有錯誤地運行。（即運行時沒有被零除或inf錯誤。）

#define N 5 //required to be 5 by hard-coded array definitions in main()

int LUPDecompose(double A[N][N], double Tol, int P[N])
{

    int i, j, k, imax, ii;//added ii here to increase scope below
    double maxA, ptr[N], absA;

    //for (i = 0; i <= N; i  )
    for (i = 0; i < N; i  )
        P[i] = i; //Unit permutation matrix, P[N] initialized with N (actually init with i)

    for (i = 0; i < N; i  ) {
        maxA = 0.0;
        imax = i;

        for (k = i; k < N; k  )
            if ((absA = fabs(A[k][i])) > maxA) {// using fabs, not abs to avoid conversion of double to int.
                maxA = absA;
                imax = k;
            }

        if (maxA < Tol) return 0; //failure, matrix is degenerate

        if (imax != i) {
            //pivoting P
            j = P[i];
            P[i] = P[imax];
            P[imax] = j;

            //pivoting rows of A
           //for (int ii = 0; ii < N; ii  )
           for ( ii = 0; ii < N; ii  )
            {
                ptr[ii] = A[i][ii];
                A[i][ii] = A[imax][ii];
                A[imax][ii] = ptr[ii];
            }

            //counting pivots starting from N (for determinant)
            //P[N]  ;//N will always overflow for array with only N elements
            P[ii-1]  ;//use index here instead
        }

        for (j = i   1; j < N; j  ) {
            A[j][i] /= A[i][i];

            for (k = i   1; k < N; k  ) {//extra brackets added for readability
                A[j][k] -= A[j][i] * A[i][k];
            }
        }
    }

    return 1;  //decomposition done
}

/* INPUT: A,P filled in LUPDecompose; b - rhs vector; N - dimension
 * OUTPUT: x - solution vector of A*x=b
 */
void LUPSolve(double A[N][N], int P[N], double b[N], double x[N])
{

    for (int i = 0; i < N; i  ) {
        x[i] = b[P[i]];
        for (int k = 0; k < i; k  ) {//extra brackets added for readability
            x[i] -= A[i][k] * x[k];
        }
    }

        for (int i = N - 1; i >= 0; i--) {
            for (int k = i   1; k < N; k  ) {//additional brackets added for readability
                x[i] -= A[i][k] * x[k];
            }

            x[i] /= A[i][i];
        }
}
    

//int main()
int main(void)//minimum signature for main includes void
{
    //Note hardcoded arrays in this code require N == 5 (#define at top)
    double Am[N][N] = {{0.6289, 0, 0.0128, 0.3184, 0.7151},
                        {0,   1,   0,   0,   0},
                        {0.0128, 0, 0.0021, 0.0045, 0.0380},
                        {0.3184, 0, 0.0045, 0.6618, 0.3371},
                        {0.7151, 0, 0.0380, 0.3371, 1.1381}};
    double bm[N] = {1.6752, 0, 0.0574, 1.3217, 2.2283};
    int Pm[N] = {0};
    double X[N] = {0};

    LUPDecompose( Am, 0.0001, Pm);
    LUPSolve(Am, Pm,  bm, X);

    printf("%f %f %f %f %f",X[0],X[1],X[2],X[3],X[4]);
    
    return 0;  //int main(void){...} requires return statement.
}  

基于

正確的解決辦法是：

-0.590174531351002
0
-19.76923076923077
1.0517711171662125
2.6772727272727272

但是上面代碼的實際輸出是：

演算法相關的除錯由您來執行。

標籤：C 求解器

上一篇：我如何在c中撰寫一個程式，顯示不允許名字不超過10個字符？

下一篇：具有內部鏈接的變數的暫定定義具有不完整的非陣列型別：符合要求的實作表現出不同的行為