Bài toán trị riêng trong phương pháp phần tử hữu hạn (PTHH) giải cho hệ dầm liên tục

tụ của lặp ngược và Lặp xuôi vectơ là trong lặp ngược λp xuất hiện ở mẫu số còn trong Lặp xuôi nó xuất hiện ở tử số. Điều này đã hạn chế tốc độ hội tụ của Lặp xuôi vectơ và bằng việc sử dụng trượt trị riêng có thể nhận được cặp giá trị (λ1,Φ1) hoặc (λn,Φn). Để đạt được tốc độ hội tụ lớn nhất tới và, ta cần chọn μ=(λ1+λn-1)/2 và μ=(λ2+λn)/2 từ đó có được tốc độ hội tụ tương ứng. Có thể đạt được tốc độ hội tụ cao hơn rất nhiều với trượt trị riêng trong lặp ngược vectơ so với Lặp xuôi vectơ. Vì lý do này và bởi trượt trị riêng có thể được chọn để hội tụ tới bất kỳ cặp giá trị nào, nên trong thực hành lặp ngược vectơ có vai trò quan trọng hơn rất nhiều so với lặp xuôi vectơ. μλ μλ − −= j p jp r max # μλ μλ − − −++ j p mjjjp max 1,...,1,# Trang 59 3.7. Phương pháp biến đổi ma trận hay chéo hóa ma trận Phương pháp này dựa trên tính chất trực giao của vectơ riêng qua ma trận K và M: Λ=ΦΦ KT IMT =ΦΦ Các bước cơ bản để giảm K và M tới dạng chéo là sử dụng biện pháp nhân trái và nhân phải với ma trận PkT và Pk trong đó k=1,2,...Đặc biệt, nếu ta định nghĩa K1=K và M1=M 1112 PKPK T= 2223 PKPK T= . . kk T kk PKPK =+1 . . Tương tự 1112 PMPM T= 2223 PMPM T= . . kk T kk PMPM =+1 . . Trang 60 Trong đó ma trận Pk được chọn để mang Mk và Kk gần tới dạng chéo. Ta cần có: Kk+1 -> Λ và Mk+1 -> I khi k->∞ Trong trường hợp l là bước lặp cuối cùng Φ=P1P2....Pl trong thực hành Mk+1 không cần thiết phải hội tụ tới I và Kk+1 không cần thiết phải hội tụ tới Λ nhưng chúng chỉ cần hội tụ tới dạng chéo Kk+1 -> diag(Kr) và Mk+1 -> diag(Mr) khi k->∞ vậy thì, với l là bước lặp cuối cùng Sử dụng các ý kiến cơ bản mô tả ở trên, người ta đã đề nghị nhiều phương pháp lặp khác nhau. Trong đó 2 phương pháp hiệu quả nhất trong phân tích PTHH là Jacobi và Householder-QR. Ở phần trước ta thấy phương pháp lặp được bắt đầu với nhân trái và phải bởi P1T và P1 đó thực ra là trường hợp trong phương pháp Jacobi. Tuy nhiên, mục đích đầu tiên của ta là chuyển bài toán trị riêng KΦ=λMΦ sang dạng kinh tế hơn để sử dụng phương pháp lặp. Đặc biệt, khi M=I, m bước chuyển đầu tiên thực hiện với ma trận K ở trên có thể dùng để đưa K sang dạng tam giác chéo mà không dùng phương pháp lặp, sau đó, ma trận Pi, i=m+1, ..., l được dùng trong phương pháp lặp để đưa Km+1 sang dạng chéo. Trong trường hợp này, )( )1( )1( + + =Λ l r l r M Kdiag )(.. )1( )1( 321 + + =Φ l r l r l M KdiagPPPP Trang 61 M ma trận đầu tiên P1, P2, .. ,Pm có thể có dạng khác so với các ma trận được sử dụng sau này Pm+1, Pm+2, .. ,Pl . Một trong các phương pháp để thực hiện công việc này là “phương pháp Householder-QR”, trong đó, ma trận Householder được dùng để chuyển K sang dạng tam giác chéo ở lần chuyển đầu tiên và sau đó ma trận xoay được dùng trong phép chuyển QR để đưa K sang dạng chéo. Thuật giải tương tự có thể được dùng để giải bài toán trị riêng tổng quát KΦ=λMΦ, M # I. 3.7.1 Phương pháp xoay Jacobi dùng cho bài toán chuẩn Phương pháp Jacobi được phát triển cho lời giải của bài toán trị riêng chuẩn (M = I). Ưu điểm lớn nhất của nó là đơn giản và ổn định. Khi tính chất của vectơ riêng Λ=ΦΦ KT IMT =ΦΦ được áp dụng cho tất cả các ma trận đối xứng K mà không có bất kỳ hạn chế gì về trị riêng, phương pháp Jacobi có thể được dùng để tính toán trị riêng âm, dương hoặc bằng 0. Với bài toán chuẩn, bước lặp thứ k định nghĩa ở trên giảm xuống thành Kk+1 = PTkKkPk Trong đó Pk là ma trận trực giao, nghĩa là: PTkPk = I Trong lời giải Jacobi ma trận Pk là ma trận xoay được lựa chọn theo cách 1 phần tử ngoài đường chéo trong Kk bằng 0. Nếu phần tử (i,j) giảm xuống 0, ma trận trực giao Pk tương ứng là: Trang 62 Trong đó θ được lấy từ điều kiện để phần tử (i,j) trong Kk bằng 0. Ký hiệu phần tử (i,j) trong Kk bằng kij(k) Với kii(k) # kjj(k) Và θ= pi()/4 Với kii(k) = kjj(k) Một điểm quan trọng phải nhấn mạnh là mặc dù việc biến đổi trên làm các phần tử ngoài đường chéo thứ (i,j) ở Kk bằng 0 tại bước thứ k, nhưng những phần tử này sẽ lại khác không trong quá trình chuyển hoá tiếp sau. Do đó, để thiết lập 1 thuật giả đúng, ta phải quyết định phần tử nào giảm tới 0. Một lựa chọn là cho phần tử ngoài đường chéo lớn nhất giảm tới 0. tuy nhiên, việc tìm kiếm phần tử lớn nhất tốn thời gian và người ta thường sử dụng chuyển Jacobi đối xứng, hàng với hàng hoặc cột với cột, quá trình này được gọi là Jacobi vòng. Quá trình chạy 1 lần qua tất cả các phần tử ngoài đường chéo gọi là 1 vòng quét. Nhược điểm của quá trình này là không để ý tới giá trị của phần tử ngoài đường chéo chính, nghĩa là phần tử đã bằng 0 hoặc gần bằng 0 và ma trận xoay vẫn được chấp nhận. )()( )(2 2tan k jj k ii k ij kk k −=θ Trang 63 Một quá trình được sử dụng hiệu quả là phương pháp xoay Jacobi có ngưỡng quét, trong phương pháp này, người ta định giá trị kiểm tra độ lớn phần tử ngoài đường chéo và gọi là ngưỡng quét, các phần tử ngoài đường chéo được kiểm tra liên tiếp về độ lớn, từ hàng đến hàng (hoặc từ cột đến cột), và phép quay được chấp nhận khi phần tử lớn hơn so với ngưỡng đã định. Ta có thể tóm tắt những bước được sử dụng trong lặp ngưỡng Jacobi. + Bước 1: Thiết lặp ngưỡng quét. Thông thường, ngưỡng quét sử dụng có thể là 10-2m. + Bước 2: Với mọi i, j với i<j tính tỷ số ((kij(k))2/kii(k)kjj(k))1/2 và chấp nhận xoay nếu phần tử này lớn hơn ngưỡng quét hiện thời. + Bước 3: sử dụng biểu thức sau để kiểm tra sự hội tụ. Nếu quan hệ trên không thoả mãn, tiếp tục quét (nghĩa là trở lại bước 1). Nếu quan hệ trên thoả mãn, kiểm tra điều kiện Với mọi i và j thỏa mãn i < j, có thoả mãn hay không, nếu có thì phép lặp hội tụ, nếu không, tiếp tục quét. 3.7.2 Phương pháp Jacobi dùng cho bài toán tổng quát Ở trên, ta đã sử dụng ma trận quay Jacobi truyền thống để giảm ma trận K về dạng đường chéo trong khi giải bài toán trị riêng chuẩn. Để giải bài toán trị riêng tổng quát, mà vẫn sử dụng phương pháp Jacobi Trang 64 chuẩn ta phải chuyển bài toán tổng quát sang bài toán chuẩn. Tuy nhiên, quá trình chuyển này có thể được bỏ qua nếu sử dụng phương pháp Jacobi tổng quát được thao tác trực tiếp trên ma trận K và M. Tương tự như phần trên, trong phương pháp lặp Jacobi, ta sử dụng ma trận Pk sau: Trong đó, hằng số α và γ được chọn theo cách sao cho phần tử (i,j) của Kk và Mk đồng thời giảm tới 0. Do đó, giá trị của α và γ là hàm của các yếu tố k(k)ij, k(k)ii, k(k)jj, m(k)ij, m(k)ii, m(k)ij trong đó chỉ số trên (k) là ký hiệu bước lặp thứ k. Tiến hành nhân PTkKkPk và PTkMkPk và sử dụng điều kiện để k(k+1)jj, m(k+1)jj,bằng 0, ta nhận được 2 phương trình sau cho α và γ. 0)1( )()()( =+++ kjjkijkii kkk γαγα 0)1( )()()( =+++ kjjkijkii mmm γαγα Nói chung, để giải α và γ ta định nghĩa Trang 65 Và Giá trị x ở trên được xác định bởi Quan hệ cho α và γ được sử dụng và ưu tiên thực hiện trong trường hợp M là ma trận xác định dương hoặc dạng dải. Trong trường hợp đó, ta có Và từ đây x luôn khác 0. Thêm nữa, det Pk # 0, là điều kiện cần thiết để thuật toán hoạt động. Nếu M là ma trận khối lượng dạng chéo, M # I và mii # 0, trong trường hợp này ta sử dụng quan hệ sau thay cho hệ phương trình trên Nếu M=I, ta có α=-γ và ta nhận thấy rằng Pk chính là ma trận xoay trong Jacobi chuẩn. Thêm nữa, cần lưu ý là quá trình giải có thể thích hợp để giải bài toán KΦ=ΛΜΦ khi M la ma trận chéo với vài phần tử trên đường chéo bằng 0. Toàn bộ quá trình giải cũng tương tự như trong lặp Jacobi của bài toán chuẩn. Sự khác nhau là tỷ số cặp khối lượng phải tính toán, trừ khi M là ma trận chéo và quá trình biến đổi được áp dụng cho Kk và Mk. Sự Trang 66 hội tụ được xác định bằng việc so sánh các trị riêng xấp xỉ liên tiếp và bởi việc kiểm tra điều kiện tất cả các phần tử ngoài đường chéo đủ nhỏ; nghĩa là với l là bước lặp cuối cùng, kết quả được chấp nhận nếu: Trong đó: Và với mọi i<j (trong đó 10-s là dung sai hội tụ) Phương pháp Jacobi là phương pháp được sử dụng rộng rãi trong lập trình để giải bài toán trị riêng có kích thước vừa và nhỏ, không có yêu cầu cao về tính chất của ma trận khối lượng và ma trận độ cứng. Phương pháp xoay Jacobi bộ phận cấu thành phương pháp lặp không gian con, dùng để giải bài toán trị riêng trong không gian con. Với ý nghĩa trên, tác giả đã thử nghiệm lập trình bằng ngôn ngữ lập trình Matlab. Cụ thể như sau: function [V,D] = xoayJACOBI(A,B); N = length(A(:,1)); X = zeros(N,N); for i = 1:N if A(i,i)>0 & B(i,i)>0, D(i) =A(i,i)/B(i,i); EIGV(i) = D(i); X(i,i)=1.; Trang 67 end end % Khoi tao vong lap NSWEEP = 0; NSMAX = 6; NR = N - 1; while NSWEEP < NSMAX NSWEEP = NSWEEP + 1; EPS = 0.01^(2*NSWEEP); for j =1:NR jj = j+1; for k = jj:N EPTOLA = (A(j,k)/A(j,j))*(A(j,k)/A(k,k)); EPTOLB = (B(j,k)/B(j,j))*(B(j,k)/B(k,k)); if (EPTOLA>EPS | EPTOLB>EPS), AKK = A(k,k)*B(j,k) - B(k,k)*A(j,k); AJJ = A(j,j)*B(j,k) - B(j,j)*A(j,k); AB = A(j,j)*B(k,k) - B(j,j)*A(k,k); SCALE = A(k,k)*B(k,k); ABCH = AB/SCALE; AKKCH = AKK/SCALE; AJJCH = AJJ/SCALE; CHECK =(ABCH*ABCH+4*AKKCH*AJJCH)/4; if CHECK > 0, SQCH = SCALE*sqrt(CHECK); Trang 68 D1 = AB/2. + SQCH; D2 = AB/2. - SQCH; DEN = D1; if (abs(D2)>abs(D1)),DEN = D2; end if (DEN == 0) CA = 0; CG = -A(j,k)/A(k,k); else CA = AKK/DEN; CG = -AJJ/DEN; end % Thuc hien xoay ma tran JP1 = j+1; JM1 = j-1; KP1 = k+1; KM1 = k-1; if JM1>0 for i =1:JM1 AJ = A(i,j); BJ = B(i,j); AK = A(i,k); BK = B(i,k); A(i,j) = AJ + CG*AK; B(i,j) = BJ + CG*BK; A(i,k) = AK + CA*AJ; B(i,k) = BK + CA*BJ; end end if (KP1<N+1) for i =KP1:N AJ = A(j,i); BJ = B(j,i); AK = A(k,i); BK = B(k,i); A(j,i) = AJ + CG*AK; B(j,i) = BJ + CG*BK; A(k,i) = AK + CA*AJ; B(k,i) = BK + CA*BJ; end end Trang 69 if (JP1<KM1+1) for i =JP1:KM1 AJ = A(j,i); BJ = B(j,i); AK = A(i,k); BK = B(i,k); A(j,i) = AJ + CG*AK; B(j,i) = BJ + CG*BK; A(i,k) = AK + CA*AJ; B(i,k) = BK + CA*BJ; end end AK = A(k,k); BK = B(k,k); A(k,k) = AK + 2.*CA*A(j,k) +CA*CA*A(j,j); B(k,k) = BK + 2.*CA*B(j,k) +CA*CA*B(j,j); A(j,j) = A(j,j) + 2.*CG*A(j,k) +CG*CG*AK; B(j,j) = B(j,j) + 2.*CG*B(j,k) +CG*CG*BK; A(j,k) = 0.0; B(j,k) = 0.0; % Tinh ma tran tri rieng for i = 1:N XJ = X(i,j); XK = X(i,k); X(i,j) = XJ + CG*XK; X(i,k) = XK + CA*XJ; end end end end end % Tinh tri rieng sau moi NSWEEP Trang 70 for i = 1:N EIGV(i) = A(i,i)/B(i,i); end % ket thuc for i = 1:N end % Lay doi xung phan duoi cac ma tran ket qua for i = 1:N for j = i:N A(j,i) = A(i,j); B(j,i) = B(i,j); end end % Tinh vec to tri rieng for j = 1:N BB = sqrt(B(j,j)); for k = 1:N X(k,j) = X(k,j)/BB; end end D = EIGV; V = X; 3.7.3 Phương pháp lặp ngược Householder – QR Phương pháp HQRI chỉ được dùng trong việc giải bài toán chuẩn. Do đó, nếu xem xét bài toán tổng quát thì trước hết phải thực hiện chuyển từ bài toán tổng quát sang bài toán chuẩn, nhưng trong thực tế việc biến đổi để chuyển từ dạng tổng quát sang dạng chuẩn chỉ có tác dụng trong một vài trường hợp. Đây là hạn chế lớn nhất của phương pháp Trang 71 này. Nhưng với nó, ta có thể giải bài toán chuẩn với K là ma trận có thể có trị riêng bằng 0 hoặc âm. Nghĩa là, không cần thiết phải sử dụng trượt trong thuật toán HQRI với mục đích giải cho những trị riêng dương. Phương pháp này dựa trên 3 bước cơ bản sau: + Biến đổi Householder để giảm ma trận K về dạng tam giác. + Lặp QR để giải ra tất cả các trị riêng. + Sử dụng lặp ngược để giải ra các trị riêng yêu cầu của ma trận tam giác ở trên. Những vectơ này sẽ được biến đổi thành vectơ riêng của ma trận K. Điểm khác biệt cơ bản với lời giải Jacobi là ma trận nhận được ở lần biến đổi đầu tiên mà không lặp là ma trận chéo ba. Ma trận này được sử dụng rất hiệu quả trong lời giải lặp QR, trong đó tất cả các vectơ riêng đều được tính toán. Cuối cùng, chỉ những vectơ riêng yêu cầu mới được tính. Nếu có nhiều vectơ riêng được yêu cầu tính, việc chuyển ma trận K sang dạng chéo ba đòi hỏi một số lượng lớn thao tác. Biến đổi Householder Với K1=K, ta tính: Pk là ma trận phản chiếu hay ma trận Householder Trang 72 Xem xét trường hợp k=1 (là trường hợp điển hình). Chia K1, P1 và w1 thành các ma trận con có hạng n-1 như sau: Trường hợp tổng quát, tại bước thứ k, ta có ma trận tương ứng có hạng n-k. Thực hiện phép nhân ta có: Ta muốn K2 có dạng chéo ba như sau: Trong đó, x là ký hiệu trị số khác không và Ta xác định 1w từ biểu thức sau Trang 73 Trong đó k21 là phần tử (2,1) của K1. Lặp QR Lặp QR được áp dụng cho ma trận chéo ba nhận được qua phép biến đổi Householder của ma trận K. Tuy nhiên, nó cũng được dùng cho ma trận gốc, nhưng để cải thiện hiệu quả của lời giải, ta biến đổi K sang dạng chéo ba trước khi lặp. Đầu tiên ta xem xét việc lặp cho ma trận K đối xứng và tổng quát. Bước cơ bản trong việc lặp là phân chia K sang dạng K=QR Trong đó: Q là ma trận có QT=Q-1 và R là ma trận dạng tam giác trên. Từ đó RQ = QTKQ Việc này có thể nhận được khi sử dụng quá trình Gram – Schmidt với các cột của ma trận K. Trong thực hành, hiệu quả hơn khi giảm ma trận K về dạng tam giác trên khi sử dụng ma trận xoay Jacobi, nghĩa là Trong đó ma trận xoay PTj,i được lựa chọn cho phần tử (j,i) bằng 0 hay nói cách khác, đó là chuyển vị của ma trận Pi,j được lựa chọn cho phần tử (i,j) bằng 0. Vậy: Thuật toán lặp QR nhận được bằng việc lặp lại quá trình ở trên. Sử dụng K1=K, ta có: Và sau đó: Trang 74 Trị riêng, vectơ riêng được tính như bình thường Tính toán các vectơ riêng Trị riêng được tính với độ chính xác cao bởi trong lặp QR có sử dụng trượt trị riêng, tốc độ hội tụ là rất nhanh. Khi trị riêng được tính với độ chính xác cao, ta có thể tính vectơ riêng yêu cầu của ma trận chéo ba T1 bằng phép lặp ngược đơn giản với trị số trượt tương ứng với trị riêng. Thông thưởng chỉ 2 bước của lặp ngược với vectơ lặp ban đầu là vectơ đơn vị là đủ. Vectơ riêng T1 cần được biến đổi Householder để nhận được vectơ riêng của K, nghĩa là nếu ký hiệu vectơ riêng thứ i của T1 là Ψi, ta có 3.8. Phương pháp lặp đa thức và phương pháp lặp với dãy Sturm 3.8.1 Lặp đa thức rõ Bước đầu tiên là viết đa thức đặc trưng p(λ) dưới dạng sau Rồi xác định các hệ số của nó. Bước thứ 2 là xác định nghiệm của đa thức. Trang 75 Khi hạng của ma trận n lớn, việc xác định các trị số của đa thức sẽ không dễ dàng, nó đòi hỏi rất nhiều thao tác. Một khó khăn khác đó là lỗi khi xác định các hệ số của đa thức. Các lỗi này là không tránh khỏi. Nhưng chỉ một lỗi nhỏ của các hệ số này cũng gây ra lỗi rất lớn trong nghiệm của đa thức. Vì thế, phương pháp này không hiểu quả trong bài toán tổng quát. 3.8.2 Lặp đa thức ẩn Trong phương pháp này ta xác định trực tiếp giá trị của p(λ) mà không tính các hệ số của nó. Giá trị của p(λ) có thể nhận được bằng cách phân tích K-λM thành ma trận tam giác đơn vị dưới L và ma trận tam giác trên S Do đó 3.8.3 Lặp dựa trên tính chất dãy Sturm Nếu ta muốn giải trị riêng nằm giữa λl và λu, trong đó λl và λu là các cận dưới và cận trên. Quá trình giải như sau: + Cấu thành K-μlM và tìm ql trị riêng nhỏ hơn λl. + Dùng tính chất dãy Sturm kiểm tra K-μuM và tìm qu trị riêng nhỏ hơn λu. Do đó, trị riêng giữa λl và λu là qu- ql. + Sử dụng sơ đồ mặt cắt kẹp để tìm các khoảng chứa tất cả các trị riêng cần tìm. + Tính trị riêng tới độ chính xác yêu cầu và sau đó nhận được vectơ riêng tương ứng bằng lặp ngược. Trang 76 3.9. Phương pháp lặp không gian con 3.9.1 Sơ bộ về phương pháp lặp không gian con Phương pháp này đặc biệt hiệu quả trong trường hợp tính toán một vài trị số trị riêng nhỏ nhất và vectơ riêng tương ứng của hệ nhiều phần tử. Phương pháp lặp không gian con được phát triển và đặt tên bởi K.J.Bathe bao gồm 3 bước chủ yếu sau: + Thiết lập q vectơ lặp ban đầu q>p (với p là số trị riêng và vectơ cần tính). + Sử dụng đồng thời lặp ngược cho q vectơ và phân tích Ritz để thu được các xấp xỉ trị riêng và vectơ riêng tốt nhất từ q vectơ lặp. + Sau khi phép lặp hội tụ, sử dụng tần số Sturm kiểm tra trị riêng yêu cầu và vectơ riêng tương ứng đã được tính. Quá trình giải đó được đặt tên là “lặp không gian con” bởi quá trình lặp là quá trình lặp với không gian con q chiều chứ không phải lặp đồng thời với q vectơ lặp riêng biệt. Đặc biệt, cần chú ý, việc chọn vectơ lặp ban đầu trong bước 1 và tần số kiểm tra Sturm trong bước 3 là rất quan trọng của quá trình lặp. Phương pháp lặp không gian con sử dụng đồng thời cả lặp vectơ, phân tích Rayleigh-Ritz, tính chất của dãy Sturm. Các ưu điểm của phương pháp lặp không gian con: lý thuyết dễ hiểu, phương pháp thực hành mạnh và dễ dàng cho việc lập trình. 3.9.2 Nội dung phương pháp lặp không gian con Mục đích chính: giải cho p trị riêng nhỏ nhất và vectơ riêng tương ứng thoả mãn KΦ = MΦΛ Trang 77 Trong đó Λ=diag (λi) và Φ = [ Φ1, Φ2, ..., Φp] Các vectơ riêng phải thoả mãn điều kiện trực giao ΦTKΦ = Λ và ΦTMΦ = I Ý tưởng cơ bản của phương pháp lặp không gian con là p vectơ riêng cần tìm tạo thành 1 không gian con gọi là ∞E . Lời giải lặp với p vectơ không phụ thuộc tuyến tính có thể được coi như lời giải lặp với không gian con. Không gian vectơ lặp ban đầu gọi là E1 và vòng lặp tiếp tục lặp cho tới khi đủ độ chính xác thì hình thành ∞E . Số bước lặp yêu cầu phụ thuộc vào việc E1 gần với ∞E ra sao chứ không phụ thuộc vào việc mỗi vectơ lặp gần với vectơ riêng ra sao. Vậy thì, hiệu quả của phương pháp này thể hiện ở việc ta dễ dàng hơn rất nhiều trong việc thiết lập không gian con ban đầu đủ gần ∞E so với tìm p vectơ mà mỗi vectơ gần với vectơ riêng yêu cầu. Thực tế chứng minh rằng, ngay cả khi ta có vectơ lặp ban đầu là tổ hợp tuyến tính của các vectơ riêng yêu cầu thì nếu sử dụng phương pháp lặp đồng thời với p vectơ cũng phải qua rất nhiều bước ta mới có thể nhận được kết quả mong muốn. Trong khi đó, nếu sử dụng phương pháp lặp không gian con, trong trường hợp này chỉ cần 1 bước ta sẽ có được kết quả. Bởi phép lặp được thực hiện trên không gian con nên sự hội tụ của không gian con là tất cả những thứ mà ta yêu cầu chứ không phải là sự hội tụ của các vectơ lặp riêng lẻ tới vectơ riêng. Sự hội tụ của phương pháp: Giả sử rằng, trong quá trình lặp, vectơ trong Xk+1 được sao cho phần tử thứ i trên đường chéo trong Λk+1 lớn hơn phần tử thứ (i-1) với i =2,..., p, vậy thì cột thứ i trong Xk+1 hội tụ tuyến tính tới Φi và tốc độ hội tụ là λi/λp+1. Mặc dù đó là tốc độ hội tụ tiệm cận, nhưng nó cũng chỉ ra rằng trị riêng nhỏ nhất hội tụ nhanh Trang 78 nhất. Thêm nữa, tốc độ lặp cao hơn có thể nhận được bằng việc sử dụng q vectơ lặp (thay vì p), trong đó q>p. Nhưng nếu sử dụng nhiều vectơ lặp thì sẽ làm tăng bộ nhớ của máy tính, nên trong thực hành, nói chung, người ta sử dụng q=min{2p,p+8}. Các bước thực hiện + Giải phương trình cân bằng tĩnh + Xây dựng hình chiếu của ma trận độ cứng K trong không gian con EK+1 k T kkq YXK .11, ++ = + Xây dựng hình chiếu của ma trận khối lượng trong không gian con EK+1 + Giải bài toán trị riêng trong không gian con EK+1 + Tính hệ véc tơ lặp cho bước lặp tiếp theo + Sau quá trình lặp trên, khi số bước lặp k→ ∞ thì Λk+1→Λ kk YXK =+1 111, 11 . +++ ++ = = k T kkq kk YXM XMY 11,1,1,1, +++++ ΛΦ=Φ kkqkqkqkq MK 1,11 +++ Φ= kqkk YY Trang 79 Yk+1→Φ 3.9.3 Một số chú ý khi chọn vectơ lặp ban đầu Nếu vectơ lặp ban đầu được sử dụng tạo thành không gian con cần tìm thì phép lặp hội tụ chỉ trong 1 bước. Đó là trường hợp ma trận khối lượng co p trị số khác 0 trên đường chéo và vectơ lặp ban đầu là vectơ đơn vị ei với các giá trị +1 tại các bậc tự do khối lượng. trong trường hợp này, phép lặp không gian con trên thực tế là phân tích quy rút tĩnh hoặc suy giảm Guyan. Trường hợp thứ 2 mà phép lặp không gian con hội tụ chỉ trong 1 bước là khi K và M đều là ma trận chéo. Đây là trường hợp khá hiếm. Khi đó, vectơ lặp nên là vectơ đơn vị với trị số +1 tương ứng ở các bậc tự do có tỷ số kii/mii nhỏ nhất. Chúng là những vectơ riêng tương ứng với trị riêng nhỏ nhất. Có thể thấy, cả 2 trường hợp trên là khá đặc biệt và trong cả 2 trường hợp, ta đều sử dụng vectơ đơn vị ei , trong đó i = r1,r2,...rp với rj (j = 1, 2, ..., p) tương ứng với p tỷ số kii/mii nhỏ nhất. Trong trường hợp tổng quát, khi việc gộp các khối lượng hoặc thuộc tính độ cứng là không thể hoặc ma trận K và M không có cùng dạng như trên thì vectơ lặp xuất phát phải được cấu tạo tương ứng với những bậc tự do có khối lượng lớn và độ cứng nhỏ, bởi những bậc tự do như thế sẽ có tần số (tức là trị riêng) nhỏ nhất. Từ đó, người ta đưa ra nguyên tắc lựa chọn vectơ lặp xuất phát: + Cột đầu tiên của tích MX1 chính là những giá trị trên đường chéo của M (Điều này đảm bảo tất cả các bậc tự do khối lượng đều được tính đến). Trang 80 + Các cột khác của tích MX1, trừ cột cuối cùng, là những vectơ đơn vị với trị số +1 tương ứng ở các bậc tự do có tỷ số kii/mii nhỏ nhất. + Cột cuối cùng của tích MX1 là vectơ tuỳ chọn. Trong phân tích hệ lớn, khoảng cách gần đúng giữa các vectơ đơn vị trong vectơ lặp xuất phát cũng rất quan trọng và phải được tính đến. Không gian xuất phát mô tả ở trên được chứng minh là khá phù hợp. Tuy nhiên, vectơ lặp xuất phát tốt hơn cũng có thể tìm được trong 1 số trường hợp cụ thể. Chẳng hạn, khi tính một số vectơ riêng khác còn lại khi biết trước 1 vài vectơ riêng hoặc trong phân tích động. Số bước lặp cần thiết phụ thuộc vào + Dạng ma trận K và M + Số lượng trị riêng cần tìm + Số lượng vectơ lặp sử dụng và độ chính xác yêu cầu. Kinh nghiệm chỉ ra rằng khi p là nhỏ, việc sử dụng vectơ lặp như mô tả ở trên và số lượng vectơ lặp q=min(2p,p+8) thì sau 10 tới 20 bước lặp ta tính được trị riêng lớn nhất với độ chính xác tới 6 số thập phân, với những trị riêng nhỏ hơn, độ chính xác còn cao hơn 3.9.4 Sự hội tụ Phương pháp lặp không gian con cũng như mọi phương pháp lặp khác trên máy là phải đánh giá được sai số trong quá trình lặp và dừng tính khi đạt độ chính xác mong muốn. Kí hiệu λi(k) và Φqi(k) trị riêng và véc tơ trong ma trận Qk tương ứng với trị riêng λi(k) tính được trong bước k-1, mức hội tụ được tính theo công thức: ( )( ) tolkqiTkqi k i ≤ ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ ΦΦ − 2/1 )()( 2)( 1 λ Trang 81 Sai số tol = 10-2s thì các trị riêng được tính với sai số đến 2s chữ số thập phân. Trị riêng càng lớn thì sai số càng lớn. Khi trị riêng cuối cùng cần tính chính xác đến 2s chữ số thập phân thì các trị riêng nhỏ hơn càng có độ chính xác cao hơn. Do xấp xỉ các trị riêng chính là thương Rayleigh nên khi các trị riêng chính xác đến 2s chữ số thập phân thì các véc tơ riêng đạt chính xác chỉ đến s chữ số thập phân (hoặc hơn). Một vấn đề quan trọng nữa trong hội tụ của phép lặp không gian con là phải xác định được đúng p trị riêng nhỏ nhất và p véc tơ riêng tương ứng. Kiến nghị xây dựng hệ véc tơ lặp xuất phát như trên sẽ dẫn đến q cặp vectơ riêng ứng với q trị riêng nhỏ nhất. Tuy nhiên điều đó mới chỉ được xác nhận bằng tính toán thử nghiệm số chứ chưa được chứng minh bằng toán học. Cách thử nghiệm số cũng tương đối đơn giản, ta có thể chọn bài toán với số bậc tự do n nhỏ để máy có thể giải lặp bằng các phương pháp lặp khác như phương pháp xoay Jacobi, khi giải theo phương pháp này ta đươc n cặp riêng, lấy p cặp ứng với p trị riêng nhỏ nhất rồi so sánh với nghiệm tìm được theo phương pháp lặp không gian con. Để chứng tỏ là ta đã tìm đúng p trị riêng nhỏ nhất, với bài toán trị riêng tổng quát mà ma trận K và ma trận khối lượng M đều xác định dương, mọi trị riêng đều hữu hạn và lớn hơn không. áp dụng tính chất dãy Sturm ta phân tích về dạng tam giác ma trận K - μM = LDLT với μ nằm ngay bên phải của λp thì trong ma trận chéo D phải chưa p giá trị âm. Trang 82 CHƯƠNG 4: THỬ NGHIỆM LẬP TRÌNH TRÊN MATLAB Trên cơ sở lý thuyết đã trình bày ở chương 2 và chương 3. Với mục tiêu hiểu biết sâu sắc hơn về các thuật toán giải các loại bài toán trị riêng chuẩn cũng như tổng quát, kích thước nhỏ cũng như kích thước lớn. Tác giả thử nghiệm lập trình trong môi trường Matlab chương trình Eigbeam. Khả năng của chương trình đạt được: 1. Phát sinh và tính toán kết cấu thành phẳng. Trong chương trình đã sử dụng phần tử thanh phẳng – tổ hợp của thanh dầm và thanh kéo nén. Vì vậy, hệ dầm liên tục là trường hợp riêng của kết cấu thành phẳng. 2. Thử nghiệm lập trình + Xây dựng phương trình trị riêng hệ dầm liên tục + Tính trị riêng nhỏ nhất và vectơ riêng bằng phương pháp lặp ngược vectơ. + Giải bài toán trị riêng kích thước vừa và nhỏ bằng thuật toán xoay Jacobi. + Giải bài toán trị riêng kích thước lớn bằng phương pháp lặp không gian con có kể đến lưu skyline của ma trận độ cứng K và ma trận khối lượng M. Trong chương này, tác giả trình bày sơ lược về tổ chức số liệu, tổ chức chương trình và một số kết quả tính minh họa. 4.1. Tổ chức số liệu trong chương trình Sơ đồ chức năng chính của chương trình bao gồm 3 khối chính Khối nhập số liệu (tên biến là SLV) Khối tính toán trung gian (tên biến là SLTG) Trang 83 Khối xuất kết quả tính (tên biến là SLR) 4.1.1 Số liệu vào SLV Việc nhập số liệu trong phương pháp phần tử hữu hạn bao gồm các bước nhập số liệu nút, phần tử, hình học, vật liệu và tải trọng − Nhóm số liệu về nhịp (SLV.Nhip) Chiều dài nhịp SLV.Nhip.KcGoi Số lượng gối SLV.Nhip.TsGoi − Nhóm số liệu về nút (SLV.Nut) Ký hiệu TsNut SLV.Nut.TsNut Mảng toạ độ nút XYZ SLV.Nut.XY Mảng bậc tự do nút JF SLV.Nut.JF − Nhóm số liệu về phần tử (SLV.Pt) Tổng số phần tử TsPhantu SLV.Pt.TsPhantu Mảng chỉ số nút LNC SLV.Pt.LNC Mảng lưu chỉ số hình học SLV.Pt.ChisoHh Mảng lưu chỉ số vật liệu SLV.Pt.ChisoVl Mảng lưu loại, kích thuớc, tiết diện SLV.Hh.HINHHOC − Nhóm số liệu về vật liệu Tổng số nhóm vật liệu SLV.Vl.TsNhomVl Mảng lưu loại, đặc trưng vật liệu SLV.Vl.VATLIEU 4.1.2 Số liệu trung gian SLTG Mảng NDF - đánh số bậc tự do NHẬP SỐ LIỆU TÍNH TOÁN XUẤT KẾT QUẢ Trang 84 Mảng ND - quảng lý sự tương ứng giữa bậc tự do phần tử và tổng quát Mảng CHT lưu chiều cao cột Mảng NDS - lưu địa chỉ Kii trên đường chéo chính trong SK Mảng SK vectơ một chiều lưu ma trận độ cứng tổng quát Mảng SM vectơ một chiều lưu ma trận khối lượng 4.1.3 Số liệu kết quả tính SLR Mảng D vectơ một chiều lưu trị riêng (hoặc các trị riêng) cần tính Mảng V mảng lưu trữ các vectơ riêng tương ứng. 4.2. Tổ chức chương trình và một số hàm cơ bản 4.2.1 Phát sinh kết cấu Chức năng này có nhiệm vụ phát sinh kết cấu dầm liên tục bao gồm: Nhận số liệu về chiều dài nhịp, số lượng gối... Tạo mảng toạ độ nút XY Tạo mảng bậc tự do nút JF (xác định các điều kiện biên) Tạo mảng LNC đánh số phần tử Hiển thị kết cấu (số thự tự nút, phần tử, điều kiện biên...) Hàm thực hiện các chức năng trên trong chương trình là hàm newFcn, hiennutFcn, hienptuFcn, hienkcFcn, hiendieukienbienFcn, hienthamsoFcn. 4.2.2 Xây dựng phương trình trị riêng 4.2.2.1 Xây dựng cấu trúc ma trận độ cứng và ma trận khối lượng Tạo các mảng NDF, ND, CHT, NDS theo SLV.JF và SLV.LNC, sơ đồ thực hiện như đã trình bày tại chương 2, mục 2.1.4 “thuật toán xây dựng và lưu trữ ma trận”. Để minh họa cho quá trình này, ta lấy ví dụ Trang 85 minh họa là 1 dầm liên tục 3 nhịp có sơ đồ 3x300m như hình dưới đây: Mảng SLV.Nut.JF - mảng bậc tự do nút JF u1 v1 qz1 1 1 1 0 2 0 1 0 3 0 1 0 4 0 1 0 Với ui,vi,fzi là chuyển vị theo phương ngang, phương đứng và góc xoay quanh trục z Mảng SLV.Pt.LNC - mảng chỉ số nút của phần tử LNC 1 2 1 1 2 2 2 3 3 3 4 Mảng NDF - Mảng đánh số bậc tự do nút NDF u1 v1 θz1 1 0 0 1 2 2 0 3 3 4 0 5 4 6 0 7 Tổng số bậc tự do nút hay tổng số phương trình cần giải là Neq=7 Trang 86 Mảng ND - mảng bậc tự do nút của phần tử ND u1 v1 θz1 u2 v2 θz2 1 0 0 1 2 0 3 2 2 0 3 4 0 5 3 4 0 5 6 0 7 Mảng CHT – mảng chiều cao cột PTử 1 2 3 CHT 1 0 0 2 1 0 1 3 2 1 2 4 2 0 2 5 3 1 3 6 2 2 7 3 3 Mảng NDS – quản lý ma trận độ cứng STT 1 2 3 4 5 6 7 CHT 0 1 2 2 3 2 3 NDS 1 2 4 7 10 14 17 21 Tính tổng số phương trình cần giải Neq(=7) Tính tổng số phần tử cần lưu của ma trận độ cứng Nsky(=20) Khai báo SK, SM 4.2.2.2 Xây dựng phương trình trị riêng Gọi thư viện phần tử, xây dựng tính chất phần tử Xếp vào phương trình trị riêng Hàm thực hiện chức năng này là XaydungPTTR 4.2.2.3 Xây dựng tính chất phần tử Chức năng chính của chức năng này là Trang 87 − Tính ma trận độ cứng [EK] Đối với phần tử dầm chịu uốn, như đã trình bày ở chương 1, ma trận độ cứng sẽ được tính theo công thức sau: [ ] . 4626 612612 2646 612612 ][][]][[][][ 22 22 3 ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎣ ⎡ − −−− − − = == ∫ ∫∫ LLLL LL LLLL LL L EJEKhay dFdxBBEdVBDBK z V TT e e Nếu xét đến cả chuyển vị đứng ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎢ ⎣ ⎡ − − − − == l EIDX l EI l EI l EF l EI l EI l EI l EI l EI l EI l EI l EI l EF EKK e 4 612 00 2604 6120612 0000 ][][ 23 2 2323 − Tính ma trận khối lượng [EM] ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎣ ⎡ −−− − − − = 22 22 422313 221561354 313422 135422156 420 ][ LLLL LL LLLL LL lEM ρ − Tính ma trận chuyển ET Trang 88 [ ] ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎣ ⎡ − − = 100000 0cossin000 0sincos000 000100 0000cossin 0000sincos ϕϕ ϕϕ ϕϕ ϕϕ T Hàm thực hiện chức năng này là: [EK,EM,ET]=Thanhdam(X,Y,E,F,IZ) 4.2.3 Giải bài toán trị riêng Bằng hàm nội tại của matlab – hàm [V,D] = eig(K,Q) Bằng phương pháp lặp ngược vectơ – hàm lapnguocFcn Bằng phương pháp xoay Jacobi – hàm xoayJacobiFcn Bằng phương pháp lặp không gian con – hàm SspaceFcn 4.3. Tính toán minh họa Với mục tiêu thử nghiệm các thuật toán giải bài toán trị riêng của hệ dầm liên tục. Tính toán minh họa được trình bày theo 2 ví dụ. Ví dụ 1 áp dụng cho bài toán nhỏ, ma trận độ cứng tổng quát K và ma trận khối lượng M được lưu đủ (biến lưu 2 ma trận này trong chương trình là KTQ và MTQ có kích thước là Neq x Neq, với Neq là tổng số bậc tự do của hệ. Bài toán trị riêng được giải theo thuật toán lặp ngược vectơ đến trị riêng nhỏ nhất và giải theo thuật toán xoay Jacobi. Kết quả được so sánh với kết quả lời giải của Matlab. Dữ liệu đầu vào: kết cấu dầm liên tục gồm 2 nhịp 2x3000 (cm), mỗi nhịp có 1 đầu ngàm 1 đầu tựa trên gối đỡ. Dầm có cùng tiết diện hình chữ nhật hxb=20x30 (cm) bằng vật liệu có moduyn đàn hồi là E=2,1 x 103 T/cm2 và mật độ khối lượng ρ=7.85 x 10-3 T/cm3. Trang 89 Phát sinh kết cấu trong matlab Kết quả (chi tiết xem tại “phụ lục 1 – ví dụ 1” phần phụ lục) Trị riêng nhận được bằng hàm nội tại của matlab: λ=60 (1/s2) Trị riêng nhận được qua phương pháp xoay Jacobi: λ =60 (1/s2) Trị riêng nhận được qua phương pháp lặp ngược vec tơ Lần lặp trị riêng sai số Lần lặp trị riêng sai số 1 127.8151 113.03% 9 59.1110 -1.48% 2 124.9819 108.30% 10 59.0250 -1.63% 3 121.7120 102.85% 11 59.0061 -1.66% 4 110.1729 83.62% 12 59.0019 -1.66% 5 86.8332 44.72% 13 59.0010 -1.67% 6 68.0412 13.40% 14 59.0008 -1.67% 7 61.2181 2.03% 15 59.0008 -1.67% 8 59.5002 -0.83% 16 59.0008 -1.67% Trang 90 Qua bảng so sánh trên có thể nhận thấy số lần lặp càng nhiều thì trị số trị riêng nhận được càng chính xác và càng gần với trị riêng mong muốn. Nếu sai số cho phép là dưới 5% thì ta có thể dừng ở bước lặp thứ 7. Nói chung, kết quả nhận được ở lần lặp >10 là chấp nhận được. Ví dụ 2 : giải cho bài toán lớn, ma trận độ cứng K và ma trận khối lượng M được lưu dưới dạng skyline (biến trong chương trình là SK và SM). Bài toán trị riêng được giải theo phương pháp lặp không gian con. Thông thường, theo kinh nghiệm của những người thiết kế, với mọi cơ hệ, ta chỉ xét đến khoảng 3 tần số dao động riêng nhỏ nhất và vectơ riêng tương ứng với nó là đủ. Nghĩa là số vectơ lặp nên sử dụng trong trường hợp này là q = min (2p;p+8) = min (6;11) = 6. Đây cũng là 1 điểm mạnh của phương pháp lặp không gian con, trong khi phương pháp lặp ngược vectơ chỉ cho ta trị riêng nhỏ nhất còn phương pháp xoay Jacobi lại cho ta tất cả các trị riêng của bài toán, thì lặp không gian con sẽ cho ta số trị riêng mong muốn, điều này sẽ làm tăng tốc độ tính toán lên đáng kể. Ta xét 3 trường hợp sau (kết quả chi tiết xem tại “phụ lục 2 – ví dụ 2” phần phụ lục): Trường hợp 1: mỗi nhịp dầm được chia làm 2 phần tử nhỏ (như trong ví dụ 1). Trang 91 Trường hợp 2: mỗi nhịp dầm được chia làm 4 phần tử nhỏ Trường hợp 3: mỗi nhịp dầm được chia làm 6 phần tử nhỏ Trang 92 Bảng so sánh kết quả tính theo phương pháp lặp không gian con và sử dụng hàm nội tại của matlab. Trường hợp 1 Trường hợp 2 Trường hợp 3 Matlab LKGC Matlab LKGC Matlab LKGC 60 59 58 57.9 130 126 122.3 122 830 831 Không chính xác 615.9 Không chính xác 609.9 Bảng so sánh kết quả tính theo phương pháp lặp không gian con và sử dụng phần mềm SAP 2000 TH 1 SAP 2000 Lặp không gian con Trị riêng 54.300 94.998 69647.644 59.000 126.000 831.000 V.tốc góc 7.369 9.747 263.908 7.681 11.225 28.827 Chu kỳ 0.853 0.645 0.024 0.818 0.560 0.218 Mô tả dạng dao động uốn phản đối xứng uốn đối xứng dọc trục phản đối xứng uốn phản đối xứng 1 bước sóng uốn đối xứng uốn phản đối xứng 2 bước sóng Trang 93 TH 2 SAP 2000 Lặp không gian con Trị riêng 58.714 123.037 594.736 58.000 122.300 615.900 V.tốc góc 7.663 11.092 24.387 7.616 11.059 24.817 Chu kỳ 0.820 0.566 0.258 0.825 0.568 0.253 Mô tả dạng dao động uốn phản đối xứng 1 bước sóng uốn đối xứng uốn phản đối xứng 2 bước sóng uốn phản đối xứng 1 bước sóng uốn đối xứng uốn phản đối xứng 2 bước sóng TH 3 SAP 2000 Lặp không gian con Trị riêng 58.822 123.721 614.335 57.900 122.000 609.900 V.tốc góc 7.670 11.123 24.786 7.609 11.045 24.696 Chu kỳ 0.819 0.565 0.253 0.826 0.569 0.254 Mô tả dạng dao động uốn phản đối xứng 1 bước sóng uốn đối xứng uốn phản đối xứng 2 bước sóng uốn phản đối xứng 1 bước sóng uốn đối xứng uốn phản đối xứng 2 bước sóng Qua ví dụ này có thể thấy, kết cấu càng được chia mịn thì trị riêng tính theo phương pháp lặp không gian con càng chính xác và ngược lại trị riêng tính theo bằng hàm nội tại của matlab càng cho sai số lớn. Tùy theo độ lớn của kết cấu mà số vòng lặp cần thiết trong phương pháp không gian con là ít hay nhiều, nhưng nói chung, chỉ sau số lần lặp nhỏ hơn rất nhiều so với phương pháp lặp vectơ ta đã nhận được kết quả mong muốn. Trang 94 KẾT LUẬN Với mục đích đặt ra ban đầu, luận văn đã giải quyết được các vấn đề cơ bản sau: − Đã khái quát được những phương pháp, những thuật toán chủ yếu sử dụng trong việc giải bài toán trị riêng, vectơ riêng. − Phân tích và xây dựng thuật toán giải bài toán trị riêng bằng phương pháp phần tử hữu hạn, đặc biệt coi trọng phương pháp tổ chức xây dựng và lưu trữ ma trận độ cứng, một vấn đề có tính chất quyết định đến tính khả thi khi giải quyết các bài toán lớn có sử dụng các chương trình tính toán. − Sử dụng ngôn ngữ lập trình Matlab để khảo sát bài toán dầm liên tục. Giải bài toán trị riêng theo các phương pháp : lặp ngược vectơ, xoay Jacobi, lặp không gian con. Trong khuôn khổ của một luận văn thạc sĩ, thời lượng có hạn, trình độ còn hạn chế, học viên cũng mới chỉ thực hiện được bước đầu trong nghiên cứu. Nếu có điều kiện tiếp tục nghiên cứu, học viên sẽ đi sâu vào giải quyết các vấn đề tính toán động lực học kết cấu cầu, dạng thanh thành mỏng và hệ dây. Trang 95 DANH MỤC TÀI LIỆU THAM KHẢO Tiếng Việt 1. NguyÔn Quèc B¶o, TrÇn NhÊt Dòng (2002), Ph−¬ng ph¸p phÇn tö h÷u h¹n lý thuyÕt vµ lËp tr×nh, NXB Khoa häc vµ Kü thuËt, Hµ Néi. 2. Ph¹m §×nh Ba, NguyÔn V¨n Hîi (1994), Gi¸o tr×nh ®éng lùc häc c«ng tr×nh, Häc viÖn kü thuËt qu©n sù, Hµ Néi. 3. Ph¹m §×nh Ba (1995), Bµi tËp ®éng lùc häc c«ng tr×nh, Häc viÖn kü thuËt qu©n sù, Hµ Néi. 4. Ph¹m Kh¾c Hïng, Lª V¨n Quý, LÒu Thä Tr×nh (1974), æn ®Þnh vµ ®éng lùc häc c«ng tr×nh, NXB §¹i häc vµ Trung häc chuyªn nghiÖp, Hµ Néi. 5. Lª Quúnh Mai (2000), TÝnh dao ®éng riªng cña kÕt cÊu d¹ng dÇm b»ng ph−¬ng ph¸p ma trËn chuyÓn tiÕp, LuËn v¨n th¹c sü khoa häc kü thuËt, §¹i häc giao th«ng vËn t¶i. 6. NguyÔn Xu©n Ngäc, NguyÔn Tµi Trung (1997), æn ®Þnh vµ ®éng lùc häc c«ng tr×nh, NXB x©y dùng. 7. NguyÔn Ngäc Quúnh, Hå ThuÇn (1978), øng dông ma trËn trong kü thuËt, NXB Khoa häc vµ Kü thuËt. 8. NguyÔn V¨n TØnh (1987), C¬ së tÝnh dao ®éng c«ng tr×nh, NXB x©y dùng. 9. Bïi §øc Vinh (2001), Ph©n tÝch vµ thiÕt kÕ kÕt cÊu b»ng phÇn mÒm SAP 2000, NXB Thèng kª. Trang 96 Tiếng Anh 1. Klaus-Jurgen Bathe, professor of Mechanical Engineering, Massachusets Institue of Technology (1992), Finite element procedures, Prentice Hall, New Jersey. 2. Ray W.Clough, Joseph Penzien (1993), Dynamics of Structures, MacGrawhill Inc. 3. CSI – computers and structures inc (July 1997-revised), SAP2000 integrated finite element analysis and design of structures, Backeley, California, USA. PHỤ LỤC PHỤ LỤC 1 – VÍ DỤ 1 Ma trận khối lượng tổng quát của hệ MTQ = 0.4788 0 0 0.1197 0 0 0 0 0 0.5335 0 0 -0.3335 0 0 0 0 0 3.0780 0 -1.1543 0 0 0 0.1197 0 0 0.4788 0 0.1197 0 0 0 -0.3335 -1.1543 0 3.0780 0 0.3335 -1.1543 0 0 0 0.1197 0 0.4788 0 0 0 0 0 0 0.3335 0 0.5335 0 0 0 0 0 -1.1543 0 0 3.0780 Ma trận độ cứng tổng quát của hệ KTQ = 1.0e+005 * 1.6800 0 0 -0.8400 0 0 0 0 0 0.0007 0 0 0.0025 0 0 0 0 0 0.0504 0 0.0126 0 0 0 -0.8400 0 0 1.6800 0 -0.8400 0 0 0 0.0025 0.0126 0 0.0504 0 -0.0025 0.0126 0 0 0 -0.8400 0 1.6800 0 0 0 0 0 0 -0.0025 0 0.0007 0 0 0 0 0 0.0126 0 0 0.0504 Trị riêng, vectơ riêng giải bằng hàm nội tại của matlab D = 1.0e+005 * 0.0006 0 0 0 0 0 0 0 0 0.0013 0 0 0 0 0 0 0 0 0.0083 0 0 0 0 0 0 0 0 0.0164 0 0 0 0 0 0 0 0 0.0590 0 0 0 0 0 0 0 0 0.7593 0 0 0 0 0 0 0 0 3.5088 0 0 0 0 0 0 0 0 9.2658 V = 0.0000 0.0000 -0.0000 -0.0000 0.0000 -0.6211 -1.0219 0.8987 0.8686 -0.9681 0.3999 -0.0000 0.4902 0.0000 0.0000 0.0000 0.0312 -0.0000 -0.2543 -0.4030 0.4185 0.0000 -0.0000 -0.0000 0.0000 0.0000 -0.0000 -0.0000 0.0000 -0.8784 -0.0000 -1.2710 -0.1142 -0.0000 0.2844 -0.0000 0.6805 -0.0000 0.0000 0.0000 -0.0000 -0.0000 -0.0000 -0.0000 0.0000 -0.6211 1.0219 0.8987 -0.8686 -0.9681 -0.3999 0.0000 -0.4902 0.0000 -0.0000 -0.0000 0.0312 0.0000 -0.2543 0.4030 0.4185 0.0000 0.0000 0.0000 Trị riêng, vectơ riêng giải bằng phương pháp lặp ngược vec tơ LANLAP = 1 LAMDA = 127.8151 VK = 0.0012 0.9179 0.0097 0.0015 0.0129 0.0012 1.0145 0.0097 LANLAP = 2 LAMDA = 124.9819 VK = 0.0000 0.8550 -0.0030 0.0000 0.0133 0.0000 1.0669 -0.0030 LANLAP = 3 LAMDA = 121.7120 VK = 0.0000 0.7181 -0.0078 0.0000 0.0286 0.0000 1.1557 -0.0078 LANLAP = 4 LAMDA = 110.1729 VK = 0.0000 0.4246 -0.0151 0.0000 0.0554 0.0000 1.2681 -0.0151 LANLAP = 5 LAMDA = 86.8332 VK = 0.0000 -0.0398 -0.0239 0.0000 0.0873 0.0000 1.2881 -0.0239 LANLAP = 6 LAMDA = 68.0412 VK = 0.0000 -0.4521 -0.0290 0.0000 0.1062 0.0000 1.1636 -0.0290 LANLAP = 7 LAMDA = 61.2181 VK = 0.0000 -0.6779 -0.0307 0.0000 0.1123 0.0000 1.0303 -0.0307 LANLAP = 8 LAMDA = 59.5002 VK = 0.0000 -0.7818 -0.0311 0.0000 0.1138 0.0000 0.9490 -0.0311 LANLAP = 9 LAMDA = 59.1110 VK = 0.0000 -0.8286 -0.0312 0.0000 0.1141 0.0000 0.9072 -0.0312 LANLAP = 10 LAMDA = 59.0250 VK = 0.0000 -0.8500 -0.0312 0.0000 0.1142 0.0000 0.8869 -0.0312 LANLAP = 11 LAMDA = 59.0061 VK = 0.0000 -0.8600 -0.0312 0.0000 0.1142 0.0000 0.8772 -0.0312 LANLAP = 12 LAMDA = 59.0019 VK = 0.0000 -0.8646 -0.0312 0.0000 0.1142 0.0000 0.8726 -0.0312 LANLAP = 13 LAMDA = 59.0010 VK = 0.0000 -0.8667 -0.0312 0.0000 0.1142 0.0000 0.8705 -0.0312 LANLAP = 14 LAMDA = 59.0008 VK = 0.0000 -0.8677 -0.0312 0.0000 0.1142 0.0000 0.8695 -0.0312 LANLAP = 15 LAMDA = 59.0008 VK = 0.0000 -0.8682 -0.0312 0.0000 0.1142 0.0000 0.8690 -0.0312 LANLAP = 16 LAMDA = 59.0008 VK = 0.0000 -0.8684 -0.0312 0.0000 0.1142 0.0000 0.8688 -0.0312 KIEMTRA2 = 1.0000 Trị riêng, vectơ riêng giải bằng phương pháp xoay Jacobi V = 0.8987 0 0 0.6211 0 -1.0219 0 0 0 0.8686 -0.3999 0 0.4902 0 0.9681 0.0000 0 0.0312 0.2543 0 0.4185 0 0.0000 -0.4030 -1.2710 0 0 0.8784 0 -0.0000 0 0 0 -0.1142 -0.2844 0 0.6805 0 -0.0000 -0.0000 0.8987 0 0 0.6211 0 1.0219 0 0 0 -0.8686 0.3999 0 -0.4902 0 0.9681 -0.0000 0 0.0312 0.2543 0 0.4185 0 0.0000 0.4030 D = 1.0e+005 * 9.2658 0.0006 0.0083 0.7593 0.0590 3.5088 0.0013 0.0164 PHỤ LỤC 2 – VÍ DỤ 2 Trường hợp 1: mỗi nhịp được chia thành 2 phần tử nhỏ Trị riêng nhận được bằng hàm nội tại của matlab D = 1.0e+005 * 0.0006 0 0 0 0 0 0 0 0 0.0013 0 0 0 0 0 0 0 0 0.0083 0 0 0 0 0 0 0 0 0.0164 0 0 0 0 0 0 0 0 0.0590 0 0 0 0 0 0 0 0 0.7593 0 0 0 0 0 0 0 0 3.5088 0 0 0 0 0 0 0 0 9.2658 Giải bằng phương pháp lặp không gian con Vectơ lặp xuất phát Xstat = 1.0e+005 * 0.4788 0 0 0 0 0.3412 0.5335 0 0 0 0 0.5341 3.0780 0 0 0 0 0.7271 0.4788 0 0 0 0 0.3093 3.0780 0 0 0 0 0.8385 0.4788 0 0 0 0 0.5681 0.5335 0 0 0 0 0.3704 3.0780 0 0 0 0 0.7027 Trị riêng trong mỗi lần lặp Buoclap = 1 D = 1.0e+004 * 7.9881 0.0059 0.5903 0.1638 0.0831 0.0126 Buoclap = 2 D = 1.0e+004 * 7.6028 0.0059 0.5902 0.1637 0.0831 0.0126 Buoclap = 3 D = 1.0e+004 * 7.5930 0.0059 0.5902 0.1637 0.0831 0.0126 Buoclap = 4 D = 1.0e+004 * 7.5926 0.0059 0.5902 0.1637 0.0831 0.0126 Buoclap = 5 D = 1.0e+004 * 7.5926 0.0059 0.5902 0.1637 0.0831 0.0126 Buoclap = 6 D = 1.0e+004 * 7.5926 0.0059 0.5902 0.1637 0.0831 0.0126 Buoclap = 7 D = 1.0e+004 * 7.5926 0.0059 0.5902 0.1637 0.0831 0.0126 Buoclap = 8 D = 1.0e+004 * 7.5926 0.0059 0.5902 0.1637 0.0831 0.0126 Buoclap = 9 D = 1.0e+004 * 7.5926 0.0059 0.5902 0.1637 0.0831 0.0126 Buoclap = 10 D = 1.0e+004 * 7.5926 0.0059 0.5902 0.1637 0.0831 0.0126 Buoclap = 11 D = 1.0e+004 * 7.5926 0.0059 0.5902 0.1637 0.0831 0.0126 Buoclap = 12 D = 1.0e+004 * 7.5926 0.0059 0.5902 0.1637 0.0831 0.0126 Buoclap = 13 D = 1.0e+004 * 7.5926 0.0059 0.5902 0.1637 0.0831 0.0126 Buoclap = 14 D = 1.0e+004 * 7.5926 0.0059 0.5902 0.1637 0.0831 0.0126 Buoclap = 15 D = 1.0e+004 * 7.5926 0.0059 0.5902 0.1637 0.0831 0.0126 Buoclap = 16 D = 1.0e+004 * 7.5926 0.0059 0.5902 0.1637 0.0831 0.0126 Trị riêng nhận được bằng phần mềm SAP 2000 PROGRAM SAP2000 - VERSION N6.11 FILE:CHIA 2.OUT MODE PERIOD FREQUENCY FREQUENCY EIGENVALUE (TIME) (CYC/TIME) (RAD/TIME) (RAD/TIME)**2 1 0.852603 1.172879 7.369415 54.308278 2 0.644647 1.551236 9.746704 94.998236 3 0.023808 42.002327 263.908402 69647.644 4 0.012885 77.610180 487.639141 237791.932 5 0.009862 101.402587 637.131242 405936.219 Trường hợp 2: mỗi nhịp được chia thành 4 phần tử nhỏ Trị riêng nhận được bằng hàm nội tại của matlab D = 1.0e+006 * 0.0001 0 0 0 0 0 0 0 0 0 0 0 0 0.0001 0 0 0 0 0 0 0 0 0 0 0 0 0.0006 0 0 0 0 0 0 0 0 0 0 0 0 0.0009 0 0 0 0 0 0 0 0 0 0 0 0 0.0028 0 0 0 0 0 0 0 0 0 0 0 0 0.0037 0 0 0 0 0 0 0 0 0 0 0 0 0.0098 0 0 0 0 0 0 0 0 0 0 0 0 0.0133 0 0 0 0 0 0 0 0 0 0 0 0 0.0260 0 0 0 0 0 0 0 0 0 0 0 0 0.0364 0 0 0 0 0 0 0 0 0 0 0 0 0.0662 0 0 0 0 0 0 0 0 0 0 0 0 0.0731 Giải bằng phương pháp lặp không gian con Vectơ lặp xuất phát Xstat = 0.2394 0 0 0 0 0.5466 0.2668 0 0 0 0 0.4449 0.3848 0 0 0 0 0.6946 0.2394 0 0 0 0 0.6213 0.2668 0 0 0 0 0.7948 0.3848 0 0 0 0 0.9568 0.2394 0 0 0 0 0.5226 0.2668 0 0 0 0 0.8801 0.3848 0 0 0 0 0.1730 0.2394 0 0 0 0 0.9797 0.3848 0 0 0 0 0.2714 0.2394 0 0 0 0 0.2523 0.2668 0 0 0 0 0.8757 0.3848 0 0 0 0 0.7373 0.2394 0 0 0 0 0.1365 0.2668 0 0 0 0 0.0118 0.3848 0 0 0 0 0.8939 0.2394 0 0 0 0 0.1991 0.2668 0 0 0 0 0.2987 0.3848 0 0 0 0 0.6614 Trị riêng trong mỗi lần lặp Buoclap = 1 D = 1.0e+003 * 0.1229 1.0225 0.0581 3.2858 0.6407 8.1754 Buoclap = 2 D = 1.0e+003 * 0.1223 0.9442 0.0580 2.7803 0.6159 3.7395 Buoclap = 3 D = 1.0e+003 * 0.1223 0.9440 0.0580 2.7692 0.6159 3.7171 Buoclap = 4 D = 1.0e+003 * 0.1223 0.9440 0.0580 2.7688 0.6159 3.7157 Buoclap = 5 D = 1.0e+003 * 0.1223 0.9440 0.0580 2.7687 0.6159 3.7156 Buoclap = 6 D = 1.0e+003 * 0.1223 0.9440 0.0580 2.7687 0.6159 3.7156 Buoclap = 7 D = 1.0e+003 * 0.1223 0.9440 0.0580 2.7687 0.6159 3.7156 Buoclap = 8 D = 1.0e+003 * 0.1223 0.9440 0.0580 2.7687 0.6159 3.7156 Buoclap = 9 D = 1.0e+003 * 0.1223 0.9440 0.0580 2.7687 0.6159 3.7156 Buoclap = 10 D = 1.0e+003 * 0.1223 0.9440 0.0580 2.7687 0.6159 3.7156 Buoclap = 11 D = 1.0e+003 * 0.1223 0.9440 0.0580 2.7687 0.6159 3.7156 Buoclap = 12 D = 1.0e+003 * 0.1223 0.9440 0.0580 2.7687 0.6159 3.7156 Buoclap = 13 D = 1.0e+003 * 0.1223 0.9440 0.0580 2.7687 0.6159 3.7156 Buoclap = 14 D = 1.0e+003 * 0.1223 0.9440 0.0580 2.7687 0.6159 3.7156 Buoclap = 15 D = 1.0e+003 * 0.1223 0.9440 0.0580 2.7687 0.6159 3.7156 Buoclap = 16 D = 1.0e+003 * 0.1223 0.9440 0.0580 2.7687 0.6159 3.7156 Trị riêng nhận được bằng phần mềm SAP 2000 PROGRAM SAP2000 - VERSION N6.11 FILE:CHIA 4.OUT MODE PERIOD FREQUENCY FREQUENCY EIGENVALUE (TIME) (CYC/TIME) (RAD/TIME) (RAD/TIME)**2 1 0.819991 1.219525 7.662501 58.713925 2 0.566450 1.765382 11.092221 123.037371 3 0.257643 3.881344 24.387206 594.735813 4 0.213410 4.685806 29.441790 866.818986 5 0.138224 7.234611 45.456402 2066.284 6 0.129933 7.696294 48.357244 2338.423 7 0.023351 42.825201 269.078673 72403.332 8 0.011904 84.004653 527.816803 278590.578 9 0.008200 121.955854 766.271230 587171.598 10 0.006442 155.220359 975.278282 951167.727 Trường hợp 3: mỗi nhịp được chia thành 6 phần tử nhỏ Trị riêng nhận được bằng hàm nội tại của matlab D = 1.0e+007 * 0.0000 0 0 0 0 0 0 0 0 0 0 0 0 0.0000 0 0 0 0 0 0 0 0 0 0 0 0 0.0001 0 0 0 0 0 0 0 0 0 0 0 0 0.0001 0 0 0 0 0 0 0 0 0 0 0 0 0.0003 0 0 0 0 0 0 0 0 0 0 0 0 0.0004 0 0 0 0 0 0 0 0 0 0 0 0 0.0008 0 0 0 0 0 0 0 0 0 0 0 0 0.0010 0 0 0 0 0 0 0 0 0 0 0 0 0.0019 0 0 0 0 0 0 0 0 0 0 0 0 0.0023 0 0 0 0 0 0 0 0 0 0 0 0 0.0045 0 0 0 0 0 0 0 0 0 0 0 0 0.0055 Giải bằng phương pháp lặp không gian con Vectơ lặp xuất phát Xstat = 0.1596 0 0 0 0 0.0498 0.1778 0 0 0 0 0.0784 0.1140 0 0 0 0 0.6408 0.1596 0 0 0 0 0.1909 0.1778 0 0 0 0 0.8439 0.1140 0 0 0 0 0.1739 0.1596 0 0 0 0 0.1708 0.1778 0 0 0 0 0.9943 0.1140 0 0 0 0 0.4398 0.1596 0 0 0 0 0.3400 0.1778 0 0 0 0 0.3142 0.1140 0 0 0 0 0.3651 0.1596 0 0 0 0 0.3932 0.1778 0 0 0 0 0.5915 0.1140 0 0 0 0 0.1197 0.1596 0 0 0 0 0.0381 0.1140 0 0 0 0 0.4586 0.1596 0 0 0 0 0.8699 0.1778 0 0 0 0 0.9342 0.1140 0 0 0 0 0.2644 0.1596 0 0 0 0 0.1603 0.1778 0 0 0 0 0.8729 0.1140 0 0 0 0 0.2379 0.1596 0 0 0 0 0.6458 0.1778 0 0 0 0 0.9669 0.1140 0 0 0 0 0.6649 0.1596 0 0 0 0 0.8704 0.1778 0 0 0 0 0.0099 0.1140 0 0 0 0 0.1370 0.1596 0 0 0 0 0.8188 0.1778 0 0 0 0 0.4302 0.1140 0 0 0 0 0.8903 Trị riêng trong mỗi lần lặp Buoclap = 1 D = 1.0e+004 * 0.0123 1.5070 0.0989 0.3781 0.0058 0.0650 Buoclap = 2 D = 1.0e+003 * 0.1220 9.5920 0.9312 2.9789 0.0579 0.6107 Buoclap = 3 D = 1.0e+003 * 0.1220 6.2618 0.9305 2.8184 0.0579 0.6099 Buoclap = 4 D = 1.0e+003 * 0.1220 4.2790 0.9304 2.7242 0.0579 0.6099 Buoclap = 5 D = 1.0e+003 * 0.1220 3.7423 0.9304 2.6846 0.0579 0.6099 Buoclap = 6 D = 1.0e+003 * 0.1220 3.6377 0.9304 2.6776 0.0579 0.6099 Buoclap = 7 D = 1.0e+003 * 0.1220 3.6178 0.9304 2.6768 0.0579 0.6099 Buoclap = 8 D = 1.0e+003 * 0.1220 3.6140 0.9304 2.6767 0.0579 0.6099 Buoclap = 9 D = 1.0e+003 * 0.1220 3.6132 0.9304 2.6767 0.0579 0.6099 Buoclap = 10 D = 1.0e+003 * 0.1220 3.6131 0.9304 2.6767 0.0579 0.6099 Buoclap = 11 D = 1.0e+003 * 0.1220 3.6130 0.9304 2.6767 0.0579 0.6099 Buoclap = 12 D = 1.0e+003 * 0.1220 3.6130 0.9304 2.6767 0.0579 0.6099 Buoclap = 13 D = 1.0e+003 * 0.1220 3.6130 0.9304 2.6767 0.0579 0.6099 Buoclap = 14 D = 1.0e+003 * 0.1220 3.6130 0.9304 2.6767 0.0579 0.6099 Buoclap = 15 D = 1.0e+003 * 0.1220 3.6130 0.9304 2.6767 0.0579 0.6099 Buoclap = 16 D = 1.0e+003 * 0.1220 3.6130 0.9304 2.6767 0.0579 0.6099 Trị riêng nhận được bằng phần mềm SAP 2000 PROGRAM SAP2000 - VERSION N6.11 FILE:CHIA 6.OUT MODE PERIOD FREQUENCY FREQUENCY EIGENVALUE (TIME) (CYC/TIME) (RAD/TIME) (RAD/TIME)**2 1 0.819240 1.220643 7.669528 58.821663 2 0.564882 1.770282 11.123008 123.721314 3 0.253500 3.944780 24.785782 614.334972 4 0.205942 4.855733 30.509473 930.827941 5 0.123037 8.127645 51.067499 2607.889 6 0.107281 9.321290 58.567390 3430.139 7 0.075271 13.285343 83.474269 6967.954 8 0.069473 14.394022 90.440308 8179.449 9 0.055982 17.862917 112.236018 12596.924 10 0.054623 18.307333 115.028367 13231.525