Giáo trình Turbo C nâng cao và C++ - Chương 9: Các vấn đề về ma trận

116 Ch−ơng 9 : Các vấn đề về ma trận Đ1.Định thức của ma trận Cho một ma trận vuông cấp n.Ta cần tìm định thức của nó.Tr−ớc hết chúng ta nhắc lại một số tính chất quan trọng của định thức: - nếu nhân tất cả các phần tử của một hàng (hay cột) với k thì định thức đ−ợc nhân với k - định thức không đổi nếu ta cộng thêm vào một hàng tổ hợp tuyến tính của các hàng còn lại. Ta sẽ áp dụng các tính chất này để tính định thức của một ma trận cấp 4 nh− sau(ph−ơng pháp này có thể mở rộng cho một ma trận cấp n) bằng ph−ơng pháp trụ: ⎟⎟ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎜⎜ ⎝ ⎛ = 44434241 34333231 24232221 14131211 aaaa aaaa aaaa aaaa A Lấy giá trị trụ là p1= a11.Ta chia các phần tử của hàng thứ nhất cho p1= a11 thì định thức sẽ là D/p1 (theo tính chất 1) và ma trận còn lại là: ⎟⎟ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎜⎜ ⎝ ⎛ ′′′ 44434241 34333231 24232221 141312 aaaa aaaa aaaa aaa1 Lấy hàng 2 trừ đi hàng 1 đã nhân với a21,lấy hàng 3 trừ đi hàng 1 đã nhân với a31 và lấy hàng 4 trừ đi hàng 1 đã nhân với a41 (thay hàng bằng tổ hợp tuyến tính của các hàng còn lại) thì định thức vẫn là D/p1 và ma trận là: ⎟⎟ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎜⎜ ⎝ ⎛ ′′′ ′′′ ′′′ ′′′ 444342 343332 242322 141312 aaa0 aaa0 aaa0 aaa1 Lấy giá trị trụ là 222 ap ′= .Ta chia các phần tử của hàng thứ hai cho p2 thì định thức sẽ là D/(p1p2) và ma trận còn lại là: ⎟⎟ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎜⎜ ⎝ ⎛ ′′′ ′′′ ′′′′ ′′′ 444342 343332 2423 141312 aaa0 aaa0 aa10 aaa1 Lấy hàng 1 trừ đi hàng 2 đã nhân với 12a′ ,lấy hàng 3 trừ đi hàng 2 đã nhân với 32a′ và lấy hàng 4 trừ đi hàng 2 đã nhân với 42a′ thì định thức vẫn là D/p1 và ma trận là: thì định thức vẫn là D/(p1p2) và ma trận là: 117 ⎟⎟ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎜⎜ ⎝ ⎛ ′′′′ ′′′′ ′′′′ ′′′′ 4443 3433 2423 1413 aa00 aa00 aa10 aa01 Tiếp tục lấy hàng 3 rồi hàng 4 làm trụ thì ma trận sẽ là: ⎟⎟ ⎟ ⎠ ⎞ ⎜⎜ ⎜ ⎝ ⎛ 1000 0100 0010 0001 Định thức của ma trận này là D/(p1p2p3p4)= D/( 44332211 aaaa ′′′′′′ ) =1 nên định thức của ma trận A là D = p1p2p3p4. Sau đây là ch−ơng trình tìm định thức của một ma trận: Ch−ơng trình 9-1 //tinh dinh thuc #include #include #include #include void main() { int i,j,k,n,ok1,ok2,t; float d,c,e,f,g,h; float a[50][50]; char tl; clrscr(); printf("** TINH DINH THUC CAP n **"); printf("\n"); printf("\n"); printf("Cho cap cua dinh thuc n = "); scanf("%d",&n); printf("Nhap ma tran a\n"); for (i=1;i<=n;i++) { printf("Dong %d:\n",i); for (j=1;j<=n;j++) { printf("a[%d][%d] = ",i,j); scanf("%f",&a[i][j]); } printf("\n"); } printf("\n"); printf("Ma tran a ma ban da nhap\n"); printf("\n"); for (i=1;i<=n;i++) { for (j=1;j<=n;j++) 118 printf("%.5f\t",a[i][j]); printf("\n"); } printf("\n"); t=1; flushall(); while (t) { printf("Co sua ma tran a 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; } printf("Ma tran a ban dau\n"); printf("\n"); for (i=1;i<=n;i++) { for (j=1;j<=n;j++) printf("%.5f\t",a[i][j]); printf("\n"); } printf("\n"); d=1; i=1; ok2=1; while ((ok2)&&(i<=n)) { if (a[i][i]==0) { ok1=1; k=k+1; while ((ok1)&&(k<=n)) if (a[k,i]!=0) { for (j=i;j<=n;j++) { c=a[i][j]; a[i][j]=a[k][j]; a[k][j]=c; } d=-d; 119 ok1=0; } else k=k+1; if (k>n) { printf("\n"); printf("** MA TRAN SUY BIEN **"); ok2=0; d=0; } } if (a[i][i]!=0) { c=a[i][i]; for (j=i+1;j<=n;j++) a[i][j]=a[i][j]/c; for (k=i+1;k<=n;k++) { c=a[k][i]; for (j=i+1;j<=n;j++) a[k][j]=a[k][j]-a[i][j]*c; } } i=i+1; } if (ok2) { for (i=1;i<=n;i++) d=d*a[i][i]; printf("\n"); printf("** GIA TRI DINH THUC D **"); printf("\n"); printf("%.3f",d); } getch(); } Đ2.Nghịch đảo ma trận Gọi A-1 là ma trận nghịch đảo của một ma trận A bậc n ta có AA-1 = E.(trong biểu thức này E là một ma trận vuông có các phần tử trên đ−ờng chéo chính bằng 1). Dạng của ma trận E,ví dụ cấp 4,là: ⎟⎟ ⎟ ⎠ ⎞ ⎜⎜ ⎜ ⎝ ⎛ = 1000 0100 0010 0001 E 120 Ph−ơng pháp loại trừ để nhận đ−ợc ma trận nghịch đảo A-1 đ−ợc thực hiện qua nhiều giai đoạn (n),mỗi một giai đoạn gồm hai b−ớc.Đối với giai đoạn thứ k: - chuẩn hoá phần tử akk bằng cách nhân hàng với nghịch đảo của nó - làm cho bằng không các phần tử phía trên và phía d−ới đ−ờng chéo cho đến cột thứ k.Khi k = n thì A(k) sẽ trở thành ma trận đơn vị và E trở thành A-1 Ví dụ: Tính ma trận nghịch đảo của ma trận ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛= 211 121 112 A Ta viết lại ma trận A và ma trận đơn vị t−ơng ứng với nó ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛=⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛= 100 010 001 E 211 121 112 A Giai đoạn 1: B−ớc a: Nhân hàng 1 với 1/a11,nghĩa là a,1j = a1j/a11 đối với dòng thứ nhất,a,ij = aij đối với các dòng khác ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ =⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ = 100 010 0021 E 211 121 21211 A B−ớc b: Trừ hàng 3 và hàng 2 cho hàng 1,nghĩa là a(1)1j = aij - ai1aij đối với i ≠ 1. ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ − −=⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ = 1021 0121 0021 E 23210 21230 21211 A Giai đoạn 2: B−ớc a: Lấy hàng 2 làm chuẩn,nhân hàng 2 với 2/3,để nguyên các hàng khác ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ − −=⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ = 1021 03231 0021 E 23210 3110 21211 A B−ớc b: Lấy hàng 1 trừ đi hàng 2 nhân 1/2 và lấy hàng 3 trừ đi hàng 2 nhân 1/2 ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ −− − − =⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ = 13131 03231 03132 E 3400 3110 3101 A Giai đoạn 3: B−ớc a: Lấy hàng 3 làm chuẩn,nhân hàng 3 với 3/4,để nguyên các hàng khác ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ −− − − =⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ = 434141 03231 03132 E 100 3110 3101 A B−ớc b: Lấy hàng 1 trừ đi hàng 3 nhân 1/3 và lấy hàng 2 trừ đi hàng 3 nhân 1/3 ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ −− −− −− =⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛= 434141 414341 414143 E 100 010 001 A Nh− vậy A-1 là: 121 ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ −− −− −− =− 434141 414341 414143 A 1 áp dụng ph−ơng pháp này chúng ta có ch−ơng trình sau: Ch−ơng trình 9-2 #include #include #include #include #include void main() { int i,j,k,n,t,t1; float c,a[50][50],b[50][50]; char tl; clrscr(); printf(" **MA TRAN NGHICH DAO** \n"); printf("Cho bac cua ma tran n = "); scanf("%d",&n); printf("Vao ma tran ban dau a\n"); for (i=1;i<=n;i++) { printf("Vao hang thu %d :\n",i); for (j=1;j<=n;j++) { printf("a[%d][%d] = ",i,j); scanf("%f",&a[i][j]); } printf("\n"); } printf("\n"); printf("Ma tran ban da nhap\n"); printf("\n"); for (i=1;i<=n;i++) { for (j=1;j<=n;j++) printf("%.5f\t",a[i][j]); printf("\n"); } t=1; flushall(); 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; } printf("\nMa tran ban dau\n"); 122 printf("\n"); for (i=1;i<=n;i++) { for (j=1;j<=n;j++) printf("%.5f\t",a[i][j]); printf("\n"); } printf("\n"); for (i=1;i<=n;i++) for (j=n+1;j<=2*n;j++) { if (j==i+n) a[i][j]=1; else a[i][j]=0; } i=1; t1=1; while (t1&&(i<=n)) { if (a[i][i]==0) { t=1; k=i+1; while (t&&(k<=n)) if (a[k][i]!=0) { for (j=1;j<=2*n;j++) { c=a[i][j]; a[i][j]=a[k][j]; a[k][j]=c; } t=0; } else k=k+1; if (k==n+1) { if (a[i][k-1]==0) { printf("MA TRAN SUY BIEN\n "); t1=0; } } } if (a[i][i]!=0) { c=a[i][i]; for (j=i;j<=2*n;j++) 123 a[i][j]=a[i][j]/c; } for (k=1;k<=n;k++) { if (k!=i) { c=a[k][i]; for (j=i;j<=2*n;j++) a[k][j]=a[k][j]-a[i][j]*c; } } i=i+1; } if (t1) { printf("\n"); printf("\nMA TRAN KET QUA\n"); printf("\n"); for (i=1;i<=n;i++) { for (j=n+1;j<=2*n;j++) printf("%.4f\t\t",a[i][j]); printf("\n"); } printf("\n"); } getch(); } Dùng ch−ơng trình tính nghịch đảo của ma trận: ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ 678 789 899 cho ta kết quả ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ −− − −− 991 9102 121 Đ3.Tích hai ma trận Giả sử ta có ma trận Amn và ma trận Bnp.Tích của Amn và Bnp là ma trận Cmp trong đó mỗi phần tử của Cmp là: ∑ = = n 1k kjikij bac Ch−ơng trình d−ới đây thực hiện nhân hai ma trận với nhau. Ch−ơng trình 9-3 #include #include #include #include #include #define max 50 void main() 124 { int n,l,m,i,j,k,t; float a[max][max],b[max][max],c[max][max]; char tl; clrscr(); printf("Cho so hang cua ma tran a : "); scanf("%d",&n); printf("Cho so cot cua ma tran a : "); scanf("%d",&l); printf("Cho so cot cua ma tran b : "); scanf("%d",&m); printf("\nNHAP MA TRAN A\n"); for (i=1;i<=n;i++) for (j=1;j<=l;j++) { printf("a[%d][%d] = ",i,j); scanf("%f",&a[i][j]); } printf("\n"); printf("Ma tran a ma ban da nhap\n"); for (i=1;i<=n;i++) { for (j=1;j<=l;j++) printf("%10.5f",a[i][j]); printf("\n"); } flushall(); t=1; while (t) { 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; } printf("Ma tran a ban dau"); printf("\n"); for (i=1;i<=n;i++) { for (j=1;j<=l;j++) 125 printf("%10.5f",a[i][j]); printf("\n"); } printf("\n"); printf("NHAP MA TRAN B\n"); for (i=1;i<=l;i++) for (j=1;j<=m;j++) { printf("b[%d][%d] = ",i,j); scanf("%f",&b[i][j]); } printf("\n"); printf("Ma tran b ban da nhap\n"); for (i=1;i<=l;i++) { for (j=1;j<=m;j++) printf("%10.5f",b[i][j]); printf("\n"); } flushall(); t=1; while (t) { 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("b[%d][%d] = ",i,j); scanf("%f",&b[i][j]); } if (toupper(tl)=='K') t=0; } printf("Ma tran b ban dau"); printf("\n"); for (i=1;i<=l;i++) { for (j=1;j<=m;j++) printf("%10.5f",b[i][j]); printf("\n"); } printf("\n"); for (i=1;i<=n;i++) for (j=1;j<=m;j++) { 126 c[i][j]=0; for (k=1;k<=l;k++) c[i][j]=c[i][j]+a[i][k]*b[k][j]; } printf("Ma tran tich c :\n"); for (i=1;i<=n;i++) { for (j=1;j<=m;j++) printf("%10.5f",c[i][j]); printf("\n"); } getch(); } Dùng ch−ơng trình tính tính hai ma trận ta nhận đ−ợc kết quả ⎟⎟ ⎟ ⎠ ⎞ ⎜⎜ ⎜ ⎝ ⎛ −− − − −=⎟⎠ ⎞⎜⎝ ⎛ − −ì⎟⎟ ⎟ ⎠ ⎞ ⎜⎜ ⎜ ⎝ ⎛− 1214 221 11148 105 343 221 35 01 31 12 Đ4.Giá trị riêng và vec tơ riêng của ma trận 1.Khái niệm chung: Trong nghiên lí thuyết và ứng dụng,ta gặp bài toán về ma trận cấp n.Cho một ma trận A cấp n,giá trị λ đ−ợc gọi là giá trị riêng và vectơ X đ−ợc gọi là vectơ riêng của ma trận A nếu: AX = λX (1) Vectơ riêng phải là vectơ khác không.T−ơng ứng với một giá trị riêng có vô số vectơ riêng.Nếu X là một véc tơ riêng t−ơng ứng với giá trị riêng λ thì cX cũng là vec t− riênh ứng với λ.Có nhiều thuật toán tìm giá trị riêng và vectơ riêng của một ma trận.Giả sử ta có ma trận A,gọi E là ma trận đơn vị thì theo (1) ta có: (A-λE)X = 0 (2) và (A - λE) là ma trận có dạng: ⎟⎟ ⎟⎟ ⎟ ⎠ ⎞ ⎜⎜ ⎜⎜ ⎜ ⎝ ⎛ − − − λ λλ aaa ....... aaa aaa nn2n1n n22221 n11211 .... .... .... (3) Nh− vậy do (2) là hệ ph−ơng trình tuyến tính thuần nhất nên điều kiện cần và đủ để λ là giá trị riêng của ma trận trên là định thức của nó bằng không: det(A - λE) = 0 (4) Ph−ơng trình (4) đ−ợc gọi là ph−ơng trình đặc tr−ng của ma trận A.Định thức det(A - λE) đ−ợc gọi là định thức đặc tr−ng của ma trận A.Định thức PA(λ) của ma trận trên đ−ợc gọi là đa thức đặc tr−ng của ma trận vuông A. Ví dụ tìm vec tơ riêng và trị riêng của ma trận: ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ − − − 022 113 313 Tr−ớc hết ta tính đa thức đặc tr−ng của ma trận A: 127 )4()4( 22 113 313 )(P 2 A + ⎟⎟ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎜⎜ ⎝ ⎛ λλ−=λ−− −λ− −λ−=λ Nghiệm của PA(λ) = 0 là λ1 = 4,λ2 = 2j và λ3 = -2j.Vì tr−ờng cơ sở là số thực nên ta chỉ lấy λ = 4.Để tìm vec tơ riêng t−ơng ứng với λ = 4 ta giải hệ 0 22 113 313 3 2 1 = ⎟⎟ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎜⎜ ⎝ ⎛ ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ ξ ξ ξ ì λ−− −λ− −λ− ta nhận đ−ợc các giá trị của ξ,chúng tạo thành vec tơ riêng ứng với λ. Nh− vậy khi khai triển định thức ta có một đa thức bậc n có dạng: Pn(λ) = λn - p1λn-1 - p2λn-2 - - pn = 0 Muốn xác định các hệ số của đa thức đặc tính này ta dùng ph−ơng pháp Fadeev-Leverrier.Ta xét ma trận A: ⎟⎟ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎜⎜ ⎝ ⎛ ⋅⋅⋅ ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ ⋅⋅⋅ ⋅⋅⋅ = nn2n1n n22221 n11211 aaa aaa aaa A Ta gọi vết của ma trận A là số: vet(A)= a11 + a22 +...+ ann Khi đó tham số pi của Pn(λ) đ−ợc các định nh− sau: p1 = vet(B1) với B1 = A p2 = (1/2)vet(B2) với B2 = A(B1-p1E) p3 = (1/3)vet(B3) với B3 = A(B2-p2E) ...... Ch−ơng trình tính các hệ số pi nh− sau: Ch−ơng trình 9-4 // Faddeev_Leverrier; #include #include #include #define max 50 void main() { int i,j,k,m,n,k1,t; float vet,c1,d; char tl; float p[max]; float a[max][max],b[max][max],c[max][max],b1[max][max]; clrscr(); printf("Cho bac cua ma tran n = "); scanf("%d",&n); printf("Cho cac phan tu cua ma tran a : \n"); 128 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; 129 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: ∑ = =+⋅⋅⋅++= n 1i iinn2211 XvXvXvXvV (5) Trong đó X1,X2,..,Xn là các vec tơ riêng t−ơng ứng với các giá trị riêng λ1,λ2,λ3,..,λn 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 = v1λ1X1 ; Av2X2 = v2AX2 = v2λ2X2 v.v. Vậy nên: AV = v1λ1X1 + v2λ2X2 +...+ vnλnXn 130 XvXAvAV ii n 1i iii n 1i i λ∑∑ == == Lại nhân biểu thức trên với A ta có: A2V = v1λ1 AX1 + v2λ2 AX2 +...+ vnλn AXn = v1λ21X1 + v2λ22 X2 +...+ vnλn2 Xn và tiếp đến lần thứ p ta có: XvXvXvXvVA n p nn2 p 211 p 11i p i n 1i i p λλλλ∑ +⋅⋅++== = Lấy λp1 làm thừa số chung ta có: ⎥⎥⎦ ⎤ ⎢⎢⎣ ⎡ ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ λ λ+⋅⋅⋅+⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ λ λ+⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ λ λ+λ= n p 1 n n3 p 1 3 32 p 1 2 211 p 1 p XvXvXvXvVA T−ơng tự ta có: ⎥⎥⎦ ⎤ ⎢⎢⎣ ⎡ ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ λ λ+⋅⋅⋅+⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ λ λ+⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ λ λ+λ= +++ ++ n 1p 1 n n3 1p 1 3 32 1p 1 2 211 1p 1 1p XvXvXvXvVA Khi p rất lớn,vì λ1 > λ2 > λ3 >,...,λn nên: ∞→→⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ λ λ pkhi0 1 i Do đó: 11 p 1 p p XvVAlim λ=∞→ (6) 11 1p 1 1p p XvVAlim ++∞→ λ= nghĩa là khi p đủ lớn thì: 11 p 1 p XvVA λ= 11 1p 1 1p XvVA ++ λ= do đó: VAVA p1 1p λ=+ hay: ( ) VAVAA p1p λ= Nh− vậy VA p là véc tơ riêng của A ứng với λ1 còn giá trị riêng λ1 sẽ là: 1p 1p p VA VA lim λ= + ∞→ Trong thực tế để tránh v−ợt quá dung l−ợng bộ nhớ khi λ1 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 p1p VAV ′=+ với vp+1,j là phần tử thứ j của Vp+1.Ta có: ⎪⎩ ⎪⎨ ⎧ λ= =′ +∞→ ∞→ 1j,1p p 1p p vlim XVlim 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: 131 ⎟⎟ ⎟ ⎠ ⎞ ⎜⎜ ⎜ ⎝ ⎛ −−−− = 26544323 68102 720138 17302417 A 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 λ 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 λ 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 λ 6.9634 6.9742 Dùng thuật toán trên ta có ch−ơng trình sau: Ch−ơng trình 9-5 #include #include #include #include #include #define max 50void 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++) 132 { 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++) 133 { 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: ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ −− − 308 649 012 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à λ2 > λ3 >...> λn.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: 134 Giả sử X1 là vec tơ riêng của ma trận A t−ơng ứng với giá trị riêng λ1 và W1 là vec tơ riêng của ma trận AT t−ơng ứng với giá trị riêng λ1.Từ định nghĩa AX1 = λ1X1 ta viết: (A - λE)X1 = 0 Ta tạo ma trận A1 dạng: WX XW AA T 11 1 T 1 1 1 λ−= (7) Ta chú ý là X1W1 T là một ma trận còn W1 TX1 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ó: 0 XAX XW XWXAX XWX XW AXXA 111 1 T 1 1 T 1 111 1 T 11 1 T 1 1 111 = −= −= −= λ λ λ (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 λ2,thì khi nhân A1 với X2 ta có: XW XWXAX XWX XW AXXA 1 T 1 2 T 1 112 2 T 11 1 T 1 1 221 λ λ −= −= (9) Theo định nghĩa vì W1 là vectơ riêng của A T nên: λ1W1 =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 = λ1WT1 Hay: W1 TA = λ1W1T (11) Khi nhân (11) với X2 ta có: λ1W1TX2 = W1TAX2 và do định nghĩa: AX2 = λ2X2 nên: λ1W1TX2 = W1T λ2X2 vậy thì: (λ1 - λ2) W1TX2 = 0 khi λ1 ≠ λ2 thì: W1 TX2 = 0 (12) Cuối cùng thay (12) vào (9) ta có: A1X2 = AX2 = λ2X2 Nh− vậy λ2 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ó λ1 và X1 ta tìm W1 là vec tơ riêng của AT ứng với giá trị riêng λ1 (ví dụ tìm W1 bằng cách giải ph−ơng trình (AT -λ1E)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 công suất 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: 135 ⎟⎟ ⎟ ⎠ ⎞ ⎜⎜ ⎜ ⎝ ⎛ −−−− = 26544323 68102 720138 17302417 A Ta đã tìm đ−ợc giá trị riêng lớn nhất λ1 = 7 và một vectơ riêng t−ơng ứng: X1 = { 1,1,0,-2} T. Ma trận AT có dạng: ⎟⎟ ⎟ ⎠ ⎞ ⎜⎜ ⎜ ⎝ ⎛ −− −−= 266717 5482030 43101324 232817 A T và theo ph−ơng trình AT -λ1E)W1 = 0 ta tìm đ−ợc vectơ W1 = {293,695,746,434}T Ta lập ma trận mới A1 theo (7): ⎟⎟ ⎟ ⎠ ⎞ ⎜⎜ ⎜ ⎝ ⎛ −−−− =λ 86814921390586 0000 434746695293 434746695293 120 7 XW WX 1 T 1 T 11 1 và: ⎟⎟ ⎟ ⎠ ⎞ ⎜⎜ ⎜ ⎝ ⎛ −−−− −−−−= 6333.240333.330833.381833.11 68102 3167.185167.235417.270917.9 3167.85167.135417.160917.0 A1 Từ ma trận A1 ta tìm tiếp đ−ợc λ2 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 9-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++) 136 { 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++) 137 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(); 138 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]; } } 139 Do (6) đúng với mọi n nên cho n = 1 , 2 , 3 , . . . ta có : | x2 - x1 | ≤ q | x1 - xo | | x3 - x2 | ≤ q | x2 - x1 | . . . . . . . . . . . . . . . . . . . | xn+1 - xn | ≤ q | xn - xn-1 | Điều này có nghĩa là dãy xi+1 - xi , một cách gần đúng,là một cấp số nhân . Ta cũng coi rằng dãy xn - y với y là nghiệm đúng của (1) , gần đúng nh- một cấp số nhân có công sai q . Nh- vậy : 1q yx yx n 1n <=− −+ (7) hay : )yx(qyx n1n −=−+ (8) T-ơng tự ta có : )yx(qyx 1n2n −=− ++ (9) Từ (8) và (9) ta có : n1n 1n2n xx xx q − −= + ++ (10) Thay giá trị của q vừa tính ở (10) vào biểu thức của q ở trên ta có : ( ) 1n1nn 2 1nn n xx2x xx xy ++ + +− −−= (11) Công thức (11) đ-ợc gọi là công thức ngoại suy Adam.Nh- vậy theo (11) tr-ớc hết ta dùng ph-ơng pháp lặp để tính giá trị gần đúng xn+2,xn+1,xn của nghiệm và sau đó theo (11) tìm đ-ợc nghiệm với sai số nhỏ hơn. Để làm ví dụ chúng ta xét ph-ơng trình : lnx - x2 + 3 = 0 Ta đ-a về dạng lặp : 3)xln(x += 3xlnx2 1 )x(f +=′ Phép lặp hội tụ trong đoạn [0.3,∝].Ta cho x1 = 1 thì tính đ-ợc : x2 = 1,7320508076 x3 = 1.883960229 x4 = 1.90614167 y = 1.909934347 Để giảm sai số ta có thể lặp nhiều lần Ch−ơng trình 8-9 //phuong phap Aitken #include #include #include #define m 5 void main() { 140 float x[m]; float epsi,n,y; int i,z; float f(float); clrscr(); printf("Cho tri so ban dau x[1] = "); scanf("%f",&x[1]); printf("Cho tri so sai so epsilon = "); scanf("%f",&epsi); printf("\n"); printf( "Ngoai suy Aitken cua ham\n"); z=0; while (z<=20) { for (i=2;i<=4;i++) x[i]=f(x[i-1]); n=x[4]-2*x[3]+x[2]; if ((fabs(n)<1e-09)||(fabs(x[1]-x[2])<epsi*fabs(x[1]))) z=20; else { y=x[2]-(x[3]-x[2])*(x[3]-x[2])/n; if (z>20) printf("Khong hoi tu sau hai muoi lan lap\n"); x[1]=y; } z=z+1; } printf("Nghiem cua phuong trinh y = %.6f",y); getch(); } float f(float x) { float s=sqrt(log(x)+3); return(s); } Với giá trị ban đầu là 1 và sai số là 1e-8,ch-ơng trình cho kết quả y = 1.9096975944 Đ10.Ph−ơng pháp Bairstow Nguyên tắc của ph-ơng pháp Bairstow là trích từ đa thức Pn(x) một tam thức Q2(x) = x2 - sx + p mà ta có thể tính nghiệm thực hay nghiệm phức của nó một cách đơn giản bằng các ph-ơng pháp đã biết. Việc chia đa thức Pn(x) cho tam thức Q2(x) đ-a tới kết quả : Pn(x) = Q2(x).Pn-2(x) + R1(x) với Pn(x) = aoxn + a1xn-1 + a2xn-2 +...+ an 141 Q2(x) = x2 - sx + p Pn-2(x) = boxn-2 + b1xn-3 + b2xn-4 +...+ bn-2 R1(x) = αx + β Để có đ-ợc một th-ơng đúng,cần tìm các giá trị của s và p sao cho R1(x) = 0 (nghĩa là α và β triệt tiêu).Với s và p đã cho,các hệ số b của đa thức Pn-2(x) và các hệ số α và β đ-ợc tính bằng ph-ơng pháp truy hồi.Các công thức nhận đ-ợc khi khai triển biểu thức Pn(x) = Q2(x).Pn-2(x) + R1(x) và sắp xếp lại các số hạng cùng bậc : aoxn + a1xn-1 + a2xn-2 +...+ an = (x2 - sx + p)( boxn-2 + b1xn-3 + b2xn-4 +...+ bn-2) Số hạng bậc Hệ số của Pn(x) Hệ số của Q2(x).Pn-2(x) xn ao bo xn-1 a1 b1 - sbo xn-2 a2 b2 - sb1 + pbo ...... ...... ..... xn-k ak bk - sbk-1 + pbk-2 x an-1 α - sbn-2 + pbn-3 xo an β + pbn-2 Nh- vậy : bo = a o (1) b1 = a1 + sbo b2 = a2 + sb1 - pbo .................. bk = ak + sbk-1 - pbk-2 α = an-1 + sbn-2 - pbn-3 β = an - pbn-2 Chúng ta nhận thấy rằng α đ-ợc tính toán xuất phát từ cùng một công thức truy hồi nh- các hệ số bk và t-ơng ứng với hệ số bn-1 bn-1 = an-1 + sbn-2 - pbn-3 = α Hệ số bn là : bn = an + sbn-1 - pbn-2 = sbn-1 + β và cuối cùng : R1(x) = αx + β = bn-1(x - s) + bn Ngoài ra các hệ số bi phụ thuộc vào s và p và bây giờ chúng ta cần phải tìm các giá trị đặc biệt s* và p* để cho bn-1 và bn triệt tiêu.Khi đó r1(x) = 0 và nghiệm của tam thức x2 - s*x + p*x sẽ là nghiệm của đa thức Pn(x).Ta biết rằng bn-1 và bn là hàm của s và p : bn-1 = f(s,p) bn = g(s,p) Việc tìm s* và p* đ-a đến việc giải hệ ph-ơng trình phi tuyến: ⎩⎨ ⎧ ==0)p,s(g 0)p,s(f Ph-ơng trình này có thể giải dễ dàng nhờ ph-ơng pháp Newton.Thật vậy với một ph-ơng trình phi tuyến ta có công thức lặp : xi+1 = xi - f(xi)/f'(xi) hay f'(xi)(xi+1 - xi) = -f(xi) Với một hệ có hai ph-ơng trình,công thức lặp trở thành: 142 J(Xi)(Xi+1 - Xi) = -F(Xi) với Xi = { si,pi}T và Xi+1 = { si+1,pi+1}T )p,s(g )p,s(f )X(F ii ii i = iJ X f s f p g s g p ( )= ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ Quan hệ : J(Xi)∆X = -F(Xi) với ∆X = {si+1 - si,pi+1 - pi}T t-ơng ứng với một hệ ph-ơng trình tuyến tính hai ẩn số ∆s = si+1 - si và ∆p = pi+1 - pi : ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ f s s f p p f s p g s s g p p g s p i i i i ∆ ∆ ∆ ∆ + = + = ⎧ ⎨ ⎪⎪⎪ ⎩ ⎪⎪⎪ − − ( ) ( ) , , Theo công thức Cramer ta có : ∆s f g p g f p= − +∂∂ ∂ ∂ δ ∆p g f s f g s= − +∂∂ ∂ ∂ δ δ ∂∂ ∂ ∂ ∂ ∂ ∂ ∂= − f s g p f p g s Để dùng đ-ợc công thức này ta cần tính đ-ợc các đạo hàm ∂∂ f s , ∂∂ f p , ∂∂ g s , ∂∂ g p .Các đạo hàm này đ-ợc tính theo công thức truy hồi. Do bo = ao nên ob s ∂ ∂ =0 ob p ∂ ∂ =0 b1 = a1 + sbo nên 1∂∂ b s bo = 1 0∂∂ b p = b2 = a2 + sb1- pbo nên 2 2 1∂∂ ∂ ∂ ∂ ∂ ∂ ∂ b s a s sb s pb s o= + − ( ) ( ) Mặt khác : 2 0∂∂ a s = 1 1 1 ∂ ∂ ∂ ∂ ( )sb s s b s b= + opb s ∂ ∂ ( ) =0 nên : 2 1 ∂ ∂ b s b sbo = + b3 = a3 + sb2- pb1 nên 3 2 2 1∂∂ ∂ ∂ ∂ ∂ b s b s b s p b s = + − Nếu chúng ta đặt : k k b s c ∂ ∂ = −1 143 thì : co = bo (2) c1 = b1 + sbo = b1 + sco c2 = b2 + sc1 - pco .................... ck = bk + sck-1 - pck-2 cn-1 = bn-1 + scn-2 - pcn-3 Nh- vậy các hệ số cũng đ-ợc tính theo cách nh- các hệ số bk.Cuối cùng với f = bn-1 và g = bn ta đ-ợc: 2n1n3n2n cs f c s f c s f c s f −−−− =∂ ∂=∂ ∂=∂ ∂=∂ ∂ 2 2n3n1n 3nn2n1n ccc cbcb s −−− −−− − −=∆ (3) 2 2n3n1n 2nn1n1n ccc cbcb p −−− −−− − −=∆ (4) Sau khi phân tích xong Pn(x) ta tiếp tục phân tích Pn-2(x) theo ph-ơng pháp trên Các b-ớc tính toán gồm : - Chọn các giá trị ban đầu bất kì s0 và p0 - Tính các giá trị bo,..,bn theo (1) - Tính các giá trị co,...,cn theo (2) - Tính ∆so và ∆po theo (3) và (4) - Tính s1 = s0 + ∆so và p1 = po+ ∆po - Lặp lại b-ớc 1 cho đến khi pi+1 = pi = p và si+1 = si = s - Giải ph-ơng trình x2 - sx + p để tìm 2 nghiệm của đa thức - Bắt đầu quá trình trên cho đa thức Pn-2(x) Ví dụ : Tìm nghiệm của đa thức P4(x) = x4 - 1.1x3 + 2.3x2 + 0.5x2 + 3.3. Với lần lặp ban đầu ta chọn s = -1 và p =1,nghĩa là tam thức có dạng x2 + x + 1 a0 a1 a2 a3 a4 1 -1.1 2.3 0.5 3.3 sbi -1 2.1 -3.4 0.8 -pbi-1 -1 2.1 -3.4 bi 1 -2.1 3.4 -0.8 = bn-1 0.7=bn sbi -1.0 3.1 -5.5 -pbi-1 -1.0 3.1 ci 1 -3.1 5.5 -3.2 11.0 5.52.3 1.35.5 5.57.0 1.38.0 s = − − − − =∆ 06.0 5.52.3 1.35.5 7.02.3 8.05.5 p = − − −=∆ 144 s* = -1 + 0.11 = -0.89 p* = 1 + 0.06 = 1.06 Tiếp tục lặp lần 2 với s1 = s* và p1 = p* ta có : a0 a1 a2 a3 a4 1 -1.1 2.3 0.5 3.3 sbi -0.89 1.77 -2.68 0.06 -pbi-1 -1.06 2.11 -3.17 bi 1 -1.99 3.01 -0.07 = bn-1 0.17=bn sbi -0.89 2.56 -4.01 -pbi-1 -1.0 3.1 ci 1 -2.88 4.51 -1.03 01.0 51.403.1 88.251.4 5.57.0 88.207.0 s −= − − − − =∆ 04.0 51.403.1 88.251.4 17.003.1 07.051.4 p = − − −−=∆ s* = -0.89 - 0.01 = -0.9 p* = 1.06 + 0.04 = 1.1 Nh- vậy P4(x) = (x2+0.9x+1.1)(x2 + 2x+3) Ch-ơng trình sau áp dụng lí thuyết vừa nêu để tìm nghiệm của đa thức. Ch−ơng trình 8-10 //phuong phap Bairstow #include #include #include #include #define m 10 void main() { float a[m],b[m],c[m]; int i,n,v; float s,e1,t,p,q,r,p1,q1; clrscr(); printf("Cho bac cua da thuc n = "); scanf("%d",&n); printf("Cho cac he so cua da thuc can tim nghiem\n"); for (i=n;i>=0;i--) { 145 printf("a[%d] = ",n-i); scanf("%f",&a[i]); } printf("\n"); e1=0.0001; if (n<=2) if (n==1) { printf("Nghiem cua he\n"); printf("%.8f",(a[0]/(-a[1]))); getch(); exit(1); } do { v=0; p=1; q=-1; b[n]=a[n]; c[n]=a[n]; do { b[n-1]=b[n]*p+a[n-1]; c[n-1]=b[n-1]+b[n]*p; for (i=n-2;i>=0;i--) { b[i]=b[i+2]*q+b[i+1]*p+a[i]; c[i]=c[i+2]*q+c[i+1]*p+b[i]; } r=c[2]*c[2]-c[1]*c[3]; p1=p-(b[1]*c[2]-b[0]*c[3])/r; q1=q-(b[0]*c[2]-b[1]*c[1])/r; if ((fabs(b[0])<e1)&&(fabs(b[1])<e1)) goto tt; v=v+1; p=p1; q=q1; } while (v<=40); if(v>40) { printf("Khong hoi tu sau 40 lan lap"); getch(); exit(1); } tt:s=p1/2; t=p1*p1+4*q1; if(t<0) { printf("Nghiem phuc\n"); 146 printf("%.8f+%.8fj\n",s,(sqrt(-t)/2)); printf("%.8f-%.8fj\n",s,(sqrt(-t)/2)); printf("\n"); } else { printf("Nghiem thuc\n"); printf("%.8f\n",(s+sqrt(t)/2)); printf("%.8f\n",(s-sqrt(t)/2)); printf("\n"); } for (i=2;i<=n;i++) a[i-2]=b[i]; n=n-2; } while ((n>2)&(r!=0.0)); s=-a[1]/(2*a[2]); t=a[1]*a[1]-4*a[2]*a[0]; if (t<0) { printf("Nghiem phuc\n"); printf("%.8f+%.8fj\n",s,(sqrt(-t)/(2*a[2]))); printf("%.8f-%.8fj\n",s,(sqrt(-t)/(2*a[2]))); printf("\n"); } else { printf("Nghiem thuc\n"); printf("%.8f\n",(s-sqrt(t)/(2*a[2]))); printf("%.8f\n",(s-sqrt(t)/(2*a[2]))); printf("\n"); } getch(); } Dùng ch-ơng trình trên để xác định nghiệm của đa thức : x6 - 2x5 - 4x4 + 13x3 - 24x2 + 18x - 4 = 0 ta nhận đ-ợc các nghiệm : x1 = 2.61903399 x2 = -2.73205081 x3 = 0.732050755 x4 = 0.381966055 x5 = 0.500011056 + i*1.3228881 x6 = 0.500011056 - i*1.3228881 Đ11.Hệ ph−ơng trình phi tuyến Ph-ơng pháp Newton có thể đ-ợc tổng quát hoá để giải hệ ph-ơng trình phi tuyến dạng : 147 ⎪⎪ ⎪ ⎩ ⎪⎪ ⎪ ⎨ ⎧ = = = = 0)x,...,x,x,x(f ............................... 0)x,...,x,x,x(f 0)x,...,x,x,x(f 0)x,...,x,x,x(f 2321n 23213 23212 23211 hay viết gọn hơn d-ới dạng : F(X) = 0 Trong đó : X = (x1,x2,x3,.....,xn) Với một ph-ơng trình một biến,công thức Newton là : i i i i x x f x f x+ = −1 ( ) '( ) hay : f'(xi).∆x = -f(xi) với ∆x = xi+1 - xi Đối với hệ,công thức lặp là : J(Xi)∆x = -F(Xi) Trong đó J(Xi) là toán tử Jacobi.Nó là một ma trận bậc n ( n - t-ơng ứng với số thành phần trong vectơ X) có dạng : i n n n n n n n J x f x f x f x f x f x f x f x f x f x f x f x f x ( ) . . .............. .............. .............. = 1 1 1 2 1 3 1 2 1 2 2 2 3 2 1 2 3 ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ và ∆X = Xi+1 - Xi Ph-ơng pháp Newton tuyến tính hoá hệ và nh- vậy với mỗi b-ớc lặp cần giải một hệ ph-ơng trình tuyến tính (mà biến là ∆xi) xác định bởi công thức lặp cho tới khi vectơ X(x1,x2,x3,.....,xn) gần với nghiệm. D-ới đây là ch-ơng trình giải hệ ph-ơng trình phi tuyến 1 3 2 3 1 2 4 1 2 3 4 1 2 3 1 2 3 4 3 8 0 5 0 25 8 4 0 2 8 0 x x x x x x x x x x x x x x x − − = + + + − = − + + = − + = −⎧ ⎨ ⎪⎪⎪ ⎩ ⎪⎪⎪ Ma trận đạo hàm riêng J(xi)là : 1 2 2 4 2 2 1 4 1 2 1 2 3 2 3 2 3 3 3 3 3 0 3 1 1 1 1 25 0 8 0 2 2 2 1 1 2 x x x x x x x x x x x x x x x x − − − − − − −⎛ ⎝ ⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜ ⎞ ⎠ ⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟ Ma trận này đ-ợc ch-ơng trình đọc vào nhờ thủ tục doc.Trong thủ tục này,các hệ số a[i,5] là các hàm fi(x).Vectơ nghiệm ban đầu đ-ợc chọn là { 0,-1,-1,1}T.Kết quả tính cho ta : x = {0.01328676,-1.94647929,-1.12499779,8.05819031 }T với độ chính xác 0.000001.Vectơ số d- r = { 0.00000536,-0.00000011,-0.00000001,-0.00000006}T. 148 Ch−ơng trình 8-11 //giai he pt phi tuyen #include #include #include #include #define n 4 float a[n+1][n+2]; float x[n+1],y[n+1]; int i,j,k,l,z,r; float e,s,t; void main() { void doc(); clrscr(); printf("Cho cac gia tri nghiem ban dau\n"); for (i=1;i<=n;i++) { printf("x[%d] = ",i); scanf("%f",&x[i]); } e=1e-6; z=30; for (r=1;r<=z;r++) { doc(); for (k=1;k<=n-1;k++) { s=0 ; for (i=k;i<=n;i++) { t=fabs(a[i][k]); if (s<=t) { s=t; l=i; } } for (j=k;j<=n+1;j++) { s=a[k][j]; a[k][j]=a[l][j]; a[l][j]=s; } if (a[1][1]==0) { 149 printf("Cac phan tu duong cheo cua ma tran bang khong"); getch(); exit(1); } else { if (fabs(a[k][k]/a[1][1])<(1e-08)) { printf("Ma tran suy bien"); goto mot; } } for (i=k+1;i<=n;i++) { if (a[k][k]==0) { printf("Cac phan tu duong cheo cua ma tran bang khong\n"); goto mot; } s=a[i][k]/a[k][k]; a[i][k]=0; for (j=k+1;j<=n+1;j++) a[i][j]=a[i][j]-s*a[k][j]; } y[n]=a[n][n+1]/a[n][n]; for (i=n-1;i>=1;i--) { s=a[i][n+1]; for (j=i+1;j<=n;j++) s=s-a[i][j]*y[j]; if (a[i][i]==0) { printf("Cac phan tu duong cheo cua ma tran bang khong\n"); goto mot; } y[i]=s/a[i][i]; } } if (r!=1) for (i=1;i<=n;i++) { if (fabs(y[i])<e*fabs(x[i])) goto ba; } for (i=1;i<=n;i++) x[i]=x[i]-y[i]; printf("\n"); } 150 printf("Khong hoi tu sau %d lan lap\n",z); goto mot; clrscr(); ba:printf("Vec to nghiem\n"); for (i=1;i<=n;i++) printf("%.5f\n",(x[i]-y[i])); printf("\n"); printf("Do chinh xac cua nghiem la %.5f: \n", e); printf("\n"); printf("Vec to tri so du :\n"); for (i=1;i<=n;i++) printf("%.5f\n",(a[i][n+1])); mot:printf("\n"); getch(); } void doc() { a[1][1]=3*x[1]*x[1]-3*x[2]*x[4]; a[1][2]=-3*x[2]*x[2]-3*x[1]*x[4]; a[1][3]=0; a[1][4]=-3*x[1]*x[2]; a[1][5]=x[1]*x[1]*x[1]-x[2]*x[2]*x[2]-3*x[1]*x[2]*x[4]-8; a[2][1]=1; a[2][2]=1; a[2][3]=1; a[2][4]=1; a[2][5]=x[1]+x[2]+x[3]+x[4]-5; a[3][1]=-x[1]/sqrt(25-x[1]*x[1]); a[3][2]=0; a[3][3]=8; a[3][4]=0; a[3][5]=sqrt(25-x[1]*x[1])+8*x[3]+4; a[4][1]=2*x[2]*x[3]; a[4][2]=2*x[1]*x[3]; a[4][3]=2*x[1]*x[2]; a[4][4]=-1; a[4][5]=2*x[1]*x[2]*x[3]-x[4]+8; }