Ý tưởng: Giải hệ A.X = B (1) đk: Det(A) <> 0.
- Dùng thuật toán phân rã LU phân tích A = L.U, với L là ma trận tam giác dưới, U là ma trận tam giác trên.
- Hệ (1) tương đương L.U.X = B (2)
- Đặt U.X = Y (3), hệ (2) tương đương L.Y = B (4).
- Giải hệ (4) tương ứng với hệ tam giác dưới ta tìm được Y.
- Có Y giải hệ (3) tương ướng với hệ tam giác trên ta tìm được X.
Bài này code cho dữ liệu cụ thể, các bạn muốn nhập dữ liệu từ bàn phím thì sửa lại 1 tí.
- Dùng thuật toán phân rã LU phân tích A = L.U, với L là ma trận tam giác dưới, U là ma trận tam giác trên.
- Hệ (1) tương đương L.U.X = B (2)
- Đặt U.X = Y (3), hệ (2) tương đương L.Y = B (4).
- Giải hệ (4) tương ứng với hệ tam giác dưới ta tìm được Y.
- Có Y giải hệ (3) tương ướng với hệ tam giác trên ta tìm được X.
Code:
#include "math.h"
#include "conio.h"
#include "iostream.h"
#define max 4
/*Nhap ma tran he so*/
void NhapMaTran(float A[max][max], int m, int n){
for(int i = 0; i<m; i++)
for(int j = 0; j<n; j++) {
cout<<"a["<<i<<"]["<<j<<"] = ";
cin>>A[i][j];
}
}
/*Xuat ma tran*/
void XuatMaTran(float A[max][max], int m, int n) {
cout<<"\n";
for(int i=0 ; i<m; i++){
cout<<endl;
for(int j=0 ; j<n; j++)
cout<<A[i][j]<<"\t";
}
}
/*Xuat nghiem*/
void XuatNghiem(float X[], int n){
cout<<"\nNghiem cua he PTTT";
for(int i=0; i<n; i++)
cout<<"\nx"<<i+1<<"="<<X[i];
}
/*He duong cheo*/
char HeDuongCheo(float A[max][max], float X[max], float B[max], int n ) {
for(int i = 0; i<n; i++) {
if (A[i][i] != 0)
X[i] = B[i]/A[i][i];
else
return 0;
}
return 1;
}
/*He tam giac duoi*/
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;
}
/*He tam giac tren*/
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;
}
/*Phan ra A = LU*/
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];
}
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];
PhanRaLU(A,L,U,n);
HeTamGiacDuoi(L,Y,B,n);
HeTamGiacTren(U,X,Y,n);
}
/*Chuong trinh chinh*/
void main(){
int m,n;
float A[4][4] ={{2,3,1,5},{6,13,5,19},{2,19,10,23},{4,10,11,31}};
float B[4] = {1,1,1,1};
float X[4];
float L[max][max], U[max][max];
clrscr();
GiaiHePTTT(A,X,B,4);
XuatNghiem(X,4);
getch();
}
Bài này code cho dữ liệu cụ thể, các bạn muốn nhập dữ liệu từ bàn phím thì sửa lại 1 tí.