Luận án Phân tích dao động phi tuyến trong hệ chịu kích động ngẫu nhiên và tuần hoàn

 VIỆN HÀN LÂM KHOA HỌC VÀ CÔNG NGHỆ VIỆT NAM 
VIỆN CƠ HỌC 
---o0o--- 
DƯƠNG NGỌC HẢO 
PHÂN TÍCH DAO ĐỘNG PHI TUYẾN TRONG HỆ 
CHỊU KÍCH ĐỘNG NGẪU NHIÊN VÀ TUẦN HOÀN 
LUẬN ÁN TIẾN SĨ KỸ THUẬT 
Hà Nội - 2015 
ii 
VIỆN HÀN LÂM KHOA HỌC VÀ CÔNG NGHỆ VIỆT NAM 
VIỆN CƠ HỌC 
---o0o--- 
DƯƠNG NGỌC HẢO 
PHÂN TÍCH DAO ĐỘNG PHI TUYẾN TRONG HỆ 
CHỊU KÍCH ĐỘNG NGẪU NHIÊN VÀ TUẦN HOÀN 
LUẬN ÁN TIẾN SĨ KỸ THUẬT 
Chuyên ngành: Cơ kỹ thuật 
Mã số: 62 52 01 01 
Người hướng dẫn khoa học: 
GS. TSKH. Nguyễn Đông Anh 
Hà Nội - 2015 
iii 
LỜI CÁM ƠN 
Tác giả chân thành cám ơn thầy hướng dẫn GS.TSKH. Nguyễn Đông Anh đã 
tận tâm hướng dẫn khoa học, luôn động viên và giúp đỡ tác giả cả về vật chất lẫn tinh 
thần để tác giả hoàn thành luận án này. 
Tác giả xin gửi lời cám ơn đến Khoa Đào tạo sau đại học và cán bộ Viện Cơ 
học, bạn bè và đồng nghiệp tại trường đại học Công nghệ thông tin, ĐHQG Tp. HCM, 
đã động viên, giúp đỡ và tạo điều kiện thuận lợi cho tác giả trong quá trình làm luận 
án. Nhân đây, tác giả cũng gửi lời cám ơn đến NCS. Nguyễn Như Hiếu, người đã lắng 
nghe và chia sẻ rất nhiều với tác giả về chuyên môn, và đặc biệt là PGS.TS. Dương 
Anh Đức, người đã tạo điều kiện tốt nhất để tác giả an tâm thực hiện nghiên cứu của 
mình. 
Sau hết, tác giả chân thành cám ơn bố mẹ, vợ con, và gửi lời cám ơn đến người 
thân đã rất kiên nhẫn động viên tác giả trong thời gian làm luận án. 
Tác giả luận án, 
Dương Ngọc Hảo 
iv 
LỜI CAM ĐOAN 
Tôi cam đoan đây là công trình nghiên cứu của riêng tôi, dưới sự hướng dẫn 
trực tiếp của GS. TSKH. Nguyễn Đông Anh. Các số liệu, kết quả nêu trong luận án là 
trung thực và chưa từng được công bố trong bất kỳ công trình nào khác. 
Tác giả luận án, 
Dương Ngọc Hảo 
v 
MỤC LỤC 
LỜI CÁM ƠN ........................................................................................................... iii 
LỜI CAM ĐOAN .................................................................................................... iv 
MỤC LỤC ................................................................................................................... v 
DANH MỤC HÌNH VẼ, ĐỒ THỊ ............................................................................ viii 
DANH MỤC BẢNG .................................................................................................... x 
CÁC KÝ HIỆU DÙNG TRONG LUẬN ÁN ........................................................... xii 
MỞ ĐẦU ..................................................................................................................... 1 
CHƯƠNG 1.TỔNG QUAN ......................................................................................... 5 
1.1. Giới thiệu .............................................................................................................. 5 
1.2. Các phương pháp nghiên cứu hệ dao động ngẫu nhiên phi tuyến ........................... 7 
1.3. Hệ dao động chịu kích động tuần hoàn và ngẫu nhiên ......................................... 13 
1.4. Mục tiêu của luận án ........................................................................................... 15 
CHƯƠNG 2. CƠ SỞ LÝ THUYẾT ........................................................................... 16 
2.1. Các khái niệm cơ bản trong giải tích ngẫu nhiên ................................................. 16 
2.1.1. Sơ lược về lý thuyết xác suất ..................................................................... 16 
2.1.1.1. Không gian xác suất ......................................................................... 16 
2.1.1.2. Biến ngẫu nhiên ............................................................................... 17 
2.1.2. Quá trình ngẫu nhiên ................................................................................. 21 
2.1.2.1. Định nghĩa ....................................................................................... 21 
2.1.2.2. Một số quá trình ngẫu nhiên thường gặp .......................................... 22 
2.1.3. Tích phân ngẫu nhiên ................................................................................. 26 
2.1.3.1. Mở đầu ............................................................................................. 26 
vi 
2.1.3.2. Tích phân Ito – Tích phân Stratonovich ............................................ 28 
2.1.3.3. Tính chất của tích phân Ito ............................................................... 29 
2.1.4. Phương trình vi phân ngẫu nhiên ............................................................... 31 
2.2. Cơ sở lý thuyết nghiên cứu hệ dao động ngẫu nhiên ............................................ 34 
2.2.1. Phương pháp trung bình ngẫu nhiên theo biên độ và pha ........................... 34 
2.2.2. Phương pháp trung bình ngẫu nhiên trong hệ tọa độ Đề-các ...................... 36 
2.2.3. Phương pháp hàm bổ trợ và lời giải phương trình Fokker-Planck (FP) ...... 39 
2.2.3.1. Phương pháp hàm bổ trợ .................................................................. 39 
2.2.3.2. Nghiệm của phương trình FP với các hệ số dịch chuyển tuyến tính .. 40 
2.2.3.3. Tuyến tính hóa tương đương- giải xấp xỉ phương trình FP ............... 46 
2.2.4. Phương pháp mô phỏng số ......................................................................... 50 
2.3. Kết luận chương 2 ............................................................................................... 52 
CHƯƠNG 3. PHÂN TÍCH DAO ĐỘNG TRONG HỆ PHI TUYẾN CHỊU KÍCH 
ĐỘNG NGẪU NHIÊN VÀ TUẦN HOÀN ................................................................ 53 
3.1. Hệ dao động Van der Pol ..................................................................................... 55 
3.1.1. Tính toán lý thuyết ..................................................................................... 56 
3.1.2. Kết quả và thảo luận .................................................................................. 58 
3.1.3. So sánh với phương pháp phi tuyến tương đương ...................................... 65 
3.2. Hệ dao động Duffing ........................................................................................... 67 
3.2.1. Tính toán lý thuyết ..................................................................................... 67 
3.2.2. Kết quả và thảo luận .................................................................................. 69 
3.3. Dao động Van der Pol – Duffing ......................................................................... 74 
3.3.1. Tính toán lý thuyết ..................................................................................... 74 
3.3.2. Kết quả và thảo luận .................................................................................. 75 
3.4. Hệ dao động Mathieu-Duffing ............................................................................. 79 
vii 
3.4.1. Tính toán lý thuyết ..................................................................................... 79 
3.4.2. Kết quả và thảo luận .................................................................................. 82 
3.5. Kết luận chương 3 ............................................................................................... 87 
CHƯƠNG 4. PHÂN TÍCH BAN ĐẦU ĐÁP ỨNG THỨ ĐIỀU HÒA TRONG HỆ 
DAO ĐỘNG PHI TUYẾN CHỊU KÍCH ĐỘNG NGẪU NHIÊN VÀ TUẦN HOÀN 89 
4.1. Giới thiệu ............................................................................................................ 89 
4.2. Kỹ thuật phân tích ............................................................................................... 90 
4.3. Kết quả và thảo luận ............................................................................................ 97 
4.4. Kết luận chương 4 ............................................................................................. 100 
KẾT LUẬN ............................................................................................................ 102 
DANH SÁCH CÔNG TRÌNH CỦA TÁC GIẢ ĐÃ ĐƯỢC CÔNG BỐ LIÊN QUAN 
ĐẾN LUẬN ÁN ..................................................................................................... 105 
TÀI LIỆU THAM KHẢO ........................................................................................ 106 
PHỤ LỤC ................................................................................................................ 112 
Phụ lục A ................................................................................................................. 112 
Phụ lục B ................................................................................................................. 116 
viii 
DANH MỤC HÌNH VẼ, ĐỒ THỊ 
Hình 1.1. Hệ một bậc tự do. a) Kết cấu toà nhà 1 tầng. b). Mô hình tương đương. ....... 8 
Hình 2.1. Một quĩ đạo của chuyển động Brown (quá trình Wiener) ........................... 23 
Hình 2.2. Quĩ đạo của phương trình vi phân thường .................................................. 27 
Hình 2.3. Quĩ đạo của một quá trình ngẫu nhiên ........................................................ 27 
Hình 3.1.1. Đồ thị trung bình theo thời gian của trung bình bình phương đáp ứng 
theo tham số Q ........................................................................................................... 61 
Hình 3.1.2. Đồ thị trung bình theo thời gian của trung bình bình phương đáp ứng 
theo tham số Q so sánh với kết quả mô phỏng số ...................................................... 62 
Hình 3.1.3. Đồ thị hàm mật độ xác suất đồng thời ( ),p x x& của hệ dao động Van 
der Pol tại thời điểm 294t s= .................................................................................... 63 
Hình 3.1.4. Đồ thị của hàm mật độ xác suất của dịch chuyển x theo các thời gian 
khác nhau .................................................................................................................. 64 
Hình 3.1.5. Đồ thị của hàm mật độ xác suất của dịch chuyển x tại thời điểm 
294t = (s) ................................................................................................................... 64 
Hình 3.1.6. Đồ thị đường cong ( )2E x của hệ Van der Pol theo n trong lân cận 
w . .............................................................................................................................. 65 
Hình 3.2.1. Kết quả tính toán ( )E x té ùë û và ( )2E x té ùë û bằng phương pháp giải tích 
và so với kết quả mô phỏng số ................................................................................... 71 
Hình 3.2.2. Đồ thị bình phương biên độ của đáp ứng trung bình theo tham số Q ....... 71 
Hình 3.2.3. Đồ thị bình phương biên độ của đáp ứng trung bình theo tham số 2s ..... 72 
ix 
Hình 3.2.4. Đồ thị trung bình theo thời gian của trung bình bình phương đáp ứng 
( )2E x theo tham số 2s .......................................................................................... 72 
Hình 3.2.5. Đồ thị đường cong cộng hưởng của hệ Duffing. ...................................... 73 
Hình 3.3.1. Đồ thị trung bình theo thời gian của ( )2E x té ùë û theo tham số phi tuyến 
g ............................................................................................................................... 78 
Hình 3.3.2. Đồ thị trung bình theo thời gian của ( )2E x té ùë û theo biên độ lực kích 
động tuần hoàn Q ...................................................................................................... 78 
Hình 3.4.1. Kết quả giải tích ( )E x té ùë û được so sánh với các kết quả số ..................... 84 
Hình 3.4.2. Kết quả giải tích ( )2E x té ùë û được so sánh với các kết quả số ................... 84 
Hình 3.4.3. Đồ thị hàm mật độ xác suất đồng thời của hệ Mathieu-Dufing tại thời 
điểm 294t s= ........................................................................................................... 85 
Hình 3.4.4. Đồ thị hàm mật độ xác suất của x tại thời điểm 294t = (s) ..................... 86 
Hình 3.4.5. Đồ thị hàm mật độ xác suất của x tại vài thời điểm (s) ........................... 86 
Hình 4.1. Đồ thị trung bình theo thời gian của trung bình bình phương đáp ứng thứ 
điều hòa theo tham số 2s ........................................................................................... 99 
Hình 4.2. Ảnh hưởng của 2s và 0Q lên trung bình bình phương đáp ứng thứ điều 
hòa ............................................................................................................................. 99 
Hình 4.3. Ảnh hưởng 2s và h lên trung bình bình phương đáp ứng thứ điều hòa ... 100 
x 
DANH MỤC BẢNG 
Bảng 3.1.1. Sai số giữa kết quả mô phỏng và các giá trị xấp xỉ của trung bình theo 
thời gian của trung bình bình phương đáp ứng ( )2E x té ùë û theo tham số e ................. 58 
Bảng 3.1.2. Sai số giữa kết quả mô phỏng và các giá trị xấp xỉ của trung bình theo 
thời gian của trung bình bình phương đáp ứng ( )2E x té ùë û theo tham số n ................. 59 
Bảng 3.1.3. Sai số giữa kết quả mô phỏng và các giá trị xấp xỉ của trung bình theo 
thời gian của trung bình bình phương đáp ứng ( )2E x té ùë û theo tham số Q ................. 60 
Bảng 3.1.4. Sai số giữa kết quả mô phỏng và các giá trị xấp xỉ của trung bình theo 
thời gian của trung bình bình phương đáp ứng ( )2E x té ùë û theo tham số 
2s ............... 61 
Bảng 3.1.5. Sai số giữa kết quả mô phỏng và các giá trị xấp xỉ của trung bình theo 
thời gian của trung bình bình phương đáp ứng ( )2E x té ùë û theo kỹ thuật của luận án 
và phương pháp phi tuyến tương đương theo tham số ................................................ 66 
Bảng 3.2.1. Sai số giữa kết quả mô phỏng và các giá trị xấp xỉ của trung bình theo 
thời gian của trung bình bình phương đáp ứng ( )2E x té ùë û theo tham số g .................. 69 
Bảng 3.2.2. Sai số giữa kết quả mô phỏng và các giá trị xấp xỉ của trung bình theo 
thời gian của trung bình bình phương đáp ứng ( )2E x té ùë û theo tham số 
2s ................ 69 
Bảng 3.2.3. Sai số giữa kết quả mô phỏng và các giá trị xấp xỉ của trung bình theo 
thời gian của trung bình bình phương đáp ứng ( )2E x té ùë û theo tham số 
2s với các 
giá trị e khác nhau .................................................................................................... 70 
Bảng 3.3.1. Sai số giữa kết quả mô phỏng và các giá trị xấp xỉ của trung bình theo 
thời gian của trung bình bình phương đáp ứng ( )2E x té ùë û theo tham số e .................. 75 
xi 
Bảng 3.3.2. Sai số giữa kết quả mô phỏng và các giá trị xấp xỉ của trung bình theo 
thời gian của trung bình bình phương đáp ứng ( )2E x té ùë û theo tham số n .................. 76 
Bảng 3.3.3. Sai số giữa kết quả mô phỏng và các giá trị xấp xỉ của trung bình theo 
thời gian của trung bình bình phương đáp ứng ( )2E x té ùë û theo tham số g .................. 76 
Bảng 3.3.4. Sai số giữa kết quả mô phỏng và các giá trị xấp xỉ của trung bình theo 
thời gian của trung bình bình phương đáp ứn ... ntroduction with 
Application, Springer. 
[65] Ramakrishnan V, Brian FF (2012), “Resonances of a forced Mathieu equation 
with reference to wind turbine blades”, J. Vib. Acoust., 134(6). 
[66] Rayleigh JWS (1877), The Theory of Sound, reprinted by Dover, New York 
1945. 
[67] Roberts J B (1986), “First passage probabilities for randomly excited systems: 
Diffusion methods”, Probab. Eng. Mech. 1, pp. 66-81 
[68] Roberts JB, Spanos PD (1999), Random Vibration and Statistical 
Linearization, Dover Publications, Inc., Mineola, New York. 
[69] Roberts JB, Spanos PD (1986), “Stochastic averaging: An approximate 
method of solving random vibration problems”, Int. J. Nonlinear mechanics; 
21(2), pp. 111-134. 
[70] Ruby L (1996), “Applications of the Mathieu equation”, Am. J. Phys., Vol. 64, 
No. 1, pp. 39-44. 
[71] Socha L, Soong TT (1991), “Linearization in analysis of nonlinear stochastic 
systems”, Appl. Mech. Rev., 44, pp. 399-422. 
[72] Socha L (1998), “Probability density equivalent linearization technique for 
nonlinear oscillator with stochastic excitations”, Z. Angew. Math. Mech., 78, 
pp. 1087-1088. 
[73] Socha L (2008), Linearization Methods for Stochastic Dynamic System, 
Lecture Notes in Physics. Springer, Berlin. 
[74] Spanos P (1981), “Monte Carlo simulations of response of nonsymmetric 
111 
dynamic system to random excitations”, Comput. Struct. 13, pp.371-376. 
[75] Spanos P, Lutes LD (1987), “A primer of random vibration techniques in 
structural Engineering”, Shock Vib. Dig., 19(4), pp. 3-9. 
[76] Stratonovich RL(1963), Topics in the Theory of Random Noise. Vol. I, II 
(1967), New York: Gordon and Breach. 
[77] Von Wagner U, Wedig WV (2000), “On the calculation of stationary solution 
of multi-dimensional Fokker-Planck equations by orthogonal functions”, 
Nonlinear Dynamics, 21, pp. 289-306. 
[78] Xie WX, Xu W, Cai L (2006), “Study of the Duffing-Rayleigh oscillator 
subject to harmonic and stochastic excitations by path integration”, Applied 
Mathematics and Computation, 172, pp. 1212-1224. 
[79] Yu JS, Lin YK (2004), “Numerical path integration of a nonlinear oscillator 
subject to both sinusoidal and white noise excitations”, Int. J. Non-Linear 
Mechanics, 37, pp. 1493-1500. 
[80] Zhu WQ, Yu JS (1987), “On the response of the Van der Pol Oscillator to 
white noise excitation”, J. Sound and Vibration, 117(3) 421-431. 
[81] Zhu WQ (1988), “Stochastic averaging methods in random vibrations”, Appl. 
Mech. Rev. 41, pp. 189-199. 
[82] Zhu WQ, Huang ZL, Suzuki Y (2001), “Response and stability of strongly 
non-linear oscillators under wide-band random excitation”, Non-Linear 
Mechanics, 36, pp. 1235-1250. 
[83] Zhu WQ, Wu YJ (2003), “First passage time of Duffing oscillator under 
combined harmonic and white noise excitations”, Nonlinear Dynamics, 32, pp. 
291-305. 
Trang web và phần mềm: 
[84] John MC (2010), Probability Density Functions , www.mne.psu.edu/me345/ 
Lectures/Probability_density_functions.pdf . 
[85] Laurence CE (2002), An introduction to stochastic differential equations 
(version 1.2), Department of Mathematics, UC Berkeley (math.berkeley.edu/ 
~evans/SDE.course.pdf). 
[86] Jonathan MB, Matthew PS (2011), An Introduction to Modern Mathematical 
Computing With Maple™, Springer. 
[87] Jaan Kiusalaas (2010), Numerical methods in engineering with Matlab 
(Second Edition), Cambridge University Press. 
112 
PHỤ LỤC 
Phụ lục A. Xây dựng và giải hệ phương trình phi tuyến cho các hệ 
số tuyến tính hoá 
Chương trình Maple tính các hệ số tuyến tính hoá theo các mô men của 1a và 2a (để 
tránh nhiều chỉ số, trong các chương trình dưới đây luận án dùng ký hiệu b và d thay 
cho 1a và 2a ). 
----------------------------------------------- 
Đoạn chương trình tính các mô men bậc cao của 1a và 2a theo trung bình, phương sai và 
hiệp phương sai của 1a và 2a . 
113 
Tính các hệ số của hàm mật độ dừng 
----------------------------------------------- 
Chương trình Matlab tính các hệ số tuyến tính hoá 
function tuyentinhhoa_vanderpol 
% chuong trinh tim cac he so tuyen tinh hoa khi phan tich he Van der 
Pol 
global a B sig2 P nu Delta omega epsilon; 
% clear all 
clc 
a = 1; % alpha 
B = 4; % beta 
epsilon = 0.2; 
omega = 1; 
nu = 1.01; 
Delta = (omega^2-nu^2)/epsilon; 
val=[1] % co the dua nhieu gia tri vao day de co day tinh 
num=length(val); 
sig2=1; 
X2=zeros(num,1); 
for m=1:1:num 
P=val(m); 
x0=[-2,-0.5,2,1.5,-1,2]; 
L1b=x0(1); L1d=x0(2); L10=x0(3); % he so eta11, 12, 13 
L2b=x0(4); L2d=x0(5); L20=x0(6); % he so eta21, 22, 23 
alpha1 = (1/2)*a+L1b; 
114 
beta1 = Delta/(2*nu)+L1d; 
lambda1 = L10; 
alpha2 = -Delta/(2*nu)+L2b; 
beta2 = (1/2)*a+L2d; 
lambda2 = P/(2*nu)+L20; 
Ab2= -(2*(alpha1^2+alpha1*beta2+alpha2^2-... 
alpha2*beta1))*nu^2*(alpha1+beta2)/(sig2*(alpha2^2- 
2*alpha2*beta1+beta1^2+alpha1^2+2*alpha1*beta2+beta2^2)); 
 % he so tau1 
Ab1 = ((2*(2*lambda1*alpha1+2*lambda1*beta2+2*alpha2*  
lambda2-2*beta1*lambda2))*nu^2*(alpha1+beta2)/  
(sig2*(alpha2^2- 2*alpha2*beta1+beta1^2+alpha1^2+  
2*alpha1* beta2+beta2^2))); % he so tau4 
Abd = ((2*(2*alpha2*beta2+2*beta1*alpha1))*nu^2*  
(alpha1+beta2)/(sig2*(alpha2^2-2*alpha2*b  
eta1+beta1^2+alpha1^2+2*alpha1*beta2+beta2^2))); 
 % he so tau3 
Ad2 = (-(2*(alpha1*beta2+beta2^2-alpha2*beta1+beta1^2))*  
nu^2*(alpha1+beta2)/(sig2*(alpha2^2-2*alpha2*  
beta1+beta1^2+alpha1^2+2*alpha1*beta2+beta2^2))); 
% he so tau2 
Ad1 = ((2*(-2*alpha2*lambda1+2*lambda2*alpha1+  
2*lambda2*beta2+2*beta1*lambda1))*nu^2* 
(alpha1+beta2)/(sig2*(alpha2^2-2*alpha2*beta1+  
beta1^2+alpha1^2+2*alpha1*beta2+beta2^2))); % he so tau5 
psb = 2*Ad2/(4*Ab2*Ad2-Abd^2); % phuong sai cua b 
if ((Ab2<0)|| (Ad2<0) || (psb<0)) 
fprintf('\n Xin cho lai gia tri khoi tao \n') 
else 
 fprintf('\n Dang tinh..... \n') 
 r=fsolve(@myfunc, x0); 
 val(m) 
L1b=r(1); 
L1d=r(2); 
L10=r(3); 
L2b=r(4); 
L2d=r(5); 
L20=r(6); 
% cac he so cua ham dich chuyen tuyen tinh 
alpha1 = (1/2)*a+L1b; 
beta1 = Delta/(2*nu)+L1d; 
lambda1 = L10; 
alpha2 = -Delta/(2*nu)+L2b; 
beta2 = (1/2)*a+L2d; 
lambda2 = P/(2*nu)+L20; 
%------------------------------------------- 
Ab2 = -(2*(alpha1^2+alpha1*beta2+alpha2^2- 
 alpha2*beta1))*nu^2*(alpha1+beta2)/(sig2*  
(alpha2^2-2*alpha2*beta1+beta1^2+alpha1^2+ 
2*alpha1*beta2+beta2^2)); 
Ab1 = ((2*(2*lambda1*alpha1+2*lambda1*beta2+  
2*alpha2*lambda2-2*beta1*lambda2))*nu^2*  
(alpha1+beta2)/(sig2*(alpha2^2-2*alpha2*beta1+  
beta1^2+alpha1^2+2*alpha1*beta2+beta2^2))); 
Abd=(2*(2*alpha2*beta2+2*beta1*alpha1))*nu^2*  
115 
(alpha1+beta2)/(sig2*(alpha2^2-2*alpha2*beta1+  
beta1^2+alpha1^2+2*alpha1*beta2+beta2^2)); 
Ad2 = (-(2*(alpha1*beta2+beta2^2-alpha2*beta1+beta1^2))  
*nu^2*(alpha1+beta2)/(sig2*(alpha2^2-  
2*alpha2*beta1+beta1^2+alpha1^2+  
2*alpha1*beta2+beta2^2))); 
Ad1 = ((2*(-2*alpha2*lambda1+2*lambda2*  
alpha1+2*lambda2*beta2+2*beta1*lambda1))*  
nu^2*(alpha1+beta2)/(sig2*(alpha2^2-2*alpha2* ... 
beta1+beta1^2+alpha1^2+2*alpha1*beta2+beta2^2))); 
hs=[-Ab2 Ab1 Abd -Ad2 +Ad1] 
x2 = (1/2)*(2*Ab1*Ad2+Abd*Ad1)^2/(4*Ab2*Ad2-
Abd^2)^2+Ad2/(4*Ab2*Ad2-
Abd^2)+(1/2)*(2*Ab2*Ad1+Abd*Ab1)^2/(4*Ab2*Ad2-
Abd^2)^2+Ab2/(4*Ab2*Ad2-Abd^2); 
 end 
X2(m,1)=x2 
 end 
 X2 
% ghi du lieu ra file 
% fpp=fopen('data_PPTB.txt','w+'); 
% fprintf(fpp,'epsilon = %f \n',epsilon); 
% %fprintf(fpp,'P=%f \n',sig2); 
% fprintf(fpp,'P=%f \n',P); 
% fprintf(fpp,'nu=%f \n',nu); 
% 
% fprintf(fpp,'\n ================== \n \n'); 
% fprintf(fpp,'sig^2 \t \t var(Z) \n'); 
% 
% for i=1:num 
% fprintf(fpp,'%2.5f \t %2.5f \n',val(i),X2_d(i)); 
% end 
% fclose(fpp); 
function y=myfunc(x) 
% xay dung he phi tuyen cho cac he so tuyen tinh hoa 
global a B sig2 P nu Delta omega epsilon; 
Delta = (omega^2-nu^2)/epsilon; 
L1b=x(1); L1d=x(2); L10=x(3); 
L2b=x(4); L2d=x(5); L20=x(6); 
alpha1 = (1/2)*a+L1b; 
beta1 = Delta/(2*nu)+L1d; 
lambda1 = L10; 
alpha2 = -Delta/(2*nu)+L2b; 
beta2 = (1/2)*a+L2d; 
lambda2 = P/(2*nu)+L20; 
Ab2 = -(2*(alpha1^2+alpha1*beta2+alpha2^2-  
 alpha2*beta1))*nu^2*(alpha1+beta2)/(sig2*(alpha2^2-  
 2*alpha2*beta1+beta1^2+alpha1^2+2*alpha1*beta2+beta2^2)); 
Ab1 = ((2*(2*lambda1*alpha1+2*lambda1*beta2+2*alpha2*lambda2-  
 2*beta1*lambda2))*nu^2*(alpha1+beta2)/(sig2*(alpha2^2-  
 2*alpha2*beta1+beta1^2+alpha1^2+2*alpha1*beta2+beta2^2))); 
Abd = ((2*(2*alpha2*beta2+2*beta1*alpha1))*  
nu^2*(alpha1+beta2)/(sig2*(alpha2^2-2*alpha2*  
beta1+beta1^2+alpha1^2+2*alpha1*beta2+beta2^2))); 
116 
Ad2 = (-(2*(alpha1*beta2+beta2^2-alpha2*beta1+beta1^2))* 
nu^2*(alpha1+beta2)/(sig2*(alpha2^2-2*alpha2*  
beta1+beta1^2+alpha1^2+2*alpha1*beta2+beta2^2))); 
Ad1 = ((2*(-2*alpha2*lambda1+2*lambda2*alpha1+ 
2*lambda2*beta2+2*beta1*lambda1))*nu^2*  
(alpha1+beta2)/(sig2*(alpha2^2-2*alpha2*beta1+  
beta1^2+alpha1^2+2*alpha1*beta2+beta2^2))); 
tbb = (2*Ab1*Ad2+Abd*Ad1)/(4*Ab2*Ad2-Abd^2); 
tbd = (2*Ab2*Ad1+Abd*Ab1)/(4*Ab2*Ad2-Abd^2); 
psb = 2*Ad2/(4*Ab2*Ad2-Abd^2); 
psd = 2*Ab2/(4*Ab2*Ad2-Abd^2); 
corr = Abd/sqrt(4*Ab2*Ad2); 
kbd = corr*sqrt(psb*psd); 
fL1b = -(1/8*(tbd^2+psd+3*tbb^2+3*psb))*B; 
fL1d = -(1/4*(tbb*tbd+kbd))*B; 
fL10 = (1/4*(tbb^2+tbd^2))*tbb*B; 
fL2b = -(1/4*(tbb*tbd+kbd))*B; 
fL2d = -(1/8*(tbb^2+3*tbd^2+3*psd+psb))*B; 
fL20 = (1/4*(tbb^2+tbd^2))*tbd*B; 
%bieu thuc cua he phuong trinh phi tuyen 
y=[L1b-fL1b; L1d-fL1d; L10-fL10; L2b-fL2b; L2d-fL2d; L20-fL20] 
Phụ lục B. Mô phỏng hàm mật độ xác suất 
Về mặt toán học, hàm mật độ ( )p x được tính gần đúng như sau (John, 2010): 
( ) 2 2
i i
i
dx dxP x x x
p x
dx
æ ö- < £ +ç ÷
è ø= 
với 
2 2i i
dx dxP x x xæ ö- < £ +ç ÷
è ø
 là xác suất để biến ngẫu nhiên x nằm trong khoảng 
,
2 2i i
dx dxx xæ ù- +ç úè û
. 
----------------------------------------------- 
Chương trình mô phỏng số hàm mật độ xác suất của x tại thời điểm t=295 
sim_pdf_x.m 
% Mo phong ham mat do xac suat của x, tai t=295 
clear all; 
clc 
%---------------------------------------------------------- 
Pval=[1] 
epsilon=0.1; 
Om= 1; 
117 
P=1; 
omega2=Om*Om; 
sig2=1; 
alpha=1; 
gamma=1; 
nu=1.02; 
siz=length(Pval); 
lim=1000; 
power=sig2*(epsilon); 
N=25000; 
noisenum = round(100000*rand(1)); 
sim('Van_der_pol'); 
L=length(tout); 
T0=L-round(2*pi/(0.1*nu))-1; 
T=tout(T0:L); 
M1=zeros(siz,L-T0+1); 
M2=zeros(siz,L-T0+1); 
%-- khoi tao mang chua cac thong so de tinh xac suat--- 
%dt(1,i)= dau khoang chia 
%dt(2,i)=cuoi khoang chia 
%dt(3,i)=so phan tu mo phong nam trong khoang de tinh xac suat. 
%dt(4,i)=xac suat 
%dt(5,i)=gia tri ham mat do simulation 
be=-3; 
en=5; 
delt=0.1; % do rong cua khoang chia 
ndt=(en-be)/delt; 
dt=zeros(4,ndt); 
dt(1,1)=be; 
dt(2,1)=dt(1,1)+delt; 
for i=1:(ndt-1) 
 dt(1,i+1)=dt(1,i)+delt; 
 dt(2,i+1)=dt(2,i)+delt; 
end 
tim=4; % thoi diem lay mo phong, t=295 
%---------------------------------------------------------- 
for m=1:siz 
 power=Pval(m)*epsilon; 
 y1=0; 
 y2=0; 
 for j=1:N % so duong mau 
 noisenum = round(100000*rand(1)); 
 sim('Van_der_pol'); 
 y=yout(T0:L); 
 for i=1:ndt 
 if (dt(1,i)<y(tim)) && (y(tim)<=dt(2,i)) 
 dt(3,i)=dt(3,i)+1; 
 end 
 end 
 if y(tim)<=dt(1,1) 
 dt(3,1)=dt(3,1)+1; 
 end 
 if y(tim)>dt(2,ndt) 
 dt(3,ndt)=dt(3,ndt)+1; 
118 
 end 
 y1=y1+y; 
 y2=y2+y.^2; 
 end 
end 
%---ket qua mo phong--- 
for i=1:ndt 
 dt(4,i)=dt(3,i)/N; 
 dt(5,i)=dt(4,i)/delt; 
end; 
%---ket qua giai tich--- 
tbb2=2.8635 
tbb= -1.5863 
tbbd= -1.4293 
tbd2= 1.7889 
tbd=1.0453 
phi=nu*T-pi/2; 
x1=be:0.01:en; 
px=0.1775255494e-3*erf(170.1488179+.5164647559*x1).*... 
 exp(-1.724656388*x1.^2+5.927337189*x1+2.550596956)-... 
 0.1775255494e-3*erf(-173.3429356+.5164647559*x1).*... 
 exp(-1.724656388*x1.^2+5.927337189*x1+2.550596956); 
%---Ve do thi cac ham mat do--- 
plot(x1,px,'r'); 
hold on 
plot(dt(2,1:ndt)-delt/2,dt(5,1:ndt),'+k') 
% ghi du lieu ra file 
fpp=fopen('sim_pdfx.txt','w+'); 
fprintf(fpp,'x(t) \t pdf x \n'); 
for i=1:ndt 
 fprintf(fpp,'%2.5f \t %2.5f \n',dt(2,i)-delt/2,dt(5,i)); 
end 
fclose(fpp); 
---------------------------------------------- 
Chương trình mô phỏng số hàm mật độ xác suất của x và dx/dt tại thời điểm t=294 
% Mo phong ham mat do xac suat dong thoi cua x, va xdot tai t=294 
clear all; 
clc 
%-------------------------------------------------------- 
Pval=[1]; 
epsilon=0.1%, 10 100]; 
Om= 1; 
P=1; 
omega2=Om*Om; 
sig2=1; 
alpha=1; 
gamma=1; 
nu=1.02; 
siz=length(Pval); 
lim=1000; 
power=sig2*(epsilon); 
119 
noisenum = round(100000*rand(1)); 
sim('Van_der_pol'); 
L=length(tout); 
T0=L-round(2*pi/(0.1*nu))-1; 
T=tout(T0:L); 
%-- khoi tao mang chua cac thong so de tinh xac suat--- 
%dt(1,i)= dau khoang chia 
%dt(2,i)=cuoi khoang chia 
bex = -5; enx = 5; bexd = -5; enxd = 5; 
delt=0.2; % do rong cua khoang chia 
nx=(enx-bex)/delt; % so diem chia theo x 
nxd=(enx-bex)/delt;% so diem chia theo xdot 
dx=zeros(2,nx); % ma tran phan hoach truc x 
dxd=zeros(2,nxd); % ma tran phan hoach truc xd 
%---------------------------------------------------------- 
dx(1,1)=bex; 
dx(2,1)=bex+delt; 
for i=1:(nx-1) 
 dx(1,i+1)=dx(1,i)+delt; 
 dx(2,i+1)=dx(2,i)+delt; 
end 
%---------------------------------------------------------- 
dxd(1,1)=bexd; 
dxd(2,1)=bexd+delt; 
for i=1:(nxd-1) 
 dxd(1,i+1)=dxd(1,i)+delt; 
 dxd(2,i+1)=dxd(2,i)+delt; 
end 
%---------------------------------------------------------- 
N=40000; % so duong mau 
tim=4; % thoi diem lay mo phong 
%kq(1,i)=so phan tu mo phong nam trong khoang de tinh xac suat. 
%pr(2,i)=xac suat 
%pdf(3,i)=gia tri ham mat do simulation 
kq=zeros(nx,nxd); 
pr=zeros(nx,nxd); 
pf=zeros(nx,nxd); 
%---------------------------------------------------------- 
for m=1:siz 
 power=Pval(m)*epsilon; 
 for k=1:N % so lan tinh va lai lay trungbinh 
 noisenum = round(100000*rand(1)); 
 [t,z]=sim('Van_der_pol'); 
 y=z(T0:L,:); 
 for i=1:nx 
 for j=1:nxd 
 if (dx(1,i)<y(tim,1)) && ... 
(y(tim,1)<=dx(2,i))... 
 && (dxd(1,j)<y(tim,2)) && ... 
(y(tim,2)<=dxd(2,j)) 
 kq(i,j)=kq(i,j)+1; 
 end 
 end 
 end 
 end 
120 
end 
for i=1:nx 
 for j=1:nxd 
 pr(i,j)=kq(i,j)/N; 
 pf(i,j)=pr(i,j)/(delt*delt); 
 end 
end; 
[x1 x2]=meshgrid(-5:0.1:5); 
wxxd=0.1754576952e-3*exp(-1.84*(-.9899209965*x1-... 
 .1388438165*x2).^2-... 
 .9186*(.1416206928*x1-.9705107809*x2).^2+(1.2117*... 
 (-.9899209965*x1-.1388438165*x2)).*... 
 (.1416206928*x1-.9705107809*x2)+... 
 7.576986402*x1-2.742882717*x2); 
mesh(x1,x2,wxxd) 
mesh(dx(2,1:nx)-delt/2,dxd(2,1:nxd)-delt/2,pf') 
xlabel('displacement (x)'); 
ylabel('velocity (dx/dt)'); 
zlabel('joint PDF of x and dx/dt'); 
axis=([-5 5 -5 5 -0.1 0.6]); 
% ghi du lieu ra file 
fpp=fopen('simpdf_t_294.txt','w+'); 
fprintf(fpp,'t \t x(t) \t xdot \t pdf x \n'); 
for i=1:ndt 
 fprintf(fpp,'%2.5f \t %2.5f \t %2.5f ... 
\t %2.5f \n',dt(2,i)-delt/2,dt(5,i)); 
end 
fclose(fpp);