Giáo trình Phương pháp tính - Chương 3: Các vấn đề về ma trận

n = ");
	scanf("%d",&n);
	printf("Cho cac phan tu cua ma tran a : \n");
	for (i=1;i<=n;i++)
	 for (j=1;j<=n;j++)
	{
	 printf("a[%d][%d] = ",i,j );
	 scanf("%f",&a[i][j]);
	}
	printf("\n");
	clrscr();
	printf("Ma tran ban da nhap");
	printf("\n");
	for (i=1;i<=n;i++)
	 {
	for (j=1;j<=n;j++)
	 printf("%10.5f",a[i][j]);
	printf("\n");
	 }
	t=1;
	flushall();
	while (t)
	 {
	printf("\n");
	printf("Co sua ma tran khong(c/k)?");
	scanf("%c",&tl);
	if (toupper(tl)=='C')
	 {
	printf("Cho chi so hang can sua : ");
	scanf("%d",&i);
	printf("Cho chi so cot can sua : ");
	scanf("%d",&j);
	printf("a[%d][%d] = ",i,j);
	scanf("%f",&a[i][j]);
	flushall();
	 }
	if (toupper(tl)=='K')
	 t=0;
	 }
	printf("Ma tran ban dau");
	printf("\n");
	for (i=1;i<=n;i++)
	 {
	for (j=1;j<=n;j++)
	 printf("%10.5f",a[i][j]);
	printf("\n");
	 }
	for (i=1;i<=n;i++)
	 for (j=1;j<=n;j++)
	b[i][j]=a[i][j];
	for (k=1;k<=n-1;k++)
	 {
	vet=0.0;
	for (i=1;i<=n;i++)
	 vet+=b[i][i];
	p[k]=vet/k;
	for (i=1;i<=n;i++)
	 for (j=1;j<=n;j++)
	{
	 if (j!=i)
	c[i][j]=b[i][j];
	 if (j==i)
	c[i][j]=b[i][j]-p[k];
	}
	for (i=1;i<=n;i++)
	 for (j=1;j<=n;j++)
	{
	 b[i][j]=0.0;
	 for (k1=1;k1<=n;k1++)
	b[i][j]+=a[i][k1]*c[k1][j];
	}
	 }
	vet=0.0;
	for (i=1;i<=n;i++)
	 vet+=b[i][i];
	p[n]=vet/n;
	printf("\n");
	printf("Cac he so cua da thuc dac trung\n");
	printf("\n");
	d=1.0;
	printf("%6.2f",d);
	for (i=1;i<=n;i++)
	 {
	c1=-p[i];
	printf("%5c%6.2f",' ',c1);
	 }
	getch();
 }
2. Phương pháp Mises: Thuật toán Mises tìm giá trị riêng lớn nhất của một ma trận A. Nếu ma trận A là thực và và mỗi trị riêng bội k có đủ k vec tơ riêng độc lập tuyến tính thì việc tính toán sẽ cho ta giá trị riêng lớn nhất.
Một vectơ V bất kì có thể được viết dưới dạng:
	(5)
	Trong đó X1, X2,.., Xn là các vec tơ riêng tương ứng với các giá trị riêng l1, l2,l3,.., ln và v1, v2, v3,...,vn là các hằng số.
	Khi nhân A với V ta có:
	AV = Av1X1 + Av2X2 +....+ AvnXn
do:	Av1X1 = v1AX1 = v1l1X1 ; Av2X2 = v2AX2 = v2l2X2 v.v.
Nên:	AV = v1l1X1 + v2l2X2 +...+ vnlnXn
Lại nhân biểu thức trên với A ta có: 
 	A2V = v1l1 AX1 + v2l2 AX2 + ··· + vnlnAXn 
 = v1l21X1 + v2l22 X2 +...+ vnln2 Xn	
và tiếp đến lần thứ p ta có:
Lấy lp1 làm thừa số chung ta có:
Tương tự ta có:
Khi p rất lớn,vì l1 > l2 > l3 >,..., ln nên:	
Do đó:	(6)
nghĩa là khi p đủ lớn thì:
do đó: 
hay:	
Như vậylà véc tơ riêng của A ứng với l1 còn giá trị riêng l1 sẽ là:
	Trong thực tế để tránh vượt quá dung lượng bộ nhớ khi l1 khá lớn, các vectơ Vk được chuẩn hoá sau mỗi bước bằng cách chia các phần tử của nó cho phần tử lớn nhất mk và nhận được vectơ V’k 
Như vậy các bước tính sẽ là:
	- cho một vec tơ V bất kì (có thể là V = { 1, 1, 1,..., 1}T)
	- tính V1 = AV và nhận được phần tử lớn nhất là m1j từ đó tính tiếp V¢1 = V1/m1j
Một cách tổng quát, tại lần lặp thứ p ta nhận được vectơ Vp và phần tử lớn nhất mpj thì V’p = Vp/ mpj.
	- tính với vp+1,j là phần tử thứ j của Vp+1. Ta có:
Ví dụ: Tìm giá trị riêng lớn nhất và vec tơ riêng tương ứng của ma trận:
Chọn V= {1, 1, 1, 1}T ta tính được 
V
V1 = AV
V’1
V2 = AV’1
V’2
1
88
-0.6027
-6.4801
-0.5578
1
48
-0.3288
-5.6580
-0.4870
1
26
-0.1781
0.0818
0.0070
1
-146
1
11.6179
1
l
11.6179
V3 = AV’2
V’3
V4 = AV’3
V’4
V5 = AV’4
-3.9594
-0.5358
-3.6823
-0.5218
-3.5718
-3.6526
-0.4942
-3.5196
-0.4987
-3.4791
0.0707
0.0096
0.0630
0.0089
0.0408
7.3902
1
7.0573
1
6.9638
l
7.3902
7.0573
6.9638
V’5
V6= AV’5
V’6
V7= AV’6
V’7
-0.5129
-3.5341
-0.5075
-3.5173
-0.5043
-0.4996
-3.4809
-0.4999
-3.4868
-0.5000
0.0059
0.0250
0.0036
0.0147
0.0021
1
6.9634
1
6.9742
1
l
6.9634
6.9742
Dùng thuật toán trên ta có chương trình sau:
Chương trình 3-5
#include 
#include 
#include 
#include 
#include 
#define max 50

void main()
 {
	int i,j,k,n,t;
	char tl;
	float t0,t1,epsi,s;
	float a[max][max];
	float x0[max],x1[max];
	clrscr();
	printf("Phuong phap lap luy thua tim tri rieng lon nhat\n");
	printf("Cho so hang va cot cua ma tran n = ");
	scanf("%d",&n);
	printf("Cho cac phan tu cua ma tran a : \n");
	for (i=1;i<=n;i++)
	 for (j=1;j<=n;j++)
	{
	 printf("a[%d][%d] = ",i,j);
	 scanf("%f",&a[i][j]);
	}
	printf("\n");
	printf("Ma tran ban da nhap\n");
 printf("\n");
	for (i=1;i<=n;i++)
	 {
	for (j=1;j<=n;j++)
	 printf("%15.5f",a[i][j]);
	 printf("\n");
	 }
	flushall();
	t=1;
	while (t)
	 {
	printf("\nCo sua ma tran khong(c/k)?");
	scanf("%c",&tl);
	if (toupper(tl)=='C')
	 {
	printf("Cho chi so hang can sua : ");
	scanf("%d",&i);
	printf("Cho chi so cot can sua : ");
	scanf("%d",&j);
	printf("a[%d][%d] = ",i,j);
	scanf("%f",&a[i][j]);
	 }
	if (toupper(tl)=='K')
	 t=0;
	 }
	epsi=1e-5;
	printf("\nMa tran ban dau\n");
	printf("\n");
	for (i=1;i<=n;i++)
	 {
	for (j=1;j<=n;j++)
	 printf("%15.5f",a[i][j]);
	printf("\n");
	 }
	printf("\n");
	for (i=1;i<=n;i++)
	 x0[i]=1;
	k=1;
	t=0;
	t1=0;
	do
	 {
	t0=t1;
	for (i=1;i<=n;i++)
	 {
	x1[i]=0;
	for (j=1;j<=n;j++)
	 x1[i]=x1[i]+a[i][j]*x0[j];
	 }
	s=0;
	j=0;
	for (i=1;i<=n;i++)
	 if (s<fabs(x1[i]))
	{
	 j=i;
	 s=fabs(x1[i]);
	}
	t1=x1[j];
	for (i=1;i<=n;i++)
	 x1[i]=x1[i]/t1;
	if (fabs(t1-t0)<epsi)
	 {
	printf("Da thuc hien %d buoc lap\n",k);
	printf("Gia tri rieng lon nhat Vmax = %15.5f\n",t1);
	printf("Vec to rieng tuong ung\n");
	for (i=1;i<=n;i++)
	 printf("%.5f\n",x1[i]);
	t=1;
	 }
	if (fabs(t1-t0)>epsi)
	 {
	for (i=1;i<=n;i++)
	 x0[i]=x1[i];
	k=k+1;
	 }
	if (k>max)
	 t=1;
	 }
	while(t==0);
	getch();
 }
Dùng chương trình này tính gía trị riêng và vec tơ riêng của ma trận:
ta nhận được giá trị riêng là 3.0000 và vec tơ riêng là x = { -0.75 ; 0.75 ; 1 }T
	Như chúng ta đã nói trước đây, phương pháp Mises (hay còn gọi là phương pháp lặp lũy thừa) chỉ cho phép tìm giá trị riêng lớn nhất và vec tơ riêng tương ứng của ma trận. Để xác định các giá trị riêng khác, ma trận A được biến đổi thành một ma trận khác A1 mà các giá trị riêng là l2 > l3 >... Phương pháp này gọi là phương pháp xuống thang. Sau đây là phương pháp biến đổi ma trận: 
Giả sử X1 là vec tơ riêng của ma trận A tương ứng với giá trị riêng l1 và W1 là vec tơ riêng của ma trận AT tương ứng với giá trị riêng l1. Từ định nghĩa AX1 = l1X1 ta viết:
(A - lE)X1 = 0
Ta tạo ma trận A1 dạng: 
	(7)
Ta chú ý là X1W1T là một ma trận còn W1TX1 là một con số.Khi nhân hai vế của biểu thức (7) với X1 và chý ý đến tính kết hợp của tích các ma trận ta có:
	(8)
A1 chấp nhận giá trị riêng bằng không.
Nếu X2 là vec tơ riêng tương ứng với giá trị riêng l2,thì khi nhân A1 với X2 ta có: 
	(9)
Theo định nghĩa vì W1 là vectơ riêng của AT nên:
l1W1 =ATW1	(10)	
Mặt khác do:
(AX)T =XTAT và (AT)T = A
Nên khi chuyển vị (10) ta nhận được:
(ATW1)T = l1WT1	
Hay:	
W1TA = l1W1T	(11)
Khi nhân (11) với X2 ta có:	
l1W1TX2 = W1TAX2 
và do định nghĩa:
AX2 = l2X2 
nên:	
l1W1TX2 = W1T l2X2 
Vậy thì:	
(l1 - l2) W1TX2 = 0
khi l1 ¹ l2 thì:
	W1TX2 = 0	(12)
Cuối cùng thay (12) vào (9) ta có:
A1X2 = AX2 = l2X2
	Như vậy l2 là giá trị riêng lớn nhất của ma trận A1 và như vậy có thể áp dụng thuật toán này để tìm các giá trị riêng còn lại của ma trận. Các bước tính toán như sau 
	- khi đã có l1 và X1 ta tìm W1 là vec tơ riêng của AT ứng với giá trị riêng l1 (ví dụ tìm W1 bằng cách giải phương trình (AT -l1E)W1 = 0). Từ đó tính ma trận A12 theo (7).
	- tìm giá trị riêng và vec tơ riêng của A1 bằng cách lặp luỹ thừa và cứ thế tiếp tục và xuống thang (n-1) lần ta tìm đủ n giá trị riêng của ma trận A.
Ví dụ: Tìm giá trị riêng và vectơ riêng của ma trận sau:
Ta đã tìm được giá trị riêng lớn nhất l1 = 7 và một vectơ riêng tương ứng:
X1 = { 1, 1, 0, -2}T.
Ma trận AT có dạng:
và theo phương trình (AT -l1E)W1 = 0 ta tìm được vectơ W1 = {293,695,746,434}T
Ta lập ma trận mới A1 theo (7):
và:
Từ ma trận A1 ta tìm tiếp được l2 theo phép lặp luỹ thừa và sau đó lại tìm ma trận A3 và tìm giá trị riêng tương ứng.
Chương trình lặp tìm các giá trị riêng và vec tơ riêng của ma trận như sau:
Chương trình 3-6
#include 
#include 
#include 
#include 
#include 
#define max 50
void main()
 {
	float a[max][max],vv[max][max],at[max][max];
	float x[max],y[max],vd[max];
	int i,j,k,n,l,t;
	float vp,v1,z,epsi,va,ps;
	char tl;
	clrscr();
 epsi=0.000001;
	printf("Cho bac cua ma tran n = ");
	scanf("%d",&n);
	printf("Cho cac phan tu cua ma tran a : \n");
	for (i=1;i<=n;i++)
	 for (j=1;j<=n;j++)
	{
	 printf("a[%d][%d] = ",i,j);
	 scanf("%f",&a[i][j]);
	}
	printf("\n");
	clrscr();
	printf("Ma tran ban da nhap");
	printf("\n");
	for (i=1;i<=n;i++)
	 {
	for (j=1;j<=n;j++)
	 printf("%15.5f",a[i][j]);
	printf("\n");
	 }
	t=1;
	flushall();
	while (t)
	 {
	printf("\n");
	printf("Co sua ma tran khong(c/k)?");
	scanf("%c",&tl);
	if (toupper(tl)=='C')
	 {
	printf("Cho chi so hang can sua : ");
	scanf("%d",&i);
	printf("Cho chi so cot can sua : ");
	scanf("%d",&j);
	printf("a[%d][%d] = ",i,j);
	scanf("%f",&a[i][j]);
	 }
	if (toupper(tl)=='K')
	 t=0;
	 }
	for (l=1;l<=n;l++)
	 {
	for (i=1;i<=n;i++)
	 x[i]=1;
	vp=1.23456789;
	k=0;
	for (k=1;k<=40;k++)
	 {
	for (i=1;i<=n;i++)
	 {
	y[i]=0;
	for (j=1;j<=n;j++)
	 y[i]=y[i]+a[i][j]*x[j];
	 }
	v1=y[1]/x[1];
	z=0;
	for (i=1;i<=n;i++)
	 if (fabs(y[i])>z)
	z=y[i];
	for (i=1;i<=n;i++)
	 x[i]=y[i]/z;
	if (fabs(vp-v1)<epsi)
	 break;
	vp=v1;
	 }
	 {
	printf("Gia tri rieng : %9.6f\n",v1);
	printf("Vec to rieng : \n");
	 for (i=1;i<=n;i++)
	printf("%.5f\n",x[i]);
	printf("\n");
	getch();
	 }
	vd[l]=v1;
	va=v1;
	for (i=1;i<=n;i++)
	 vv[l][i]=x[i];
	for (i=1;i<=n;i++)
	 for (j=1;j<=n;j++)
	at[i][j]=a[j][i];
	for (i=1;i<=n;i++)
	 x[i]=1;
	vp=1.23456;
	k=0;
	for (k=1;k<=40;k++)
	 {
	for (i=1;i<=n;i++)
	 {
	y[i]=0;
	for (j=1;j<=n;j++)
	 y[i]=y[i]+at[i][j]*x[j];
	 }
	v1=y[1]/x[1];
	z=0;
	for (i=1;i<=n;i++)
	 if (fabs(y[i])>z)
	z=y[i];
	for (i=1;i<=n;i++)
	 x[i]=y[i]/z;
	if (fabs(vp-v1)<epsi)
	 break;
	vp=v1;
	 }
	if (fabs(vp-v1)>epsi)
	 {
	printf("Khong hoi tu sau 40 lan lap\n");
	getch();
	exit(1);
	 }
	if (fabs(va-v1)>3*epsi)
	 {
	printf("Co loi\n");
	getch();
	exit(1);
	 }
	ps=0;
	for (i=1;i<=n;i++)
	 ps=ps+x[i]*vv[l][i];
	ps=v1/ps;
	for (i=1;i<=n;i++)
	 for (j=1;j<=n;j++)
	a[i][j]=a[i][j]-ps*vv[l][i]*x[j];
	 }
 }
Dùng chương trình này tìm giá trị riêng của ma trận :
ta nhận được kết quả :
	giá trị riêng 	3.00000	vec tơ riêng 	
0.529411
	1.000000
	-0.411765
	giá trị riêng	3.000000	vec tơ riêng 	
-0.833336
	-0.166678
	1.000000
	giá trị riêng	-1.000000	vec tơ riêng	 
0.500000
	1.000000
	-0.500000
§5. PHÂN TÍCH MA TRẬN
1. Phương pháp Crout: Khi giải hệ phương trình tuyến tính nếu ta gặp một ma trận tam giác thì việc giải hệ sẽ rất dễ dàng. Vì vậy chúng ta tìm cách phân tích ma trận A thành tích của hai ma trận L và R sao cho : A = L.R . Để phân tích được, ma trận A phải có các giá trị trụ khác 0. Các ma trận L và R là các ma trận tam giác dưới (L) và tam giác trên (R).Các hệ số lkk = 1 . Ma trận L và R bậc 3 có dạng :
	Chúng ta nhắc lại quy tắc nhân hai ma trận A.B :
với 	c11= a11b11 + a12b21 + a13b31
	c12= a11b12 + a12b22 + a13b32
	c13= a11b13 + a12b23 + a13b33
	c21= a21b11 + a22b21 + a23b31
	Tổng quát :
	Dùng quy tắc này cho hai ma trận L và R và cho đồng nhất các hệ số của chúng với ma trận A ta có :
	a11 = 1. r11 + 0.0 + 0.0 = r11 ;
a12 = r12 ; a13 = r13 
	a21 = l21r11 ; 
a22 = l21r12 + r22 ; a23 = l31r11
	a31 = l31r11 ; a32 = l31r12 ; 
a33 = l31r13 + l32r23 + r33 
Một cách tổng quát ta có :
	với j > i : 	lij = rji = 0
	với i = 1 : 	r1j = a1j (j = 1 tới n)
 	lj1 = aj1/r11 (j = 1 tới n)
	với i = 2 tới n 	
 ( j = i tới n)
	 (j = i tới n)
	Chương trình phân tích ma trận thành 2 ma trận như sau :
Chương trình 3-7
#include 
#include 
#include 
#include 
#define max 6
void main()
 {
	float a[max][max],r[max][max],l[max][max];
	int i,j,k,n;
	float tr,tl;
	clrscr();
	printf("Cho bac cua ma tran n = ");
	scanf("%d",&n);
	printf("Cho cac phan tu cua ma tran can phan tich a\n");
	for (i=1;i<=n;i++)
	 for (j=1;j<=n;j++)
	{
	 printf("a[%d][%d] = ",i,j);
	 scanf("%f",&a[i][j]);
	}
	for (i=1;i<=n;i++)
	 for (j=1;j<=n;j++)
	{
	 l[i][j]=0.0;
	 r[i][j]=0.0;
	}
	for (i=1;i<=n;i++)
	 {
	r[1][i]=a[1][i];
	l[i][i]=1.0;
	l[i][1]=a[i][1]/a[1][1];
	 }
	for (k=2;k<=n;k++)
	 {
	for (j=k;j<=n;j++)
	 {
	tr=0.0;
	for (i=1;i<=k;i++)
	 tr=tr+l[k][i]*r[i][j];
	r[k][j]=a[k][j]-tr;
	 }
	if (k!=n)
	 {
	for (i=1;i<=n;i++)
	 {
	tl=0.0;
	for (j=1;j<=k-1;j++)
	 tl=tl+l[i][j]*r[j][k];
	l[i][k]=(a[i][k]-tl)/r[k][k];
	 }
	 }
	else
	 printf("\n");
	 }
	printf("Ma tran l :\n");
	for (i=1;i<=n;i++)
	 {
	for (j=1;j<=n;j++)
	 printf("%15.5f",l[i][j]);
	printf("\n");
	 }
	printf("Ma tran r :\n");
	for (i=1;i<=n;i++)
	 {
	for (j=1;j<=n;j++)
	 printf("%15.5f",r[i][j]);
	printf("\n");
	 }
	getch();
 }
Dùng chương trình này phân tích ma trận ta được :
2. Phương pháp phân tích Cholesky : Phương pháp Cholesky dùng để phân tích một ma trận đối xứng sao cho A = RTR với R là một ma trận tam giác trên. Cách phân tích cũng tương tự như phương pháp Crout . Ta xét các ma trận A và R bậc 3 như sau :
Tích hai ma trận RT và R là :
Ta tính được :
	r112 = a11
	r11r12 = a12
	r11r13 = a13
	r11r12 = a21
	r122 + r22r12 = a22
	r222 + r12r13 = a23
	r11r13 = a31
r13r12+ r23r21 = a32
r332 + r22r23 + r132 = a23
Tổng quát ta có : 
	rij = 0 (i > j )
Dưới đây là chương trình:
Chương trình 3-8
#include 
#include 
#include 
#include 
#define max 6
void main()
 {
	float a[max][max],r[max][max],b[max][max];
	int i,j,k,n,l;
	clrscr();
	printf("Cho bac cua ma tran n : ");
	scanf("%d",&n);
	printf("Cho cac phan tu cua ma tran can phan tich a :\n");
	for (i=1;i<=n;i++)
	 for (j=1;j<=n;j++)
	{
	 printf("a[%d][%d] = ",i,j);
	 scanf("%f",&a[i][j]);
	}
	for (i=1;i<=n;i++)
	 for (j=1;j<=n;j++)
	r[i][j]=0.0;
	for (i=1;i<=n;i++)
	 {
	if (a[i][i]<0.0)
	 {
	 printf("Ma tran khong duong");
	 getch();
	 exit(1);
	 }
	else
	 {
	r[i][i]=sqrt(a[i][i]);
	for (j=1+i;j<=n;j++)
	 r[i][j]=a[i][j]/r[i][i];
	for (k=i+1;k<=n;k++)
	 for (l=k;l<=n;l++)
	a[k][l]=a[k][l]-r[i][k]*r[i][l];
	 }
	}
	printf("\n");
	printf("Ma tran chuyen vi cua r\n");
	for (i=1;i<=n;i++)
	 for (j=1;j<=n;j++)
	b[j][i]=r[i][j];
	for (i=1;i<=n;i++)
	 {
	for (j=1;j<=n;j++)
	 printf("%15.5f",b[i][j]);
	printf("\n");
	 }
	printf("\n");
	printf("Ma tran r\n");
	for (i=1;i<=n;i++)
	 {
	for (j=1;j<=n;j++)
	 printf("%15.5f",r[i][j]);
	printf("\n");
	 }
	getch();
 }
Dùng chương trình này để phân tích ma trận
ta có :