用C++解決線性方程組的程序

Question

我正在用這個簡單的2方程系統（矩陣形式爲「Mat [2] [3]」）來測試這個代碼來求解線性系統，但是當我執行它時，我得到下面的結果，不與我在系統中引入了矩陣的係數同意：

CODE：

//Gauss Elimination 
#include <iostream> 
#include <cmath> 
#include <vector> 
#include <fstream> 
#include <iterator> 
#include <string> 
#include <iomanip> 
using namespace std; 

int main() { 

    double Mat[2][3]; 
    vector<double> q(2); 
    int Nx = 2, ii, jj, kk; 

    Mat[1][1] = 2.0; 
    Mat[1][2] = 3.0; 
    Mat[1][3] = 2.5; 
    Mat[2][1] = 1.3; 
    Mat[2][2] = 3.0; 
    Mat[2][3] = 2.5;   

    cout << "Matrix: " << endl;  
    for (ii = 0; ii < Nx; ii++) { 
     for (jj = 0; jj < Nx + 1; jj++) { 
      cout << Mat[ii][jj] << " "; 
      cout << endl; 
     } 
    } 
    // Triangularization 
    for (ii = 0; ii < Nx - 1; ii++) 
     for (kk = ii + 1; kk < Nx; kk++) 
     { 
      double t = Mat[kk][ii]/Mat[ii][ii]; 
      for (jj = 0; jj <= Nx; jj++) 
       Mat[kk][jj] = Mat[kk][jj] - t * Mat[ii][jj]; 
     } // Resolution 
    for (ii = Nx - 1; ii >= 0; ii--) 
    { 
     q[ii] = Mat[ii][Nx]; 
     for (jj = Nx - 1; jj > ii; jj--)       
      q[ii] = q[ii] - Mat[ii][jj] * q[jj]; 
     q[ii] = q[ii]/Mat[ii][ii]; 
    } 

    cout << "Solution of the system: " << endl; 
    cout << q[1] << endl; 
    cout << q[2] << endl; //////////////////////////////////////////////////////////////////////// 


    return 0; 
} 

結果：

Matrix: 
0 2.07496e-317 6.95314e-310 
0 2 3 
Solution of the system: 
-nan 
0 

Answer 1

首先，你應該始終格式化您的代碼（這種類型現在爲你做了）。其次，C和C++中的索引從0開始，因此最大允許索引爲D-1，其中D是維度。這與例如MATLAB，其中索引從1開始。您的代碼中有無限制的訪問權限，例如Mat[1][3]=2.5;，因爲Mat被宣佈爲double Mat[2][3];，所以最大的行/列索引分別是1和2。 q顯示時相同，q[1]應該是q[0]和q[2]應該是q[1]。你的代碼會導致未定義的行爲。在開啓所有警告的情況下編譯代碼，即gcc上的-Wall -Wextra很可能會遇到這類錯誤。此外，請確保您的for循環也不會超出範圍。

作爲一個側面說明，你也可以直接初始化矩陣：

double Mat[2][3] = { {a,b,c}, {d,e,f} }; // where a, b etc are the coefficients 

Answer 2

我會強烈建議您通過這個網頁讀取。它使用高斯消元對線性方程進行了更深入的分析。https://martin-thoma.com/solving-linear-equations-with-gaussian-elimination/

您的代碼的主要問題是您的矩陣超出邊界訪問範圍。數組元素的索引開始於0，這意味着您的可訪問元素的索引爲0,1和2。

我也因此固定，在此編輯您的代碼：

#include <iostream> 

using namespace std; 

int main() 
{ 

const int Nx = 2; 
const int Ny = 3; 
double Mat[Nx][Ny]; 
double q[2]; 

Mat[0][0] = 2.0; 
Mat[0][1] = 3.0; 
Mat[0][2] = 2.5; 
Mat[1][0] = 1.3; 
Mat[1][1] = 3.0; 
Mat[1][2] = 2.5; 

cout << "Matrix: " << endl; 
for (int i = 0; i < Nx; i++) 
{ 
    for (int j = 0; j < Ny; j++) 
    { 
     cout << Mat[i][j] << " "; 
    } 
    cout << endl; 
} 
cout << endl; 

// Triangularization 
for (int i = 0; i < Nx - 1; i++) 
    for (int h = i + 1; h < Nx; h++) 
    { 
     double t = Mat[h][i]/Mat[i][i]; 
     for (int j = 0; j <= Nx; j++) 
     { 
      Mat[h][j] = Mat[h][j] - t * Mat[i][j]; 
     } 
    } 

// Resolution 
for (int i = Nx - 1; i >= 0; i--) 
{ 
    q[i] = Mat[i][Nx]; 
    for (int j = Nx - 1; j > i; j--) 
    { 
     q[i] = q[i] - Mat[i][j] * q[j]; 
    } 
    q[i] = q[i]/Mat[i][i]; 
} 

cout << "Solution of the system: " << endl; 
cout << q[0] << endl; 
cout << q[1] << endl; 


return 0; 
} 

使用Wolfram Alpha的計算給定矩陣的結果，結果是： this
而5/6約等於0.8(3)在Y軸上。我還沒有設法找到X軸的算法錯誤，因爲它返回-2.22045e-016，這是另一個界限或數學錯誤，但希望這會給你一個開始。