Luận án Phân tích dao động ngẫu nhiên phi tuyến bằng phương pháp tuyến tính hóa tương đương

VIỆN HÀN LÂM KHOA HỌC VÀ CÔNG NGHỆ VIỆT NAM 
VIỆN CƠ HỌC 
--o0o-- 
NGUYỄN NGỌC LINH 
PHÂN TÍCH DAO ĐỘNG NGẪU NHIÊN PHI TUYẾN 
BẰNG PHƯƠNG PHÁP TUYẾN TÍNH HÓA 
TƯƠNG ĐƯƠNG 
LUẬN ÁN TIẾN SĨ 
HÀ NỘI – 2015 
VIỆN HÀN LÂM KHOA HỌC VÀ CÔNG NGHỆ VIỆT NAM 
VIỆN CƠ HỌC 
--o0o-- 
NGUYỄN NGỌC LINH 
PHÂN TÍCH DAO ĐỘNG NGẪU NHIÊN PHI TUYẾN 
BẰNG PHƯƠNG PHÁP TUYẾN TÍNH HÓA 
TƯƠNG ĐƯƠNG 
Chuyên ngành: Cơ kỹ thuật 
Mã số: 62 52 01 01 
LUẬN ÁN TIẾN SĨ 
 NGƯỜI HƯỚNG DẪN KHOA HỌC 
 1. GS.TSKH Nguyễn Đông Anh 
 2. TS Lưu Xuân Hùng 
HÀ NỘI – 2015 
II
LỜI CÁM ƠN 
Tôi xin chân thành cám ơn các thầy hướng dẫn khoa học, đặc biệt là 
GS.TSKH. Nguyễn Đông Anh, đã tận tâm hướng dẫn khoa học, truyền niềm say 
mê nghiên cứu và giúp đỡ tôi hoàn thành luận án này. 
Tôi xin bày tỏ sự cám ơn tới Phòng Cơ học Công Trình, Khoa Đào tạo sau 
đại học và bạn bè, đồng nghiệp trong Viện Cơ học đã giúp đỡ tôi ngay từ những 
ngày đầu ở Viện Cơ học. 
Tôi xin bày tỏ sự cám ơn tới đơn vị công tác là Trường Cao đẳng Xây 
dựng số 1 đã ủng hộ, tạo điều kiện thuận lợi cho tôi trong quá trình làm nghiên 
cứu sinh. 
Tôi cũng xin bày tỏ sự cám ơn tới TS Lã Đức Việt, chủ nhiệm đề tài “Tối 
ưu hóa tham số các hệ tiêu tán hoặc tích trữ năng lượng trong điều khiển và 
giám sát kết cấu”, mã số 107.04-2011.14- Nafosted đã tạo điều kiện cho tôi 
được tham gia đề tài và có những hỗ trợ tài chính giúp ích cho quá trình làm 
nghiên cứu sinh. 
Cuối cùng tôi xin bày tỏ sự biết ơn đến gia đình đã động viên ủng hộ tôi 
trong thời gian làm luận án. 
III 
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 số liệu, 
kết quả nêu trong luận án là trung thực và chưa từng được ai công bố trong bất 
kỳ công trình nào khác. 
 Tác giả luận án 
 Nguyễn Ngọc Linh 
IV
MỤC LỤC 
LỜI CÁM ƠN .................................................................................................... II 
LỜI CAM ĐOAN ........................................................................................... III 
MỤC LỤC ...................................................................................................... IV 
DANH MỤC MỘT SỐ KÝ HIỆU VÀ CHỮ VIẾT TẮT ................................. VI 
DANH MỤC HÌNH VẼ, ĐỒ THỊ ................................................................ VIII 
DANH MỤC BẢNG VÀ SƠ ĐỒ KHỐI .......................................................... IX 
MỞ ĐẦU ........................................................................................................... 1 
CHƯƠNG 1. CƠ SỞ LÍ THUYẾT XÁC SUẤT VÀ MỘT SỐ PHƯƠNG PHÁP 
PHÂN TÍCH DAO ĐỘNG NGẪU NHIÊN PHI TUYẾN .................................. 4 
1.1 Các khái niệm cơ bản về xác suất ............................................................. 4 
1.1.1 Xác suất của sự kiện ngẫu nhiên ........................................................ 4 
1.1.2. Biến ngẫu nhiên và các đặc trưng xác suất ....................................... 5 
1.2 Quá trình ngẫu nhiên ................................................................................ 7 
1.2.1 Các đặc trưng của quá trình ngẫu nhiên ............................................. 8 
1.2.2 Các quá trình ngẫu nhiên đặc biệt .................................................... 12 
1.3 Phương trình Fokker-Planck-Kolgomorov (FPK) ................................... 16 
1.4 Dao động ngẫu nhiên chịu kích động ồn trắng Gauss ............................. 17 
1.4.1 Dao động ngẫu nhiên tuyến tính chịu kích động ồn trắng Gauss ..... 19 
1.4.2 Dao động ngẫu nhiên phi tuyến chịu kích động ồn trắng Gauss ...... 20 
1.5 Một số phương pháp xấp xỉ trong phân tích dao động ngẫu nhiên phi 
tuyến ............................................................................................................ 22 
1.5.1 Phương pháp nhiễu .......................................................................... 22 
1.5.2 Phương pháp trung bình hóa ngẫu nhiên .......................................... 22 
1.5.3 Phương pháp tuyến tính hóa tương đương ngẫu nhiên ..................... 23 
1.5.4 Phương pháp phi tuyến hóa tương đương ngẫu nhiên ...................... 24 
1.5.5 Phương pháp sử dụng hàm mật độ phổ ............................................ 24 
1.5.6 Phương pháp xấp xỉ cho phương trình FPK ..................................... 25 
1.5.7 Phương pháp mô phỏng số Monte Carlo .......................................... 26 
1.5.8 Nhận xét .......................................................................................... 26 
Kết luận chương 1 ............................................................................................ 28 
CHƯƠNG 2. PHƯƠNG PHÁP TUYẾN TÍNH HÓA TƯƠNG ĐƯƠNG NGẪU 
NHIÊN ............................................................................................................. 29 
2.1 Tiêu chuẩn kinh điển .............................................................................. 33 
V
2.2 Tiêu chuẩn cực tiểu sai số thế năng ........................................................ 35 
2.3 Tiêu chuẩn tuyến tính hóa tương đương có điều chỉnh............................ 35 
2.4 Tuyến tính hóa tương đương dựa trên phân bố khác Gauss ..................... 37 
2.5 Tiêu chuẩn tuyến tính hóa từng phần ...................................................... 38 
Kết luận chương 2 ............................................................................................ 39 
CHƯƠNG 3. TIÊU CHUẨN ĐỐI NGẪU CỦA PHƯƠNG PHÁP TUYẾN 
TÍNH HÓA TƯƠNG ĐƯƠNG NGẪU NHIÊN ............................................... 40 
3.1. Ý tưởng cơ bản của tiêu chuẩn đối ngẫu có trọng số tổng quát .............. 40 
3.2. Tiêu chuẩn đối ngẫu .............................................................................. 42 
3.2.1 Khái niệm về tiêu chuẩn đối ngẫu .................................................... 42 
3.2.2 Mức độ phụ thuộc tuyến tính trong tiêu chuẩn đối ngẫu .................. 44 
3.2.3 Ý nghĩa hình học của tiêu chuẩn đối ngẫu........................................ 46 
3.2.4 Áp dụng tiêu chuẩn đối ngẫu để phân tích mô men bậc hai của dao 
động ngẫu nhiên phi tuyến ........................................................................ 50 
3.3 Các ví dụ áp dụng của phương pháp tuyến tính hóa tương đương theo tiêu 
chuẩn đối ngẫu ............................................................................................. 52 
3.3.1. Dao động Van der pol ..................................................................... 52 
3.3.2 Dao động có cản phi tuyến bậc ba.................................................... 54 
3.3.3 Dao động Duffing ............................................................................ 55 
3.3.4 Dao động có cản và đàn hồi phi tuyến ............................................. 58 
3.3.5 Dao động Lutes Sarkani ................................................................... 60 
CHƯƠNG 4. TIÊU CHUẨN ĐỐI NGẪU CÓ TRỌNG SỐ CỦA PHƯƠNG 
PHÁP TUYẾN TÍNH HÓA TƯƠNG ĐƯƠNG NGẪU NHIÊN ...................... 66 
4.1 Tiêu chuẩn đối ngẫu có trọng số ............................................................. 66 
4.1.1 Khái niệm về tiêu chuẩn đối ngẫu có trọng số .................................. 66 
4.1.2 Xác định dạng giải tích của trọng số ................................................ 68 
4.1.3 Một số tính chất khác của tiêu chuẩn đối ngẫu có trọng số được đề 
xuất .......................................................................................................... 74 
4.1.4 Áp dụng tiêu chuẩn đối ngẫu có trọng số để phân tích mô men bậc hai 
của dao động ngẫu nhiên phi tuyến ........................................................... 76 
4.2 Các ví dụ áp dụng của phương pháp tuyến tính hóa tương đương theo tiêu 
chuẩn đối ngẫu có trọng số ........................................................................... 78 
4.2.1 Dao động đàn hồi phi tuyến không cản ............................................ 78 
4.2.2 Dao động đàn hồi phi tuyến theo qui luật mũ ................................... 80 
4.2.3 Dao động đàn hồi phi tuyến bậc 3 và bậc 5 ...................................... 83 
4.2.4 Dao động có cản phi tuyến phụ thuộc năng lượng ........................... 85 
VI
4.2.5 Dao động tự do ............................................................................... 88 
Kết luận chương 4 ............................................................................................ 91 
KẾT LUẬN VÀ KIẾN NGHỊ .......................................................................... 93 
DANH SÁCH CÔNG TRÌNH ĐÃ ĐƯỢC CÔNG BỐ CỦA LUẬN ÁN ......... 95 
TÀI LIỆU THAM KHẢO ................................................................................ 96 
PHỤ LỤC ...................................................................................................... 103 
VII
 DANH MỤC MỘT SỐ KÝ HIỆU VÀ CHỮ VIẾT TẮT 
A véc tơ, hàm phi tuyến 
 ,a x t véc tơ hệ số dịch chuyển 
, ,c  hệ số hằng số 
B véc tơ, hàm tuyến tính tương đương 
b, k 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 
,ttb 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 
3 5,c c hệ số độ cứng phi tuyến 
 1 2 12, ,xxD t t D hiệp phương sai 
 kd  tỉ số các hệ số tuyến tính hóa theo các tiêu chuẩn 
 đối ngẫu có trọng số và kinh điển 
 _ _,S ts S kdd d  tỉ số không thứ nguyên của minS 
 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 
FPK phương trình Fokker-Planck-Kolgomorov 
 ,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 
 hệ số trở về 
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 
VIII 
 mức độ phụ thuộc tuyến tính 
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 
p, p  trọng số, hàm trọng số 
 , ,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 
r hệ số tương quan 
r2 hệ số tương quan bình phương 
S biểu thức của sai số phương trình 
 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 
 góc giữa hai véc tơ 
 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 HÌNH VẼ, ĐỒ THỊ 
Hình 1.1 Hàm mật độ của tải trọng tác dụng lên dầm ........................................ 6 
Hình 1.2. Một tổng thể các hàm ngẫu nhiên theo thời gian (các hàm mẫu) ........ 7 
Hình 1.3 Hàm mật độ xác suất ........................................................................... 8 
Hình 1.4 Tương quan dương và tương quan âm ............................................... 11 
Hình 1.5 Hàm mật độ phổ và hàm tự tương quan của quá trình ồn trắng.......... 14 
Hình 1.6 Mô hình hệ cơ học một bậc tự do ...................................................... 17 
Hình 2.1 Các cách tiếp cận chủ yếu của các phương pháp tuyến tính hóa và phi 
tuyến hóa tương đương ngẫu nhiên .................................................................. 31 
Hình 3.1 Phép chiếu véc tơ của tiêu chuẩn đối ngẫu ........................................ 49 
Hình 3.2 Hàm thế năng dạng giếng đơn và giếng đôi ....................................... 56 
Hình 3.3 Mô hình hệ một bậc tự do chuyển động có ma sát ............................. 60 
Hình 4.1 Mức độ phụ thuộc tuyến tính mạnh nhất và yếu nhất ........................ 68 
Hình 4.2 Hàm nội suy tuyến tính p(µ) ............................................................ 72 
Hình 4.3 Đồ thị hàm tuyến tính từng đoạn p(µ) ............................................... 73 
Hình 4.4 Tỉ số kd  ....................................................................................... 75 
Hình 4.5 Các tỉ số _ _,S ts S kdd d  ................................................................. 76 
X
DANH MỤC BẢNG VÀ SƠ ĐỒ KHỐI 
Bảng 1.1 Phân loại mức độ tương quan tuyến tính của Cohen ......................... 12 
Bảng 3.1 Đáp ứng trung bình bình phương với 1, 0.5, 2o h  ;  thay 
đổi .................................................................................................................... 42 
Bảng 3.2 Đáp ứng trung bình bình phương của dao động Van der Pol với 
α=0.2;  =1; =2; σ2 thay đổi ......................................................................... 54 
Bảng 3.3 Đáp ứng trung bình bình phương của dao động có cản phi tuyến bậc 
ba với 0.05, 1, 4oh h  , và γ thay đổi ....................................................... 55 
Bảng 3.4 Đáp ứng trung bình bình phương của dao động Duffing với 1o , 
0.5h , 12, 1c ; 3c thay đổi ....................................................................... 57 
Bảng 3.5 Đáp ứng trung bình bình phương của dao động Duffing với 1o , 
0.1h , 12, 1c ; 3c thay đổi ..................................................................... 57 
Bảng 3.6 Đáp ứng trung bình bình phương của dao động có cản và đàn hồi phi 
tuyến, với 00.1, 1, 1h   và  thay đổi ..................................................... 60 
Bảng 3.7 Đáp ứng trung bình bình phương của dao động Lutes Sarkani với 
0 1; 1S  và a thay đổi ................................................................................. 63 
Bảng 3.8 Sai số của tiêu chuẩn đối ngẫu và tiêu chuẩn kinh điển theo giá trị của 
mức độ phụ thuộc tuyến tính ............................................................................ 65 
Bảng 4.1 Các giá trị chính xác của i và ( )cx ip của dao động Lutes Sarkani với 
mức độ phụ thuộc tuyến tính yếu ( 0 1/ 3 ) ................................................. 72 
Bảng 4.2 Các giá trị chính xác của j và ( )cx jp của dao động Lutes Sarkani với 
mức độ phụ thuộc tuyến tính mạnh ( 2 / 3 1 ) ............................................... 72 
Bảng 4.3 Trọng số và hệ số tuyến tính hóa của tiêu chuẩn đối ngẫu có trọng số
 ......................................................................................................................... 74 
Bảng 4.4 Đáp ứng trung bình bình phương của dao động đàn hồi phi tuyến 
không cản với 1 , 2 và  thay đổi ....................................................... 79 
Bảng 4.5 Đáp ứng trung bình bình phương của dao động đàn hồi phi tuyến theo 
qui luật mũ với 0.5h ; 0 1 ; 2 , a = 1/5 và γ thay đổi ............................ 81 
Bảng 4.6 Đáp ứng trung bình bình phương của dao động đàn hồi phi tuyến theo 
qui luật mũ với 0.5h ; 0 1 ; 2 , a = 1/3 và γ thay đổi ............................ 82 
Bảng 4.7 Đáp ứng trung bình bình phương của dao động đàn hồi phi tuyến theo 
qui luật mũ với 0.5h ; 0 1 ; 2 , a = 5 ... ------------------------------------------------ 
Na=max(size(n)); 
for i=1:1:Na 
%Nghiem chinh xac 
EBB_e(i)=fsolve(@(Ex)MyFunexact(Ex,n(i),So,gamma0),1); 
sigmax_e(i)=sqrt(EBB_e(i)); 
EAB_e(i)=gamma0*sigmax_e(i).^(n(i)+1)*2.^(n(i)/2)*n(i)*gamma(n(i)/2)/ 
sqrt(2*pi); 
EAA_e(i)=gamma0.^2*sigmax_e(i).^(2*n(i))*2.^(n(i)+1/2)*gamma(n(i)+1/)/
sqrt(2*pi); 
miu(i)=(n(i)^2*gamma(n(i)/2)^2)/(2*pi^(1/2)*gamma(n(i) + 1/2)); 
k_e(i)=pi*So/(EBB_e(i)); 
p_e(i)=(k_e(i)- EAB_e(i)/EBB_e(i))/(miu(i)*k_e(i)-EAB_e(i)/EBB_e(i)); 
%Tinh anpha1, anpha2 
anpha11(i)= miu(i)*(1- p_e(i)); 
anpha12(i)= miu(i).^2; 
anpha21(i)= miu(i)* p_e(i); 
anpha22(i)= miu(i)*(1- miu(i)); 
end 
%DAU RA ------------------------------------------------------------------- 
anpha1TS=anpha11(1,10:18)'; 
anpha1MS=anpha12(1,10:18)'; 
anpha1=-sum(anpha1TS)/sum(anpha1MS); 
109
anpha2TS=anpha11(1,1:9)'; 
anpha2MS=anpha12(1,1:9)'; 
anpha2=-sum(anpha2TS)/sum(anpha2MS); 
%VE ----------------------------------------------------------------------- 
plot(miu,p_e,'*-','LineWidth',2) 
grid on 
xlabel('miu') 
ylabel('p(miu)') 
%title('Lutes-Sarkani oscillator – The exact weighting coefficient') 
%HAM CON ------------------------------------------------------------------ 
function F=MyFunexact(Ex,n,So,gamma0) 
F=Ex-(pi.*So./gamma0).^(2./(n+1)).*(n+1).^(2./(n+1)).*gamma(3./(n+1)) 
./gamma(1./(n+1)); 
%----------------------------------END------------------------------------% 
7. Tính toán mô men bậc hai của dao động đàn hồi phi tuyến không cản 
function first_oder_nonlinear_spring 
clear; 
%clc; 
%THONG SO DAU VAO --------------------------------------------------------- 
c1=1; 
sigma=sqrt(2); 
n=3; 
gamma=[0.1 0.5 1 5 10 50 100]; 
%GIAI --------------------------------------------------------------------- 
Ngamma=max(size(gamma)); 
for i=1:1:Ngamma 
%nghiem chinh xac 
TS(i)=quad(@(z)Fun1(z,h,c1,sigma,gamma(i),n),0,10); 
MS(i)=quad(@(z)Fun2(z,h,c1,sigma,gamma(i),n),0,10); 
Exx_e(i)=TS(i)/MS(i); 
%nghiem kinh dien 
x01=x_e(i); %gia tri ban dau 
Exx_kd(i)=fsolve(@(Exx)MyFun_kd(Exx,c1,sigma,gamma(i),n),x01); 
e_kd(i)=abs(Exx_kd(i)- Exx_e(i))*100/x_e(i); %sai so 
%nghiem doi ngau 
x02=x_e(i); %gia tri ban dau 
Exx_dn=fsolve(@(Exx)MyFun_dn (Exx,c1,sigma,gamma(i),n),x02); 
e_dn(i)=abs(Exx_dn (i)-x_e(i))*100/Exx_e(i); %sai so 
end 
%DAU RA ------------------------------------------------------------------- 
result=[gamma' Exx_e' Exx_kd' e_kd' Exx_dn ' e_dn'] 
%VE ----------------------------------------------------------------------- 
plot(gamma,Exx_e,'ko') 
hold on 
plot(gamma,Exx_kd,'kx') 
hold on 
plot(gamma,Exx_dn,'k-*') 
hold on 
grid on 
xlabel(gamma') 
ylabel('Exx') 
%HAM CON ------------------------------------------------------------------ 
110
%Tieu chuan doi ngau 
function F_dn =MyFun_dn(Exx,c1,sigma,gamma,n) 
Exg=4*gamma*Exx^(3/2)/sqrt(2*pi); 
Egg=3*gamma^2*Exx^2; 
miu=Exg^2/(Exx*Egg); 
k=1/(2-miu)*Exg/Exx; 
F_dn=Exx-sigma^2/(2*(c1+k)); 
%Tieu chuan kinh dien 
function F_kd=MyFun_kd(Exx,c1,sigma,gamma,n) 
Exg=4*gamma*Exx^(3/2)/sqrt(2*pi); 
k=Exg/Exx; 
F_cv= Exx-sigma^2/(2*(c1+k)); 
%Nghiem chinh xac 
function FunTS=Fun1(z,c1,sigma,gamma,n) 
FunTS=z^2*exp(-2/(sigma^2)*(anpha*z^2/2+gamma*abs(z)^3/3)); 
function FunMS=Fun2(z,c1,sigma,gamma,n) 
FunMS=exp(-2/(sigma^2)*(anpha*z^2/2+gamma*abs(z)^3/3)); 
%----------------------------------END------------------------------------% 
8. Tính toán mô men bậc hai của dao động đàn hồi phi tuyến theo qui luật mũ 
function Nonlinear_spring 
clear; 
%clc; 
%THONG SO DAU VAO --------------------------------------------------------- 
h=0.5; 
omega=1; 
sigma=sqrt(2); 
n=1/5; %1/3; 5; 7 
gamma0=[0.1 0.5 1 5 10 50 100]; 
%GIAI --------------------------------------------------------------------- 
Ngamma0=max(size(gamma0)); 
for i=1:1:Ngamma 
%nghiem chinh xac 
 TS(i)=quad(@(z)Fun1(z,h,omega,sigma,gamma0(i),n),0,10); 
 MS(i)=quad(@(z)Fun2(z,h,omega,sigma,gamma0(i),n),0,10); 
 Exx_e(i)=TS(i)/MS(i); 
%muc do phu thuoc tuyen tinh 
 miu(i)=(n.^2*gamma(n/2).^2)/(2*pi.^(1/2)*gamma(n + 1/2)); 
 if miu(i)>=0 && miu(i)<=1 
 if miu(i)<1/3 
 p(i)=-6.*miu(i)./5+1; 
 elseif miu(i)==1/3 
 p(i)=27./49; 
 elseif miu(i)>1/3 && miu(i)<=2/3 
 p(i)=1./2; 
 elseif miu(i)>2/3 
 p(i)=-3.*miu(i)./2+3./2; 
 end 
 end 
%nghiem doi ngau trong so 
 x01=Exx_e(i); %gia tri ban dau 
Exx_ts(i)=fsolve(@(Exx)MyFun_ts(Exx,h,omega,sigma,gamma0(i),n,p(i)) 
,x01); 
 e_ts(i)=abs(Exx_ts(i)-Exx_e(i))*100/Exx_e(i); %sai so 
%nghiem doi ngau 
 x02=Exx_e(i); %gia tri ban dau 
 Exx_dn(i)=fsolve(@(Ex1)MyFun_dn(Ex1,h,omega,sigma,gamma0(i),n),x02); 
111
 e_dn(i)=abs(Exx_dn(i)-Exx_e(i))*100/Exx_e(i); %sai so 
%nghiem kinh dien 
 x03=Exx_e(i); %gia tri ban dau 
 Exx_kd(i)=fsolve(@(Ex_kd)MyFun_kd(Ex_kd,h,omega,sigma,gamma0(i),n) 
,x03); 
 e_kd(i)=abs(Exx_kd(i)-Exx_e(i))*100/Exx_e(i); 
end 
%DAU RA ------------------------------------------------------------------- 
result=[gamma0' Exx_e' Exx_kd' e_kd' Exx_dn ' e_dn' Exx_ts' e_ts' miu'] 
%VE ----------------------------------------------------------------------- 
plot(gamma0,Exx_e,'ko') 
hold on 
plot(gamma0,Exx_kd,'kx') 
hold on 
plot(gamma0,Exx_dn,'k-*') 
hold on 
plot(gamma0,Exx_ts,'k+') 
hold on 
grid on 
xlabel(gamma0') 
ylabel('Exx') 
%HAM CON ------------------------------------------------------------------ 
%Tieu chuan doi ngau trong so 
function F_ts =MyFun_ts(Exx,h,omega,sigma,gamma0,n,p) 
Exg= gamma0*sigmax^(n+1)*2^(n/2)*n*gamma(n/2)/sqrt(2*pi); 
Egg= gamma0^2*sigmax^(2*n)*2^(n+1/2)*gamma(n+1/2)/sqrt(2*pi); 
miu=Exg^2/(Exx*Egg); 
k=(1-p)/(1-miu*p)*Exg/Exx; 
F_ts=Exx-sigma^2/(2*(2*h*(omega+k)); 
%Tieu chuan doi ngau 
function F_dn =MyFun_dn(Exx,h,omega,sigma,gamma0,n) 
Exg= gamma0*sigmax^(n+1)*2^(n/2)*n*gamma(n/2)/sqrt(2*pi); 
Egg= gamma0^2*sigmax^(2*n)*2^(n+1/2)*gamma(n+1/2)/sqrt(2*pi); 
miu=Exg^2/(Exx*Egg); 
k=1/(2-miu)*Exg/(Exx); 
F_dn=Exx-sigma^2/(2*(2*h*(omega+k)); 
%Tieu chuan kinh dien 
function F_kd=MyFun_kd(Exx,h,omega,sigma,gamma0,n) 
Exg= gamma0*sigmax^(n+1)*2^(n/2)*n*gamma(n/2)/sqrt(2*pi); 
k=Exg/(Exx); 
F_kd= Exx-sigma^2/(2*(2*h*(omega+k)); 
%Nghiem chinh xac 
function FunTS=Fun1(z,h,omega,sigma,gamma0,n) 
FunTS=(z^2*exp(-4*h*(omega^2*z^2/2 + gamma0*z^(n+1)/(n+1))/(sigma^2))); 
function FunMS=Fun2(z,h,omega,sigma,gamma0,n) 
FunMS=(exp(-4*h*(omega^2*z^2/2 + gamma0*z^(n+1)/(n+1))/(sigma^2))); 
%----------------------------------END------------------------------------% 
9. Tính toán mô men bậc hai của dao động đàn hồi phi tuyến bậc 3 và bậc 5 
function Nonlinear_spring_3th_plus_5th 
clear; 
clc; 
%THONG SO DAU VAO --------------------------------------------------------- 
h=0.5; 
sigma=sqrt(2); 
n1=3; 
n2=5; 
112
c1=1; 
c3=1; 
c5=[0.5 1 5 10 50 100]; 
%GIAI --------------------------------------------------------------------- 
Nc5=max(size(c5)); 
for i=1:1:Nc5 
%nghiem chinh xac 
 TS(i)=quad(@(z)Fun1(z,h,sigma,c1,c3,c5(i)),0,5); 
 MS(i)=quad(@(z)Fun2(z,h,sigma,c1,c3,c5(i)),0,5); 
 Exx_e(i)=TS(i)/MS(i); 
%muc do phu thuoc tuyen tinh 
 miu1(i)=(n1.^2*gamma(n1/2).^2)./(2*pi.^(1/2)*gamma(n1 + 1/2)); 
 if miu1(i)>=0 && miu1(i)<=1 
 if miu1(i)<1/3 
 p1(i)=-6.*miu1(i)./5+1; 
 elseif miu1(i)==1/3 
 p1(i)=27./49; 
 elseif miu1(i)>1/3 && miu1(i)<=2/3 
 p1(i)=1./2; 
 elseif miu1(i)>2/3 
 p1(i)=-3.*miu1(i)./2+3./2; 
 end 
 end 
 miu2(i)=(n2.^2*gamma(n2/2).^2)./(2*pi.^(1/2)*gamma(n2 + 1/2)); 
 if miu2(i)>=0 && miu2(i)<=1 
 if miu2(i)<1/3 
 p2(i)=-6.*miu2(i)./5+1; 
 elseif miu2(i)==1/3 
 p2(i)=27./49; 
 elseif miu2(i)>1/3 && miu2(i)<=2/3 
 p2(i)=1./2; 
 elseif miu2(i)>2/3 
 p2(i)=-3.*miu2(i)./2+3./2; 
 end 
 end 
%nghiem doi ngau trong so 
 x01=Exx_e(i); %gia tri ban dau 
Exx_ts(i)=fsolve(@(Exx)MyFun_ts(Exx,h,sigma,c1,c3,c5(i),p1(i),p2(i)) ,x01); 
 e_ts(i)=abs(Exx_ts(i)-Exx_e(i))*100/Exx_e(i); %sai so 
%nghiem doi ngau 
 x02=Exx_e(i); %gia tri ban dau 
 Exx_dn(i)=fsolve(@(Exx)MyFun_dn(Exx,h,sigma,c1,c3,c5(i)),x02); 
 e_dn(i)=abs(Exx_dn(i)-Exx_e(i))*100/Exx_e(i); %sai so 
%nghiem kinh dien 
 x03=Exx_e(i); %gia tri ban dau 
 Exx_kd(i)=fsolve(@(Exx)MyFun_kd(Exx,h,sigma,c1,c3,c5(i)),x03); 
 e_kd(i)=abs(Exx_kd(i)-Exx_e(i))*100/Exx_e(i); 
end 
%DAU RA ------------------------------------------------------------------- 
result=[c5' Exx_e' Exx_kd' e_kd' Exx_dn' e_dn' Exx_ts' e_ts'] 
%VE ----------------------------------------------------------------------- 
plot(c5,Exx_e,'ko') 
hold on 
plot(c5,Exx_kd,'kx') 
hold on 
plot(c5,Exx_dn,'k-*') 
113
hold on 
plot(c5,Exx_ts,'k+') 
hold on 
grid on 
xlabel('c_5') 
ylabel('E{x^2}') 
%HAM CON ------------------------------------------------------------------ 
%Tieu chuan doi ngau trong so 
function F_ts =MyFun_ts(Exx,h,sigma,c1,c3,c5,p1,p2) 
Exg1=-3*c3*Exx^2; 
Exg2= 15*c5*Exx^3; 
Egg1=15*c3^2*Exx^3; 
Egg2=945*c5^2*Exx^5; 
miu1=Exg1^2/(Exx*Egg1); 
miu2=Exg2^2/(Exx*Egg2); 
k10=(1-p1)/(1-miu1*p1)*Exg1/Exx; 
k20=(1-p2)/(1-miu2*p2)*Exg2/Exx; 
k=k10+k20; 
F_ts=Exx-sigma^2/(2*(2*h*(c1+k))); 
%Tieu chuan doi ngau 
function F_dn =MyFun_dn(Exx,h,sigma,c1,c3,c5) 
Exg1=-3*c3*Exx^2; 
Exg2= 15*c5*Exx^3; 
Egg1=15*c3^2*Exx^3; 
Egg2=945*c5^2*Exx^5; 
miu1=Exg1^2/(Exx*Egg1); 
miu2=Exg2^2/(Exx*Egg2); 
k10=1/(2-miu1)*Exg1/Exx; 
k20=1/(2-miu2)*Exg2/Exx; 
k=k10+k20; 
F_dn=Exx-sigma^2/(2*(2*h*(c1+k))); 
%Tieu chuan kinh dien 
function F_kd=MyFun_kd(Exx,h,sigma,c1,c3,c5) 
Exg1=-3*c3*Exx^2; 
Exg2= 15*c5*Exx^3; 
Egg1=15*c3^2*Exx^3; 
Egg2=945*c5^2*Exx^5; 
k10=Exg1/Exx; 
k20=Exg2/Exx; 
k=k10+k20; 
F_kd=Exx-sigma^2/(2*(2*h*(c1+k))); 
%Nghiem chinh xac 
function FunTS=Fun1(z,h,sigma,c1,c3,c5) 
FunTS=z.^2.*exp(-(2*h).*(c1*z.^2 - c3*z.^4./2 + c5*z.^6./3)./(sigma.^2)); 
function FunMS=Fun2(z,h,sigma,c1,c3,c5) 
FunMS=exp(-(2*h).*(c1*z.^2 - c3*z.^4./2 + c5*z.^6./3)./ (sigma.^2)); 
%----------------------------------END------------------------------------% 
10. Tính toán mô men bậc hai của dao động có cản phi tuyến phụ thuộc năng 
lượng 
function NL_damping_general_WDC 
clear; 
%clc; 
%THONG SO DAU VAO---------------------------------------------------------- 
ci=1; 
omega=1; 
sigma=sqrt(2); 
114
a= [0 0.039935549 0.129760939 0.210599248 0.307413177 0.386448653 
0.456765361 0.58358966 0.700836619 0.813764905 0.925334967 1.037614507 
1.152308841 1.271032783 1.395497683 1.527689501 1.67008072 1.82596453 
1.999999049 2.199229457 2.435332206 2.729777741 3.130242591 5.227723429 
7.11067138]; 
%GIAI---------------------------------------------------------------------- 
Na=max(size(a)); 
for i=1:1:Na 
 %nghiem chinh xac 
 h_e(i)=(1+a(i)).^(1./(1+a(i)))*gamma(2./(1+a(i)))/gamma(1./(1+a(i))); 
 Exx_e(i)=(1./omega.^2)*(sigma.^2./(2*ci)).^(1./(1+a(i)))*h_e(i); 
 %muc do phu thuoc tuyen tinh 
 miu(i)=gamma(a(i)+2).^2./gamma(2*a(i)+2); 
 if miu(i)>=0 && miu(i)<=1 
 if miu(i)<1/3 
 p(i)=-6.*miu(i)./5+1; 
 elseif miu(i)==1/3 
 p(i)=27./49; 
 elseif miu(i)>1/3 && miu(i)<=2/3 
 p(i)=1./2; 
 elseif miu(i)>2/3 
 p(i)=-3.*miu(i)./2+3./2; 
 end 
 end 
 % Nghiem doi ngau trong so 
 x01=Exx_e(i); %gia tri ban dau 
 Exx_ts(i)=fsolve(@(Exx)MyFun_ts(Exx,ci,omega,a(i),sigma,p(i)),x01); 
 e_ts(i)=abs(Exx_ts(i)-Exx_e(i))*100/Exx_e(i); %sai so 
 % Nghiem doi ngau 
 x02=Exx_e(i); %gia tri ban dau 
 Exx_dn(i)=fsolve(@(Exx)MyFun_dn(Exx,ci,omega,a(i),sigma),x02); 
 e_dn(i)=abs(Exx_dn(i)-Exx_e(i))*100/Exx_e(i); %sai so 
 %Nghiem kinh dien 
 x03=Exx_e(i); %gia tri ban dau 
 Exx_kd(i)=fsolve(@(Exx)MyFun_kd(Exx,ci,omega,a(i),sigma),x03); 
 e_kd(i)=abs(Exx_kd(i)-Exx_e(i))*100/Exx_e(i); %sai so 
end 
%DAU RA-------------------------------------------------------------------- 
format short g 
result=[a' Exx_e' Exx_kd' e_kd' Exx_dn' e_dn' Exx_ts' e_ts' miu' p'] 
%VE------------------------------------------------------------------------ 
plot(a,Exx_e,'k-') 
hold on 
plot(a,Exx_kd,'g-.') 
hold on 
plot(a,Exx_dn,'b-') 
hold on 
plot(a,Exx_ts,'r-') 
hold on 
grid on 
xlabel('a') 
ylabel('') 
%HAM CON------------------------------------------------------------------ 
115
%Tieu chuan doi ngau trong so 
function F_ts=MyFun_ts(Exx,ci,omega,a,sigma,p) 
sigmax=sqrt(Exx); 
sigmay=sigmax*omega; 
EAB=ci*gamma(a+2)*(sigmay^2)^(a+1); 
EAA=ci^2*gamma(2*a+2)*(sigmay^2)^(2*a+1); 
miu=EAB^2/(EAA*sigmay^2); 
b=(1-p)/(1-miu*p)*EAB/(sigmay^2); 
F_ts=Exx-sigma^2/(2*b*(omega^2)); 
%Tieu chuan doi ngau 
function F_dn=MyFun_dn(Exx,ci,omega,a,sigma) 
sigmax=sqrt(Exx); 
sigmay=sigmax*omega; 
EAB=ci*gamma(a+2)*(sigmay^2)^(a+1); 
EAA=ci^2*gamma(2*a+2)*(sigmay^2)^(2*a+1); 
miu=EAB^2/(EAA*sigmay^2); 
b=1/(2-miu)*EAB/(sigmay^2); 
F_dn=Exx-sigma^2/(2*b*(omega^2)); 
%Tieu chuan kinh dien 
function F_kd=MyFun_kd(Exx,ci,omega,a,sigma) 
sigmax=sqrt(Exx); 
sigmay=sigmax*omega; 
EAB=ci*gamma(a+2)*(sigmay^2)^(a+1); 
EAA=ci^2*gamma(2*a+2)*(sigmay^2)^(2*a+1); 
b=EAB/(sigmay^2); 
F_kd=Exx-sigma^2/(2*b*(omega^2)); 
%----------------------------------END------------------------------------% 
11. Tính toán tần số của dao động tự do 
function Nonlinear_spring_periodic 
clear; 
%clc; 
%THONG SO DAU VAO---------------------------------------------------------- 
gamma0=1; 
a=0; 
b=2*pi; 
n=[3 5 7 9 11 13 15]; 
%GIAI---------------------------------------------------------------------- 
Nn=max(size(n)); 
for i=1:1:Nn 
 T_e(i)=Fun_e(gamma0,n(i)); 
 omega_e(i)=2*pi./T_e(i); %tan so chinh xac 
 k_e(i)=omega_e(i).^2; 
 %tinh toan cac mo men 
 EAB(i)=1/(b-a)*quad(@(x)gamma0*cos(x).^(n(i)+1),a,b); 
 EB2(i)=1./(b-a)*quad(@(x)cos(x).^(2),a,b); 
 EA2(i)=1./(b-a)*quad(@(x)gamma0.^2.*cos(x).^(2*n(i)),a,b); 
 %muc do phu thuoc tuyen tinh 
 miu(i)=EAB(i).^2./(EA2(i).*EB2(i)); 
 if miu(i)>=0 && miu(i)<=1 
 if miu(i)<1/3 
 p(i)=-6.*miu(i)./5+1; 
 elseif miu(i)==1/3 
 p(i)=27./49; 
116
 elseif miu(i)>1/3 && miu(i)<=2/3 
 p(i)=1./2; 
 elseif miu(i)>2/3 
 p(i)=-3.*miu(i)./2+3./2; 
 end 
 end 
 %Nghiem kinh dien 
 k_kd(i)=EAB(i)./EB2(i); 
 omega_kd(i)=sqrt(k_kd(i)); 
 e_kd(i)=abs(omega_kd(i)-omega_e(i))*100/omega_e(i); 
 % Nghiem doi ngau 
 k_dn(i)=1./(2-miu(i)).*k_kd(i); 
 omega_dn(i)=sqrt(k_dn(i)); 
 e_dn(i)=abs(omega_dn(i)-omega_e(i))*100./omega_e(i); 
 % Nghiem doi ngau trong so 
 k_ts(i)=(1-p(i))./(1-miu(i).*p(i)).*k_kd(i); 
 omega_ts(i)=sqrt(k_ts(i)); 
 e_ts(i)=abs(omega_ts(i)-omega_e(i))*100./omega_e(i); 
end 
%DAU RA-------------------------------------------------------------------- 
format short g 
result=[n' omega_e' omega_kd' e_kd' omega_dn' e_dn' omega_ts' e_ts' miu'] 
%HAM CON------------------------------------------------------------------
function FunT_e=Fun_e(gamma0,n) 
FunT_e=4*quad(@(z)1./(sqrt(2*(gamma0*(1-z.^(n+1))./(n+1)))),0,1); 
%----------------------------------END------------------------------------%