Các bài toán về ma trận và véc tơ

Trong mục này sẽ xét các ma trận thực vuông cấp n và các véc tơ thực cấp n. Chúng được biểu diễn thông qua các kiểu cấu trúc MT và VT:

struct MT
    {
    double a[20][20] ; // Mang a chứa các phần tử ma trận
    int n ; // Cấp ma trận
    } ;
    struct VT
    {
    double b[20]; // Mang chua cac phan tu cua vec to
    int n ; // Cap vec to
    } ;
    

Để xử lý ma trận và véc tơ, chúng ta xây dựng 9 hàm toán tử:

ostream& operator<< (ostream& os, const MT& x); // In ma trận
    ostream& operator<< (ostream& os, const VT& v); // In véc tơ
    istream& operator>> (istream& is,MT& x); // Nhập ma trận
    istream& operator>> (istream& is, VT &v); // Nhập véc tơ 
    MT operator+(const MT& x1, const MT& x2); // Cộng 2 ma trận
    MT operator-(const MT& x1, const MT& x2); // Trừ 2 ma trận
    MT operator*(const MT& x1, const MT& x2); // Nhân 2 ma trận
    VT operator*(const MT& x, const VT& v); // Nhân ma trận véc tơ
    MT operator!(MT x); // Nghịch đảo ma trận
    

Thuật toán cho 8 hàm toán tử đầu tương đối quen thuộc không có gì phải bàn. Để nghịch đảo ma trận có nhiều cách, ở đây chúng ta dùng phương pháp Jordance như sau. Giả sử cần nghịch đảo ma trận x cấp n. Ta dùng thêm ma trận đơn vị y. Sau đó thực hiện đồng thời các phép tính trên cả x và y sao cho x trở thành đơn vị. Kết quả y chính là nghịch đảo của x. Thuật toán được tiến hành trên n bước. Nội dung của bước k (k = 1,...,n) như sau:

Tìm chỉ số r ( k <= r <= n) sao cho

abs(x[r,k]) = max { abs(x[i,k] với i = k,...,n }

Nếu abs(x[r,k]) = 0 thì ma trận không có nghịch đảo và thuật toán kết thúc giữa chừng.

Hoán vị hàng k với hàng r trong cả 2 ma trận x và y.

Chia hàng k của cả x và y cho tg = x[k,k] (mục đích làm cho x[k,k] = 1).

Biến đổi để cột k của x trơ thành véc tơ đơn vị bằng cách làm cho các phần tử x[i,k] = 0 (với i khác k). Muốn vậy ta thực hiện các phép tính sau trên cả x và y:

(hàng i) = (hàng i) - x[i,k]*(hàng k) , với mọi i khác k

Nội dung chương trình là nhập 4 ma trận X, Y, R, S và véc tơ u. Sau đó tính véc tơ v theo công thức:

v = ((X + Y)*(R - S))-1u 

Như sẽ thấy trong hàm main() dưới đây, nhờ các hàm toán tử mà câu lệnh tính v được viết gần giống như công thức toán học nêu trên.

/* Chương trình */
    #include <conio.h>
    #include <iostream.h>
    #include <iomanip.h>
    #include <math.h>
    struct MT
    {
    double a[20][20]; // Mang chua cac phan tu ma tran
    int n ; // Cap ma tran
    } ;
    struct VT
    {
    double b[20]; // Mang chua cac phan tu cua vec to
    int n ; // Cap vec to
    } ;
    ostream& operator<< (ostream& os, const MT& x);
    ostream& operator<< (ostream& os, const VT& v);
    istream& operator>> (istream& is,MT& x);
    istream& operator>> (istream& is, VT &v);
    MT operator+(const MT& x1, const MT& x2);
    MT operator-(const MT& x1, const MT& x2);
    MT operator*(const MT& x1, const MT& x2);
    VT operator*(const MT& x, const VT& v);
    MT operator!(MT x); // Tinh ma tran nghich dao
    ostream& operator<< (ostream& os, const MT& x)
    {
    os << setprecision(2) << setiosflags(ios::showpoint);
    for (int i=1 ; i<= x.n ; ++i)
    {
    os << "
" ;
    for (int j=1; j<=x.n; ++j)
    os << setw(6) << x.a[i][j] ;
    }
    os << "
" ;
    return os;
    }
    ostream& operator<< (ostream& os, const VT& v)
    {
    os << setprecision(2) << setiosflags(ios::showpoint);
    for (int i=1 ; i<= v.n ; ++i)
    os << setw(6) << v.b[i] ;
    os << "
" ;
    return os;
    }
    istream& operator>> (istream& is, MT& x)
    {
    cout << " - Cap ma tran: " ;
    is >> x.n;
    cout << "Nhap cac phan tu :
" ;
    for (int i=1 ; i<= x.n ; ++i)
    for (int j=1; j<=x.n; ++j)
    {
    cout << "PT hang " << i << " cot " << j << " = " ;
    is >> x.a[i][j] ;
    }
    return is;
    }
    istream& operator>> (istream& is, VT& v)
    {
    cout << " - Cap vec to: " ;
    is >> v.n;
    cout << "Nhap cac phan tu :
" ;
    for (int i=1 ; i<= v.n ; ++i)
    {
    cout << "Phan tu thu " << i << " = " ;
    is >> v.b[i] ;
    }
    return is;
    }
    MT operator+(const MT& x1, const MT& x2)
    {
    if (x1.n!=x2.n)
    {
    cout << "
Khong thuc hien duoc phep cong vi 2 MT khong cung cap";
    getch();
    return x1;
    }
    else
    {
    MT x;
    int i, j, n;
    n = x.n = x1.n ;
    for (i=1; i<=n; ++i)
    for (j=1; j<=n ;++j)
    x.a[i][j] = x1.a[i][j] + x2.a[i][j] ;
    return x;
    }
    }
    MT operator-(const MT& x1, const MT& x2)
    {
    if (x1.n!=x2.n)
    {
    cout << "
Khong thuc hien duoc phep tru vi 2 MT khong cung cap";
    getch();
    return x1;
    }
    else
    {
    MT x;
    int i, j, n;
    n = x.n = x1.n;
    for (i=1; i<=n; ++i)
    for (j=1; j<=n ;++j)
    x.a[i][j] = x1.a[i][j] - x2.a[i][j] ;
    return x;
    }
    }
    MT operator*(const MT& x1, const MT& x2)
    {
    if (x1.n!=x2.n)
    {
    cout << "
Khong thuc hien duoc phep nhan vi 2 MT khong cung cap";
    getch();
    return x1;
    }
    else
    {
    MT x;
    int n, i, j,k;
    n = x.n = x1.n;
    for (i=1; i<=n; ++i)
    for (j=1; j<=n ;++j)
    {
    x.a[i][j] = 0.0 ;
    for (k=1 ; k<=n; ++k)
    x.a[i][j] += x1.a[i][k]*x2.a[k][j] ;
    }
    return x;
    }
    }
    VT operator*(const MT& x, const VT& v)
    {
    if (x.n != v.n)
    {
    cout << "
 Cap ma tran khac cap vec to, phep nhan vo nghia";
    getch();
    return v;
    }
    else
    {
    VT u; int n;
    n = u.n = v.n ;
    for (int i=1; i <=n ; ++i)
    {
    u.b[i] = 0;
    for (int j=1; j<=n; ++j)
    u.b[i] += x.a[i][j]*v.b[j];
    }
    return u;
    }
    }
    MT operator!(MT x)
    {
    MT y;
    int i,j,k,r,n;
    double tg;
    n = y.n = x.n ;
    for (i=1 ; i<=n ; ++i)
    for (j=1 ; j<=n ; ++j)
    if (i==j) y.a[i][j] = 1;
    else y.a[i][j] = 0;
    for (k=1; k<=n; ++k)
    {
    r=k;
    for (i=k+1; i<=n; ++i)
    if (abs(x.a[i][k]) > abs(x.a[r][k]) ) r = i;
    if (abs(x.a[r][k]) < 1.0E-8)
    {
    cout << "
 Ma tran suy bien, khong co nghich dao" ;
    getch();
    return x;
    }
    /* Hoan vi hang r va hang k */
    for (j=1 ; j<=n ; ++j)
    {
    tg = x.a[k][j];
    x.a[k][j] = x.a[r][j];
    x.a[r][j] = tg;
    tg = y.a[k][j];
    y.a[k][j] = y.a[r][j];
    y.a[r][j] = tg;
    }
    /* Chia hang k cho a[k,k] */
    tg = x.a[k][k] ;
    for (j=1 ; j<=n ; ++j)
    {
    x.a[k][j] /= tg;
    y.a[k][j] /= tg;
    }
    /* Khu cot k : lam cho a[i,k] = 0 voi i != k */
    for (int i=1; i<= n ; ++i)
    if (i != k)
    {
    tg = x.a[i][k] ;
    for (j=1 ; j<=n ; ++j)
    {
    x.a[i][j] -= tg*x.a[k][j] ;
    y.a[i][j] -= tg*y.a[k][j] ;
    }
    }
    }
    return y;
    }
    void main()
    {
    MT x,y,r,s;
    VT u,v;
    clrscr();
    cout <<"
Nhap ma tran X " ; cin >> x;
    cout <<"
Nhap ma tran Y " ; cin >> y;
    cout <<"
Nhap ma tran R " ; cin >> r;
    cout <<"
Nhap ma tran S " ; cin >> s;
    cout <<"
Nhap vec to u " ; cin >> u;
    v = !((x+y)*(r-s))*u ;
    cout << "
Vec to v = xu " << v ;
    getch();
    }