[C\C++] Giải hệ phương trình tuyến tính dựa vào phân rã LU [Xử lý ma trận - mảng 2 chiều]

Ý tưởng giải thuật: cho hệ phương trình tuyến tính tổng quát A.X = B . Ta tiến hành phân rã A = L.U . Trong đó, L là ma trận tam giác dưới và U là ma trận tam giác trên.
Khi đó,

[Code Turbo C++]

#include "conio.h"
#include "iostream.h"
#define max 100

/*Nhap ma tran he so*/
void Nhap(float A[max][max],int n) {
 for(int i = 0; i<n; i++)
   for(int j = 0; j<n; j++) {
    cout<<"a["<<i<<"]["<<j<<"] = ";
    cin>>A[i][j];
  }
}

/*Nhap ma tran he so tu do*/
void Nhap(float B[max],int n) {
  for(int i = 0; i<n; i++) {
   cout<<"b["<<i<<"] = ";
   cin>>B[i];
 }
}

/*Xuat ma tran he so tu do*/
void Xuat(float B[max],int n) {
    cout<<"(";
    for(int i = 0; i<n-1; i++)
     cout<<B[i]<<",";
   cout<<B[n-1]<<")";
}

/*Xuat ma tran*/
void Xuat(float A[max][max], int n) {
  cout<<"\n";
  for(int i=0 ; i<n; i++){
   cout<<endl;
   for(int j=0 ; j<n; j++)
    cout<<A[i][j]<<"\t";
 }
}

/*Xuat nghiem*/
void XuatNghiem(float X[], int n, char * s) {
 cout<<"\nNghiem cua he PTTT";
 for(int i=0; i<n; i++)
   cout<<s<<i+1<<"="<<X[i];
 }

char HeTamGiacDuoi (float A[max][max], float X[max], float B[max], int n ) {
  for(int i = 0; i<n; i++) {
   if (A[i][i]!=0) {
     if (i==0)
      X[i] = B[i]/A[i][i];
    else {
     X[i] = B[i];
     for(int j=0; j<i; j++)
    X[i]=X[i]-A[i][j]*X[j];
    X[i] = X[i]/A[i][i];
   }
  } else
  return 0;
 }
 return 1;
}

char HeTamGiacTren (float A[max][max], float X[max], float B[max], int n ) {
  for(int i = n-1; i>=0; i--) {
    if (A[i][i]!=0) {
     if (i==n-1)
      X[i] = B[i]/A[i][i];
     else {
      X[i] = B[i];
      for(int j=i+1; j<n; j++)
        X[i]=X[i]-A[i][j]*X[j];
     X[i] = X[i]/A[i][i];
   }
  } else
  return 0;
 }
 return 1;
}

void PhanRaLU(float A[max][max], float L[max][max], float U[max][max], int n) {
  for(int k =0; k<n; k++) {
    U[k][k] = A[k][k];
    L[k][k] = 1;
    for(int i=k+1; i<n; i++) {
       L[i][k] = A[i][k]/U[k][k];
       U[k][i] = A[k][i];
       U[i][k] = 0;
       L[k][i] = 0;
   }
  for(i = k+1; i<n; i++)
    for(int j = k+1; j<n; j++)
      A[i][j] = A[i][j]-L[i][k]*U[k][j];
  }
}

/*Giai he phuong trinh tong quat LUX=B*/
void GiaiHePTTT(float A[max][max], float X[max], float B[max], int n) {
  float L[max][max],U[max][max], Y[max];
  cout<<"Phan ra A = L.U\n";
  PhanRaLU(A,L,U,n);
  cout<<"Ma tran L";
  Xuat(L,n);
  cout<<"\nMa tran U";
  Xuat(U,n);
  cout<<"\nGiai LY = B. Nghiem Y";
  if(HeTamGiacDuoi(L,Y,B,n)) {
    XuatNghiem(Y,n,"\ny");
    cout<<"\nGiai UX = Y. Nghiem X";
    if(HeTamGiacTren(U,X,Y,n))
      XuatNghiem(X,n,"\nx");
    else
      cout<<"\nHe pttt k co nghiem duy nhat";
   } else
     cout<<"\nHe pttt k co nghiem duy nhat";
  }
//////////////////////////
void main() {
 int n;
 float A[max][max], B[max], X[max];
 clrscr();
 cout<<"Nhap so he phuong trinh n = ";
 cin>>n;
 cout<<"Nhap ma tran he so A\n";
 Nhap(A,n);
 cout<<"Nhap ma tran he so tu do B\n";
 Nhap(B,n);
 GiaiHePTTT(A,X,B,n);
 getch();
}


[Tải code chương trình tại đây - Lưu ý: Sau 5s, Click Bỏ qua quảng cáo (Skin Ad)]