Luận án Nghiên cứu dao động ngẫu nhiên phi tuyến bằng tiêu chuẩn sai số bình phương trung bình địa phương – Tổng thể

BỘ GIÁO DỤC VÀ ĐÀO TẠO VIỆN HÀN LÂM KHOA HỌC 
VÀ CÔNG NGHỆ VIỆT NAM 
HỌC VIỆN KHOA HỌC VÀ CÔNG NGHỆ 
----------------------------- 
Nguyễn Cao Thắng 
NGHIÊN CỨU DAO ĐỘNG NGẪU NHIÊN PHI TUYẾN 
BẰNG TIÊU CHUẨN SAI SỐ BÌNH PHƯƠNG TRUNG BÌNH 
ĐỊA PHƯƠNG – TỔNG THỂ 
LUẬN ÁN TIẾN SĨ KỸ THUẬT CƠ KHÍ VÀ CƠ KỸ THUẬT 
Hà Nội – 2019 
BỘ GIÁO DỤC VÀ ĐÀO TẠO VIỆN HÀN LÂM KHOA HỌC 
VÀ CÔNG NGHỆ VIỆT NAM 
HỌC VIỆN KHOA HỌC VÀ CÔNG NGHỆ 
----------------------------- 
Nguyễn Cao Thắng 
NGHIÊN CỨU DAO ĐỘNG NGẪU NHIÊN PHI TUYẾN 
BẰNG TIÊU CHUẨN SAI SỐ BÌNH PHƯƠNG TRUNG BÌNH 
ĐỊA PHƯƠNG – TỔNG THỂ 
 Chuyên ngành: Cơ kỹ thuật 
 Mã số: 9 52 01 01 
LUẬN ÁN TIẾN SỸ KỸ THUẬT CƠ KHÍ VÀ CƠ KỸ THUẬT 
 NGƯỜI HƯỚNG DẪN KHOA HỌC: 
 1. TS. Lưu Xuân Hùng 
 2. GS. TSKH. Nguyễn Đông Anh 
Hà Nội – 2019 
I 
LỜI CAM ĐOAN 
Tôi xin cam đoan đây là công trình nghiên cứu của riêng tôi, các kết quả nghiên cứu 
được trình bày trong luận án là trung thực, khách quan và chưa từng được bảo vệ ở 
bất kỳ học vị nào. 
Tôi xin cam đoan rằng mọi sự giúp đỡ cho việc thực hiện luận án đã được cảm ơn, 
các thông tin trích dẫn trong luận án này đều được chỉ rõ nguồn gốc. 
 Tác giả luận án 
 Nguyễn Cao Thắng 
II 
LỜI CÁM ƠN 
Tôi xin gửi lời cám ơn chân thành và sâu sắc tới các thầy hướng dẫn khoa học, TS. 
Lưu Xuân Hùng và GS.TSKH. Nguyễn Đông Anh, đã tận tâm hướng dẫn khoa học, 
chỉ bảo tôi nghiên cứu và giúp đỡ tôi hoàn thành luận án. 
Tôi xin bày tỏ sự cám ơn tới Phòng Cơ học Công Trình, Viện Cơ học, Học Viện 
Khoa học và Công nghệ, các thầy cô đã tạo điều kiện đã giúp đỡ tôi ngay từ những 
ngày đầu làm luận án. 
Tôi xin cám ơn các đồng nghiệp, TS. Nguyễn Như Hiếu, TS. Nguyễn Văn Hải và 
nhiều người khác đã hỗ trợ, động viên tôi trong nhiều lúc khó khăn. 
Cuối cùng tôi xin bày tỏ sự biết ơn đến gia đình đã giúp đỡ, ủng hộ tôi trong suốt 
thời gian làm luận án. 
III 
MỤC LỤC 
LỜI CAM ĐOAN ............................................................................................................. I 
LỜI CÁM ƠN ................................................................................................................... II 
MỤC LỤC ......................................................................................................................... III 
DANH MỤC CÁC KÝ HIỆU, CÁC CHỮ VIẾT TẮT .................................................... VI 
DANH MỤC CÁC HÌNH VẼ, ĐỒ THỊ .......................................................................... IX 
DANH MỤC CÁC BẢNG................................................................................................ X 
MỞ ĐẦU ..................................................................................................................... 1 
CHƯƠNG 1. GIỚI THIỆU TỔNG QUAN MỘT SỐ PHƯƠNG PHÁP PHÂN 
TÍCH DAO ĐỘNG NGẪU NHIÊN PHI TUYẾN .......................................................... 6 
1.1. Đại lượng ngẫu nhiên và các đặc trưng xác suất ....................................................... 6 
1.2 Quá trình ngẫu nhiên ................................................................................................... 8 
1.3 Một số quá trình ngẫu nhiên đặc biệt .......................................................................... 11 
1.4 Một số phương pháp giải tích gần đúng phân tích dao động ngẫu nhiên ................... 16 
1.5 Phương pháp phương trình Fokker-Planck-Kolmogorov (FPK) và phương pháp 
trung bình hóa ngẫu nhiên ................................................................................................. 19 
1.6. Tổng quan một số nghiên cứu về dao động ngẫu nhiên ............................................ 25 
Kết luận chương 1 ............................................................................................................. 28 
CHƯƠNG 2. PHƯƠNG PHÁP TUYẾN TÍNH HÓA TƯƠNG ĐƯƠNG VÀ TIÊU 
CHUẨN SAI SỐ BÌNH PHƯƠNG TRUNG BÌNH ĐỊA PHƯƠNG – TỔNG THỂ ....... 29 
2.1. Tiêu chuẩn tuyến tính hóa tương đương kinh điển .................................................... 29 
2.2. Một số tiêu chuẩn tuyến tính hóa tương đương cải tiến ............................................ 36 
2.2.1 Tiêu chuẩn tuyến tính hóa dựa trên sai số thế năng ................................................. 37 
2.2.2 Tiêu chuẩn tuyến tính hóa tương đương có điều chỉnh ............................................ 38 
2.2.3 Tiêu chuẩn tuyến tính hóa tương đương đối ngẫu ................................................... 39 
IV 
2.3 Tiêu chuẩn sai số bình phương trung bình địa phương-tổng thể (GLOMSEC) ......... 40 
Kết luận chương 2 ............................................................................................................. 49 
CHƯƠNG 3. ỨNG DỤNG TIÊU CHUẨN GLOMSEC TRONG PHÂN TÍCH 
CÁC HỆ DAO ĐỘNG NGẪU NHIÊN PHI TUYẾN MỘT BẬC TỰ DO .................... 50 
3.1. Phân tích miền tập trung đáp ứng của hệ dao động phi tuyến ................................... 50 
3.1.1. Hệ dao dộng Duffing chịu kích động ồn trắng ....................................................... 50 
3.1.2.Hệ dao dộng có cản phi tuyến chịu kích động ồn trắng ........................................... 52 
3.1.3. Hệ dao dộng có lực cản và đàn hồi phi tuyến chịu kích động ồn trắng .................. 54 
3.2. Các ví dụ ứng dụng tiêu chuẩn sai số bình phương trung bình địa phương-tổng 
thể (GLOMSEC) ............................................................................................................... 57 
3.2.1 Dao động có cản phi tuyến bậc ba ........................................................................... 57 
3.2.2. Dao động trong hệ Van der Pol với kích động ngẫu nhiên ..................................... 60 
3.2.3 Dao động trong hệ Duffing với kích động ngẫu nhiên ........................................... 63 
3.2.4. Hệ dao động Duffing với cản phi tuyến chịu kích động ngẫu nhiên ...................... 66 
3.2.5. Dao động của tàu thủy ............................................................................................ 70 
Kết luận chương 3 ............................................................................................................. 72 
CHƯƠNG 4. ỨNG DỤNG TIÊU CHUẨN GLOMSEC TRONG PHÂN TÍCH 
CÁC HỆ DAO ĐỘNG NGẪU NHIÊN PHI TUYẾN NHIỀU BẬC TỰ DO ................. 73 
4.1. Hệ dao động phi tuyến hai bậc tự do ......................................................................... 73 
4.2. Hệ dao động phi tuyến chịu kích động ngẫu nhiên ồn màu ....................................... 80 
4.2.1. Mở rộng GLOMSEC cho trường hợp chịu kích động ngẫu nhiên ồn màu ............ 80 
4.2.2. Hệ dao động Duffing chịu kích động ngẫu nhiên ồn màu ...................................... 83 
4.2.3. Hệ dao động Duffing có cản phi tuyến chịu kích động ngẫu nhiên ồn màu ........... 86 
Kết luận Chương 4 ............................................................................................................ 88 
KẾT LUẬN VÀ KIẾN NGHỊ .......................................................................................... 90 
DANH SÁCH CÁC CÔNG TRÌNH ĐÃ CÔNG BỐ CỦA LUẬN ÁN........................... 92 
V 
TÀI LIỆU THAM KHẢO ................................................................................................. 93 
PHỤ LỤC .......................................................................................................................... 100 
VI 
 DANH MỤC CÁC KÝ HIỆU, CÁC CHỮ VIẾT TẮT 
GLOMSEC (GL) tiêu chuẩn sai số bình phương trung bình địa 
 phương tổng thể (Global-Local Mean Square Error 
 Criterion) 
LOMSEC (L) tiêu chuẩn sai số bình phương trung bình địa 
 phương (Local Mean Square Error Criterion) 
TTH tuyến tính hóa 
FPK phương trình Fokker-Planck-Kolgomorov 
ENL phi tuyến hóa tương đương (Equivalent Non – 
 Linearization) 
kd kinh điển 
MC mô phỏng Monte Carlo 
PDF hàm mật độ xác xuất (Probability Density 
 Function) 
SDOF hệ một bậc tự do 
MDOF hệ nhiều bậc tự do 
NL năng lượng 
M ma trận khối lượng 
K ma trận hệ số độ cứng 
C ma trận hệ số cản 
( )α  ma trận đáp ứng tần số 
( )wS  ma trận mật độ phổ của véc-tơ w(t) 
a, r biến không thứ nguyên dương 
VII 
, , , ,     hệ số dương 
b, k, c hệ số tuyến tính hóa tương đương 
,i jb k hệ số tuyến tính hóa thành phần thứ i, thứ j 
h hệ số cản tuyến tính 
C hệ số chuẩn hóa 
1 , ttc k hệ số độ cứng tuyến tính 
 1 2 12, ,xxD t t D hiệp phương sai 
 x hàm Delta Dirac 
 ,E   kỳ vọng toán 
 ,e x x sai số phương trình 
 F x hàm phân phối xác suất 
 ,f t u t kích động ngoài 
 ,g x x hàm phi tuyến của dịch chuyển và vận tốc 
 ,H x x hàm tổng năng lượng 
 ,K x t ma trận hệ số khuyếch tán 
 1 2,R t t hàm tương quan 
m khối lượng 
xm trung bình xác suất 
minS giá trị cực tiểu của tiêu chuẩn tuyến tính hóa 
n mô men trung tâm 
nm mô men liên kết trung tâm 
 P  xác suất của một sự kiện 
VIII 
 , ,p x p x x hàm mật độ xác suất một chiều, hai chiều 
 0 0, ,p x t x t mật độ xác suất chuyển tiếp 
 xS  hàm mật độ phổ 
0S mật độ phổ hằng số 
T chu kỳ dao động 
0 1 2, , ,t t t t thời gian 
 độ trễ 
 U x hàm thế năng 
u, v véc tơ 
 ,v t x t vận tốc 
X, Y biến ngẫu nhiên 
 x t dịch chuyển 
 x t gia tốc 
 t quá trình Wiener 
 t quá trình ồn trắng 
 cường độ của ồn trắng 
x độ lệch chuẩn 
2
x phương sai 
 tần số của kích động 
0 tần số dao động tự do 
IX 
DANH MỤC CÁC HÌNH VẼ, ĐỒ THỊ 
Hình 1.1. Các hàm mẫu theo thời gian của quá trình ngẫu nhiên .................................... 9 
Hình 1.2. Các hàm mật độ xác suất chuẩn ...................................................................... 14 
Hình 3.1. Đồ thị hàm PDF p(x) của hệ Duffing, ( =0.1, 1, 10, 100) ............................ 52 
Hình 3.2. Đồ thị hàm PDF ,p x x
của hệ cản phi tuyến, (=0.1, 1, 10, 100) ............... 54 
Hình 3.3. Đồ thị hàm PDF ,p x x của hệ cản và đàn hồi phi tuyến, (=0.1, 1, 10, 
100) ................................................................................................................................. 56 
X 
DANH MỤC CÁC BẢNG 
Bảng 3.1. Các giá trị của a phụ thuộc theo  ......................................................... 51 
Bảng 3.2. Các giá trị của a phụ thuộc theo  ......................................................... 53 
Bảng 3.3. Các giá trị của a phụ thuộc theo  ......................................................... 55 
Bảng 3.4. Momen bậc hai của đáp ứng của hệ dao động cản phi tuyến với 
0.05, 1, 4oh h  , và γ thay đổi..................................................................... 60 
Bảng 3.5 Đáp ứng bình phương trung bình của dao động Van der Pol với 
α*ε=0.2; 0 =1; *  =2; σ
2 thay đổi ................................................................ 63 
Bảng 3.6 Đáp ứng bình phương trung bình của hệ dao động Duffing với 
1, 0.25, 1o h  ; hệ số đàn hồi phi tuyến  thay đổi .................................... 65 
Bảng 3.7. Đáp ứng bình phương trung bình của hệ dao động Duffing với cản phi 
tuyến ( 1.0, 0.1, 2   ;  thay đổi) ........................................................ 69 
Bảng 4.1. Các mô men bậc hai của đáp ứng của 1 2,x x theo 1 2 với 
1 2 1 2 0 1a b S    . .................................................................................. 79 
Bảng 4.2. Các mô men bậc hai của đáp ứng của 1 2,x x theo b với 
1 2 1 2 1 2 0 1a S    .......................................................................... 79 
Bảng 4.3. Mô men bậc hai của đáp ứng với 1,,,, 22  fS   và  thay đổi ......... 85 
Bảng 4.4. Đáp ứng bình phương trung bình của hệ dao động Duffing có cản với 
2, , , , 1,fS   hệ số đàn hồi phi tuyến  thay đổi ............................................. 87 
Bảng 4.5. Đáp ứng bình phương trung bình của hệ dao động Duffing có cản với 
1,,,, 2 fS   hệ số cản phi tuyến  thay đổi ..................................................... 88 
1 
MỞ ĐẦU 
1. Tính cấp thiết, ý nghĩa lý luận và thực tiễn của đề tài: 
Việc tính toán, thiết kế dao động và điều khiển dao động có vai trò rất quan 
trọng nhằm duy trì hiệu năng, hiệu quả, cũng như kéo dài tuổi thọ của các công 
trình, máy móc. Hạn chế hay loại bỏ những dao động không mong muốn, hoặc tạo 
ra một mức độ dao động mong muốn là những mục tiêu chung của kỹ thuật dao 
động. Để nghiên cứu các hiện tượng dao động, cần phải mô tả chúng bằng các 
phương trình toán học, theo đó, các hệ dao động được gọi là tuyến tính hay phi 
tuyến tùy theo phương trình vi phân mô tả dao động ấy là tuyến tính hay phi tuyến. 
Hầu hết các bài toán dao động đều liên quan đến việc xác định chuyển động 
của một hệ gây ra bởi tải trọng kích động. Nếu tải trọng thay đổi điều hòa hay chu 
kỳ, hoặc theo một cách nào đó mà được mô tả đầy đủ như là một hàm của thời gian 
(và vị trí) và nếu vị trí hay chuyển động ban đầu của hệ được biết thì đáp ứng của hệ 
ở một thời điểm tùy ý sẽ hoàn toàn được xác định. Loại dao động như vậy được gọi 
là dao động tiền định. 
Trong dao động ngẫu nhiên, kích động thay đổi có tính chất ngẫu nhiên (bất 
thường) và gây ra các đáp ứng cũng thay đổi bất thường làm cho các công trình, 
thiết bị máy móc làm việc không hiệu quả, chóng hư hỏng, hoặc bị phá hủy đột 
ngột. Chẳng hạn như: nhiễu động của gió lên máy bay, tên lửa; tải trọng gió lên 
công trình cao tầng, tháp, cầu; tác động của sóng biển lên tầu biển, công trình biển; 
lực do nhấp nhô của mặt đường tác động lên ô tô; tác động của động đất đến các 
công trình, v.v. là các dạng kích động ngẫu nhiên. Về bản chất, giá trị của đại lượng 
ngẫu nhiên ở một thời điểm bất kỳ là không thể dự báo được, thậm chí nếu có mối 
liên hệ nào đó giữa cường độ của đại lượng ngẫu nhiên và thời gian mà đo được 
trong một khoảng thời gian nhất định thì nó cũng không bao giờ lặp lại một cách 
chính xác ở khoảng thời gian khác. Do đó, các công cụ quen thuộc trong phân tích 
dao động tiền định sẽ không thể áp dụng được cho dao động ngẫu nhiên. Bởi vậy, 
việc nghiên cứu dao động ngẫu nhiên bằng các phương pháp của cơ học ngẫu nhiên 
có ý nghĩa khoa học và ứng dụng quan trọng trong kỹ thuật. 
Trong hơn nửa thế kỷ qua, đã có những phát triển quan trọng về cơ học ngẫu 
nhiên và dao động ngẫu nhiên, đặc biệt cùng với sự phát triển và đóng góp quan 
2 
trọng của kỹ thuật máy tính. Việc xây dựng mô hình cho các bài toán nêu trên 
thường dẫn tới việc thiết lập và giải các phương trình vi phân ngẫu nhiên phi tuyến. 
Đối với dao động ngẫu nhiên mà tính chất động lực của cấu trúc được mô hình hóa 
bằng các phương trình tuyến tính thì các thống kê đáp ứng có thể được phân tích 
một cách thỏa đáng. Tuy nhiên, trong thực tế, những trường hợp như vậy chỉ là 
thiểu số, phổ biến hơn cả là các cấu trúc phi tuyến. Việc phân tích dao động ngẫu 
nhiên phi tuyến là rất khó khăn, phức tạp và chỉ một số trường hợp đặc biệt cho 
phép nhận được lời giải chính xác. Bởi vậy, cần xây dựng những phương pháp phân 
tích xấp xỉ dễ ứng dụng hơn và cho phép nhận được kết ... . Journal of Ship Research 29, 
112-126 (1985) 
[57] David C.P., James L.B., Costas P. A new stationary PDF approximation for 
nonlinear oscillators. Int. J. Nonlinear Mech. 35, 657–673 (2000) 
[58] Krylov N. M., Bogolyubov N. N.: Introduction to non-linear mechanics (in 
Russian). Kiev: Publisher AN SSSR (1937). 
[59] Atalic, T.S., Utku, S.: Stochastic linearization of multi-degree of freedom 
nonlinear systems. Earthq. Eng. Struct. Dyn.4,411–420 (1976). 
98 
[60] Spanos, P.D.: Formulation of stochastic linearization for symmetric or 
asymmetric MDOF nonlinear systems. J. App. Mech47(1), 209–211 (1980) 
[61] Faravelli, L., Casciati, F., Singh, M.P.: Stochastic equivalent linearization 
algorithms and their applicability to hysteretic systems. Meccanica 23, 107–112 
(1988). 
[62] Casciati, F., Faravelli, L.: Stochastic equivalent linearization for 3-D frames. J. 
Eng. Mech.114, 1760–1771 (1988). 
[63] Casciati, F., Faravelli, L., Venini, P.: Frequency analysis in stochastic 
linearization. J. Eng. Mech.120, 2498–2518 (1994). 
[64] Di Paola, M., Loppolo, M., Muscolino, G.: Stochastic seismic analysis of 
multi-degree of freedom systems. Eng. Struct. Elsevier6(2), 113–118 (1984). 
[65] Falsone, G.: Stochastic linearization of MDOF systems under parametric 
excitations. Int. J. Nonlinear Mech.27(6), 1025–1037 (1992). 
[66] Bellizzl, S., Bouc, R.: Analysis of multi-degree of freedom strongly nonlinear 
mechanical systems with random input. Probab. Eng. Mech. Elsevier 14(3), 229–
244 (1999). 
[67] Anh, N.D., Zakovorotny, V.L., Hieu, N.N., Diep, D.V.: A dual criterion of 
stochastic linearization method for multi-degreeof-freedom systems subjected to 
random excitation. Acta Mech.223, 2667–2684. Doi:10.1007/s00707-012-0738-5. 
[68] Anh ND, Hieu NN and Linh NN.: A dual criterion of equivalent linearization 
method for nonlinear systems subjected to random excitation. Acta Mechanica 
(2012); 223(3): 645–654. 
[69] Papoulis A. (1984), Probability, Random Variables and Stochastic Process, 
McGraw-Hill, NewYork, 2nd Edition. 
[70] Lin,Y.K., Cai,G.Q.: Probabilistic Structural Dynamics: Advanced Theory and 
Applications. McGraw-Hill, NewYork (1995) 
[71] Dimentberg, M.F.: Oscillations of a system with nonlinear cubic characteristic 
under narrow ban random excitation, Mechaics of Solids 6(2), p.142-146 (1971) 
99 
[72] Richard, K. and Anand, G.V.: Nonlinear resonance in strings under narrow 
band random excitation, Part I: Planar response and stability, Journal of Sound and 
Vibration, 86, p. 85-98 (1983) 
[73] Davies, H.G. and Nandlall, D.: Phase plane for narrow band random excitation 
of a Duffing oscillator, Journal of Sound and Vibration, 104, p. 277-283 (1986) 
[74] Iyengar, R.N.: Response of nonlinear systems to narrow band excitation, 
Structural Safety, Vol. 6, Issues. 2-4, p. 177-185 (1989) 
[75] Zhu, W.Q., Huang, C.D., Soong, T.T.: Narrow band excitation of hysteretic 
systems, Sock and Vibration, Vol. 4, N. 4, p.241-250 (1997) 
[76] Hai-Wu, R., Xiang-Dong, W., Guang, M., et al.: Response of nonlinear 
oscillators under narrow band random excitation, Appl Math Mech, Vol. 24, Issue. 
7, p. 817-825 (2003) 
[77] Cho, W.S.To, Nonlinear random vibration, CRC Press (2012) 
[78] Lutes, L.D., Sarkani, S.: Random vibration: Analysis of structural and 
mechanical systems, Elsevier, Amsterdam (2004). 
[79] Rubinstein, R. Y.: Simulation and Monte Carlo method, John Wiley and Sons, 
Inc., New York (1981) 
[80] Wojtkiewicz S.F., Johnson E.A., Bergman L. A., Spencer B.F.Jr.(1997), 
“Finite element and finite difference solutions to the transient Fokker–Planck 
equation”, DESY 161: 290–306 
[81] Zhao L., Chen Q.(1997), “An equivalent nonlinearization method for 
99nalyzing response of nonlinear systems to random excitations”, Applied 
Mathematics and Mechanics, June , Volume 18, Issue 6, pp 551-561. 
[82] Zhu W, Cai G.(2002), “Nonlinear stochastic dynamics: A survey of recent 
developments”, Acta Mechanica Sinica (English Series), Vo1.18, No.6, December. 
100 
PHỤ LỤC 
1. Các công thức tính mô men theo tiêu chuẩn LOMSEC 
Đối với quá trình chuẩn vô hướng y có trung bình bằng không, tất cả các mô men 
bậc cao ny 2 có thể được biểu diễn theo mô men bậc hai: 
nn yny 22 12...5.3.1 
Tương tự, tất cả các mô men bậc cao  ny2 trong LOMSEC cũng có thể được thể 
hiện dưới dạng các mô men bậc hai 2y
theo công thức có thể dễ dàng chứng minh 
sau khi thay thế biến 
yty  
   n
yn
y
y
n yTy y
y
2
,
2
0
0
0 2 


, n=1,2, 
trong đó 
nn
y y
22  ; 
0
0
0
2
,
y
n
yn
dttntT ; 2
2
2
1
t
etn
Gán các giá trị cụ thể cho 0, yn , ta sẽ thu được 0,ynT 
là một giá trị dương. Ngoài ra, 
tất cả các mô men bậc lẻ đều bằng 0. 
Giả sử rằng x và x là các quá trình chuẩn có trung bình bằng không, ký hiệu [.] là 
giá trị trung bình địa phương của các biến ngẫu nhiên được lấy như sau 
  xdxdxxP
x
x
x
x



0
0
0
0
,..
 (a.1) 
trong đó 00, xx  là các giá trị dương đã cho; xxP , là hàm mật độ xác suất chung 
của chúng, có thể được phân thành hai hàm mật độ xác suất đơn độc lập: 
.
2
1
)(,
2
1
)(),()(),(
2222 22 xx x
x
x
x
exPexPxPxPxxP 

 
  
 (a.2) 
trong đó x và x là độ lệch chuẩn của các biến ngẫu nhiên x và x . Các tích phân 
trong (a.1) có thể được chuyển qua các biến không thứ nguyên bằng cách 
xx rxrx   00 , với r là một giá trị dương xác định: 
101 
  xdxdxxP
x
x
x
x
r
r
r
r







,.. (a.3) 
Như đã biết, với một biến ngẫu nhiên chuẩn có trung bình zero, tất cả các mô men 
bậc lẻ đều bằng không, tất cả các mô-men bậc chẵn có bậc cao hơn bậc hai có thể 
được biểu diễn theo mô men bậc hai. Bằng cách thay thế các biến xtx  , xtx   và 
sử dụng các công thức (a.2), (a.3), khi sử dụng LOMSEC có thể được biểu diễn 
   nn xx 22 ,  theo các mô men bậc hai 22 , xx  . 
Đối với biến x: 
 
r
t
r
tnn
x
r
r
x
t
x
r
r
x
t
x
n
x
n
r
r
r
r
nn
dtedtet
dtedtet
dxxPdxxPxx
xxxx
x
x
x
x
0
2
0
222
2222
22
22
222222
2
1
2
2
1
2
2
1
2
1
)()(


 

 
 










 (a.4) 
Ký hiệu 
nn
x x
22  và 
.)(,
2
1
)(
0
2
,
22
r
n
rn
t dtttTet 

(a.5) 
Công thức (a.4) có thể viết ở dạng 
2 2 2
, 0,
2 2
, 0, 0,
0
( ) ( ) 2 2 ,
( ) 2 , ( ) 2 , ( ) .
x x
x x
x x
x x
r r
nn n
n r r
r r
r r r
nn
n r r r
r r
x x P x dx P x dx T x T
x P x dx T x P x dx T T t dt
 
 
 
 





 
 
 (a.6) 
Tương tự với biến x : 
2 2 2
, 0,
2 2
, 0, 0,
0
( ) ( ) 2 2 ,
( ) 2 , ( ) 2 , ( ) .
x x
x x
x x
x x
r r
nn n
n r r
r r
r r r
nn
n r r r
r r
x x P x dx P x dx T x T
x P x dx T x P x dx T T t dt
 
 
 
 





    
   
 (a.7) 
Cho r , các công thức (a.6) và (a.7) sẽ trở thành các công thức quen thuộc: 
102 
  nnnnnn xnxTxdxPdxxPxxx 22,222 !!122)()( 
  (a.8) 
trong đó: 
.1
2
1
2)(,2
2
1
2)(
0
22
,
0
2222 22 
dtexdxPxTdtetdxxPx t
n
n
tnn
x
n
  
  nnnnnn xnxTdxxPxdxPxxx 22,222 !!122)()(  
 (a.9) 
2
2
2 2 2 2 2
,
0
2
0
1
( ) 2 2 ,
2
1
( ) 2 1.
2
nn n n t
x n
t
x P x dx t e dt T x
P x dx e dt

   
2. Các chương trình MATLAB tính toán mô phỏng hệ dao động ngẫu 
nhiên phi tuyến 
2.1 Chương trình Matlab mô phỏng Monte Carlo cho hệ Duffing có cản chịu 
kích động ồn trắng 
% CHUONG TRINH TINH KICH DONG ON TRANG CHO HE DUFFING + can 
bac 3 = MP MONTE CARLO. 
% XDOTDOT+beta*XDOT+ga*XDOT^3+omega^2*X+epxilon*X^3=W(T) 
% Clear memory 
clear all; 
format long; 
% Pho tan so S 
S=1; 
% Buoc chia 
time=0.1; 
% POWER=Cuong do on trang mu 2=2*PI*S 
power=S*2*pi; 
103 
% c = BETA he so bo can tuyen tinh 
c=0.1; 
ga=1; 
% ga = [0.1 1 10 100] he so do can phi tuyen 
ep=100; 
% epxilon = [0.1 1 10 100] he so do cung phi tuyen 
om=0; 
Numtri=10000; 
z2=zeros(Numtri,1); 
% MOPHONG MONTE CARLO 
% So lan lap bang 10000 lan 
for i=1:Numtri 
% gieo mau ngau nhien trong (0,1) 
noise=round(100000*rand(1)); 
%goi so do simulink co ten ontrang.mdl trong thu muc C:\MATLAB7\work 
[T X Y1]=sim('ontrang', [0 300]); 
%tinh gia tri vecto dap ung chuyen dich Y1 
L=length(Y1); 
y1=Y1(L); 
%z1(i)=mean(y1); 
%plot(T,Y1) 
% Tinh gia tri binh phuong cua y1 
y2=y1*y1; 
z2(i)=y2; 
% hien thi thoi gian tinh 
104 
if mod(i,100)<1 
display('processing %') 
display(num2str(i/10000*100)) 
end; 
end; 
%Tinh gia tri binh phuong trung binh 
X2=mean(z2); 
mean_square_X=round(X2*100000)/100000 
2
Out2
1
Out1
Product1
Product
1
s
Integrator3
1
s
Integrator2
-ga
Gain6
-ep
Gain5
-c
Gain4
Band-Limited
White Noise
Add1
Hình 1. Sơ đồ Simulink Ontrang.mdl 
2.2 Chương trình Matlab tính mô men bậc 2 của hệ 2 bậc tự do chịu kích động 
ồn trắng bằng phương pháp lặp 
function mdofontrang146 
% chuong trinh dung de tinh he on trang 2 bac tu do trang 146 book of W.S.To 
% bang P.P Caughey & GLOMSEC su dung phuong phap lap {x1^2} tich phan 
trong mien tan so. 
105 
% x1'' - lamda1*x1' + alpha1*x1dot^3+omega1^2*x1 + a*x2 + b*(x1-x2)^3 = 
w1(t) 
% x2'' -(lamda1-lamda2)*x2' + alpha2*x2dot^3+omega2^2*x2 + a*x1 + b*(x2-
x1)^3 = w2(t) 
% clear all 
clear all; 
format long; 
S11=1; 
S22=1; 
Sw=[S11 0;0 S22]; 
alpha1=0.1; 
alpha2=0.1; 
omega1=1; 
omega2=1; 
b=1; 
a=1; 
lamda1=1; 
lamda2=1; 
K=[omega1^2 a;a omega2^2]; 
C=[-lamda1 0;0 -lamda1+lamda2]; 
M=[1 0;0 1]; 
GL=1; % GL = 1 tinh nghiem GLOMSEC; GL = 0 tinh nghiem Caughey 
switch GL 
 case 1 
T2rT1r=2.41189; 
T1rT0r=0.83706; 
106 
 case 0 
T2rT1r=3; 
T1rT0r=1; 
end; 
% Gia tri ban dau cua Ke, Ce 
k11=zeros(100,1); 
k12=zeros(100,1); 
k21=zeros(100,1); 
k22=zeros(100,1); 
c11=zeros(100,1); 
c12=zeros(100,1); 
c21=zeros(100,1); 
c22=zeros(100,1); 
j=1; 
k11(j)=rand(1)*1e-3; 
k22(j)=rand(1)*1e-3; 
k12(j)=-k22(j); 
k21(j)=-k11(j); 
c11(j)=rand(1)*1e-3; 
c22(j)=rand(1)*1e-3; 
c12(j)=0; 
c21(j)=0; 
j=j+1; 
k11(j)=rand(1)*1e-3; 
k22(j)=rand(1)*1e-3; 
107 
k12(j)=-k22(j); 
k21(j)=-k11(j); 
c11(j)=rand(1)*1e-3; 
c22(j)=rand(1)*1e-3; 
c12(j)=0; 
c21(j)=0; 
m=3000; 
s=1; 
dx=0.01; 
N=m/dx; 
N1=N+1; 
omega=[-m/2:dx:m/2]; 
maxerror=1e3; 
%Vong lap While tinh lap cac gia tri mean square y1 y2 y1dot y2dot 
%dk dung vong lap while 
epxilon=1e-4; 
j=j+1; 
w1=1;%doi trong cua k(j-1), c(j-1) Note: Co the bat dau bang w1=1, neu khong hoi 
tu thi giam w1 tu dong giam xuong. 
w2=1-w1;%doi trong cua k(j-2), c(j-2) Note: w1+w2=1 
% VONG LAP WHILE TINH MOMEN BAC 2 CUA DAP UNG BANG 
PHUONG PHAP TAN SO + TINH GAN DUNG BANG PHUONG PHAP LAP 
while maxerror>epxilon; 
 k11(j)=w1*k11(j-1)+w2*k11(j-2); 
 k12(j)=w1*k12(j-1)+w2*k12(j-2); 
 k21(j)=w1*k21(j-1)+w2*k21(j-2); 
108 
 k22(j)=w1*k22(j-1)+w2*k22(j-2); 
 c11(j)=w1*c11(j-1)+w2*c11(j-2); 
 c12(j)=w1*c12(j-1)+w2*c12(j-2); 
 c21(j)=w1*c21(j-1)+w2*c21(j-2); 
 c22(j)=w1*c22(j-1)+w2*c22(j-2); 
 Kej=[k11(j) k12(j);k21(j) k22(j)]; 
 Cej=[c11(j) c12(j);c21(j) c22(j)]; 
 % Tinh tich phan y1 y2 y12 y21 y1dot y2dot (Kej,Cej=Ke3,4,5,...;Ce3,4,5...) 
 s=1; 
 Sxj11=zeros(N1,1); 
 Sxj12=zeros(N1,1); 
 Sxj21=zeros(N1,1); 
 Sxj22=zeros(N1,1); 
 Sxdotj11=zeros(N1,1); 
 Sxdotj12=zeros(N1,1); 
 Sxdotj21=zeros(N1,1); 
 Sxdotj22=zeros(N1,1); 
 y1=0; 
 y1dot=0; 
 y2=0; 
 y2dot=0; 
 y21=0; 
 y12=0; 
 y12dot=0; 
 y21dot=0; 
109 
 while s<=N1; 
 Aj=zeros(2); 
 Bj=zeros(2); 
 invAj=zeros(2); 
 invBj=zeros(2); 
 Sxj=zeros(2); 
 Aj=-M*omega(s)^2+K+Kej-omega(s)*(C+Cej)*i; 
 Bj=-M*omega(s)^2+K+Kej+omega(s)*(C+Cej)*i; 
 invAj=inv(Aj); 
 invBj=inv(Bj); 
 Sxj=invAj*Sw*(invBj).'; 
 Sxj11(s)=Sxj(1,1); 
 Sxj12(s)=Sxj(1,2); 
 Sxj21(s)=Sxj(2,1); 
 Sxj22(s)=Sxj(2,2); 
 Sxdotj11(s)=Sxj(1,1)*omega(s)^2; 
 Sxdotj12(s)=Sxj(1,2)*omega(s)^2; 
 Sxdotj21(s)=Sxj(2,1)*omega(s)^2; 
 Sxdotj22(s)=Sxj(2,2)*omega(s)^2; 
 s=s+1; 
end; 
 y1=sum(Sxj11)*dx; 
 y2=sum(Sxj22)*dx; 
 y12=sum(Sxj12)*dx; 
 y21=sum(Sxj21)*dx; 
110 
 y1dot=sum(Sxdotj11)*dx; 
 y2dot=sum(Sxdotj22)*dx; 
 y12dot=sum(Sxdotj12)*dx; 
 y21dot=sum(Sxdotj21)*dx; 
 k11(j+1)=b*T2rT1r*y1+3*b*T1rT0r*y2; 
 k22(j+1)=3*b*T1rT0r*y1+b*T2rT1r*y2; 
 k12(j+1)=-k22(j+1); 
 k21(j+1)=-k11(j+1); 
 c11(j+1)=alpha1*T2rT1r*y1dot; 
 c22(j+1)=alpha2*T2rT1r*y2dot; 
 c12(j+1)=0; 
 c21(j+1)=0; 
 error1=abs(k11(j)-k11(j+1)); 
 error2=abs(k22(j)-k22(j+1)); 
 error3=abs(c11(j)-c11(j+1)); 
 error4=abs(c22(j)-c22(j+1)); 
 maxerror=max([error1,error2,error3,error4]) 
 j=j+2 
 w1=1-j/100; 
 w2=1-w1; 
end; 
switch GL 
 case 1 
y1GL=y1 
y2GL=y2 
111 
y1dotGL=y1dot 
y2dotGL=y2dot 
y12GL=y12 
y21GL=y21 
y12dotGL=y12dot 
y21dotGL=y21dot 
 case 0 
y1C=y1 
y2C=y2 
y1dotC=y1dot 
y2dotC=y2dot 
y12C=y12 
y21C=y21 
y12dotC=y12dot 
y21dotC=y21dot 
end; 
k11=b*T2rT1r*y1+3*b*T1rT0r*y2 
k22=3*b*T1rT0r*y1+b*T2rT1r*y2 
c11=alpha1*T2rT1r*y1dot 
c22=alpha2*T2rT1r*y2dot 
figure(1) 
plot(omega(1:N1),Sxdotj11(1:N1)) 
xlabel('Sxdot1'); 
figure(2) 
plot(omega(1:N1),Sxdotj22(1:N1)) 
112 
xlabel('Sxdot2'); 
figure(3) 
plot(omega(1:N1),Sxj11(1:N1)) 
xlabel('Sx1'); 
figure(4) 
plot(omega(1:N1),Sxj22(1:N1)) 
xlabel('Sx2'); 
maxSx11=max(Sxj11) 
maxSxdot11=max(Sxdotj11) 
minSx11=min(Sxj11) 
minSxdot11=min(Sxdotj11) 
b 
j 
N 
dx 
alpha1 
lamda1 
lamda2 
clear all; 
2.3 Chương trình Matlab mô phỏng Monte Carlo cho hệ Duffing chịu kích động 
ồn màu bậc 2 
% CHUONG TRINH TINH KICH DONG ON MAU BAC 2 CHO HE DUFFING 
% CAN PHI TUYEN MONTE CARLO. 
% XDOTDOT+BETA*XDOT+GAMMA*XDOT^3+EPXILON*X^3=F(T) 
% FDOTDOT+ALFA*FDOT+OMEGAF^2*F=OMEGAF^2*W(T) 
% Clear memory 
113 
clear all; 
format long; 
% Pho tan so S 
S=1; 
% Buoc chia 
time=0.1; 
% Cuong do on trang mu 2 
power=S*2*pi; 
% He so do can, do cung bo loc on mau 
afa=1; 
OmegaF=1; 
OmegaF2=OmegaF^2; 
% c = He so can = beta 
c=1; 
% Neu khong co thanh phan k*X thi lay k = 0 
k=0; 
% gamma=[0.1 1 10 100] he so do can phi tuyen 
ga=1; 
% k*ep = [0.1 1 10 100] he so do cung phi tuyen 
ep=0.1; 
% So lan lap bang 10000 lan 
for i=1:10000 
% gieo mau ngau nhien trong (0,1) 
noise=round(100000*rand(1)); 
%goi so do simulink co ten onmau2.mdl 
[T X Y1]=sim('onmau2', [0 300]); 
%tinh gia tri vecto dap ung chuyen dich Y1 
L=length(Y1); 
y1=Y1(L); 
% Tinh gia tri binh phuong cua y1 
y2=y1*y1; 
%Tinh gia tri binh phuong cua y1 
114 
z2(i)=y2; 
if mod(i,100)<1 
display('processing %') 
display(num2str(i/10000*100)) 
display('') 
display(num2str(mean(z2))) 
end; 
end; 
%Tinh gia tri binh phuong trung binh 
X2=mean(z2); 
mean_square_X=round(X2*100000)/100000 
2
Out2
1
Out1
Product1
Product
1
s
Integrator3
1
s
Integrator2
1
s
Integrator1
1
s
Integrator
-ga
Gain6
-ep
Gain5
-c
Gain4
-k
Gain3
-K-
Gain2
-K-
Gain1
-K-
Gain
Band-Limited
White Noise
Add1
Add
Hình 2. Sơ đồ Simulink onmau2.mdl