Ứng dụng phần mềm cad/cam cimatron trong thiết kế, chế tạo khuôn mẫu

nh nghĩa trước độ dốc của hàm mục tiêu. Giá trị của thông số này là 'on' và 'off' (mặc định). Hessian: Bật, tắt chế độ định nghĩa trước đạo hàm riêng bậc hai của hàm mục tiêu. Giá trị của thông số này là 'on' và 'off' (mặc định). GradConstr: Bật, tắt chế độ định nghĩa trước độ dốc của hàẩìng buộc. Giá trị của thông số này là 'on' và 'off' (mặc định). TolCon: Sai số cho phép đối với các điều kiện ràng buộc. Giá trị của nó là số nguyên dương. TolFun: Sai số cho phép của hàm mục tiêu (đối với bài toán tối ưu) và đối với các hàm vế trái của hệ. Giá trị của nó là số nguyên dương. 200 TolPCG: Sai số cho phép để dừng bước lặp PCG. Giá trị của nó là số nguyên dương. Mặc định là 0.1. TolX: Quy định sai cho phép của các biến trạng thái. Giá trị của nó là một số nguyên dương. e. Các ví dụ Ví dụ 1: Tìm điểm 0 của hệ hai phương trình 2 ẩn Chuyển hệ về dạng chuẩn: Điểm ước lượng ban đầu là : x0 = [-5 -5]. Bước 1: Trước hết chúng ta phải viết hàm M-file để tính toán ra F (giá trị của các phương trình tại x) : function F = myfun(x) F = [2*x(1) - x(2) - exp(-x(1)); -x(1) + 2*x(2) - exp(-x(2))]; Bước 2: Gọi ra lệnh fsolve để giải hệ: x0 = [-5; -5]; % ước lượng nghiệm ban đầu options=optimset('Display','iter'); % Tuỳ chọn hiển thị ra màn hình [x,fval] = fsolve('myfun',x0,options) % Gọi công cụ tối ưu Sau 28 bước lặp máy sẽ tìm ra nghiệm: 201 Bước lặp Hàm đếm f(x) Chuẩn của bước Tuỳ chọn bậc nhất Các bước lặp CG 1 4 47071.2 1 2.29e+004 0 2 7 6527.47 1.45207 3.09e+003 1 3 10 918.372 1.49186 418 1 4 13 127.74 1.55326 57.3 1 5 16 14.9153 1.57591 8.26 1 6 19 0.779051 1.27662 1.14 1 7 22 0.00372453 0.484658 0.0683 1 8 25 9.21617e-008 0.0385552 0.000336 1 9 28 5.66133e-017 0.000193707 8.34e-009 1 Ví dụ 2: Tìm ma trận X thoả mãn hệ phương trình Điểm ước lượng ban đầu: x= [1,1; 1,1]. Bước 1: Viết M-file mô tả hệ phương trình: function F = myfun(x) F = x*x*x-[1,2;3,4]; Bước 2: Gọi lệnh fsolve: x0 = ones(2,2); % ước lượng ban đầu của nghiệm options = optimset('Display','off'); % Tắt chế độ hiển thị quá trình tính [x,Fval,exitflag] = fsolve('myfun',x0,options) 202 Nghiệm của phương trình là: x = -0.1291 0.8602 1.2903 1.1612 Fval = 1.0e-03 * 0.1541 -0.1163 0.0109 -0.0243 exitflag = 1 Chú ý: Nếu hệ phương trình là tuyến tính, toán tử \ sẽ cho kết quả chính xác hơn và nhanh hơn lệnh fsolve. Ví dụ để giải hệ sau: Chúng ta sẽ thực hiện các câu lệnh của Matlab sau đây: A = [ 3 11 -2; 1 1 -2; 1 -1 1]; b = [ 7; 4; 19]; x = A\b Kết quả sẽ nhận được nghiệm: x = 13.2188 -2.3438 3.4375 3. Tìm cực tiểu không ràng buộc của hàm nhiều biến Để tìm giá trị cực tiểu không ràng buộc của hàm nhiều biến trước hết chúng ta phải viết một M-file định nghĩa hàm số. Hàm này có đối số là một véc tơ có các phần tử là các biến độc lập của hàm. Sau đó chúng ta sử dụng hàm fminsearch của Matlab để tìm ra giá trị cực tiểu của hàm ở lân cận điểm ước lượng ban đầu. 203 fminsearch sử dụng thuật toán tìm kiếm đơn hình. Đó là phương pháp tìm kiếm trực tiếp không sử dụng phương pháp xác định độ dốc bằng phương pháp giải tích hay phương pháp số. Nếu n là chiều dài của véc tơ biến độc lập x, thuật đơn hình trong không gian n chiều được đặc trưng bởi n+1 véc tơ khác biệt nằm trên đỉnh của chúng. Trong không gian 2 chiều, đó là một tam giác, trong không gian 3 chiều đó là tứ diện ... Tại mỗi bước tìm kiếm, một điểm mới nằm trong đơn hình hoặc gần nó sẽ được tạo ra. Giá trị của hàm tại điểm mới được so sánh với giá trị của hàm tại các đỉnh của đơn hình. Thông thường một trong các đỉnh của đơn hình được thay thế bằng điểm mới và hình thành nên đơn hình mới. Quá trình này sẽ lặp lại cho tới khi bán kính của đơn hình nhỏ hơn giá trị cho phép. a. Sử dụng lệnh fminsearch fminsearch có cú pháp cơ sở như sau: x = fminsearch('fun',x0) Trong đó: o fun là chuỗi xác định tên của M-file định nghĩa hàm cần tối ưu o X0 là véc tơ các giá trị ước lượng ban đầu của các biến độc lập Chú ý: 1. Tương tự như ở lệnh fsolve, để có thêm các thông tin về quá trình tính toán ta có thể dùng cú pháp: [x,fval,exitflag,output] = fminsearch ('fun',x0) Trong đó: fval,exitflag,output có ý nghĩa giống như ở lệnh fsolve. 2. Để điều khiển quá trình tính toán của lệnh chúng ta sử dụng cú pháp : x = fminsearch('fun',x0, options) Trong đó: options được xác lập bằng lệnh optimset (xem ở phần trên) với các thuộc tính được sử dụng ở đây là: Display, MaxFunEvals, MaxIter, TolFun, TolX b. Các ví dụ Ví dụ 1: Xác định cực tiểu của hàm một biến: f(x) = sin(x) + 3. Bước 1: Viết M-file mô tả hàm số: function f = myfun(x) f = sin(x) + 3; Bước 2: Gọi hàm fminsearch: 204 x = fminsearch('myfun',2) Ví dụ 2: Xác định cực trị của hàm số Rosenbrock: Bằng giải tích ta có thể xác định được hàm có giá trị cực tiểu bằng không tại điểm (1,1). Bước 1: Định nghĩa M-file có tên là banana.m sẽ dùng để định nghĩa hàm cần tối ưu function f = banana(x) f = 100*(x(2)-x(1)^2)^2+(1-x(1))^2; Bước 2: gọi hàm fminsearch: [x,out] = fminsearch ('banana',[-1.2, 1.2 ]); Sẽ cho kết quả là: x= 1.0000 1.0000 out = 5.4009e-010 4. Tìm cực tiểu của hàm một biến trong vùng xác định Bài toán xác định cực tiểu của hàm một biến trong vùng xác định được viết ở dạng toán học như sau: 21x xxx:)x(fMin << Trong đó : x, x1, x2 là các giá trị vô hướng, f là hàm cho kết quả là giá trị vô hướng. a. Sử dụng hàm fminbnd Hàm fminbnd được dùng để xác định cực tiểu của hàm một biến trong một khoảng xác định. fminbnd sử dụng thuật toán "tìm kiếm mặt cắt vàng" và phép nội suy bậc hai. Các cú pháp của lệnh fminbnd bao gồm: x = fminbnd('fun',x1,x2) x = fminbnd('fun',x1,x2,options) x = fminbnd('fun',x1,x2,options,P1,P2,...) [x,fval] = fminbnd(...) [x,fval,exitflag] = fminbnd(...) [x,fval,exitflag,output] = fminbnd(...) Trong đó: 205  'fun', fval, exitflag, output, options: có ý nghĩa giống như ở lệnh fminsearch  x1, x2 : là giá trị cận dưới và trên của biến tìm kiếm b. Ví dụ Xác định cực tiểu của hàm f= X3-2X-5 ở trong khoảng [ 0 2]. Bước1: viết M.file định nghĩa hàm f: function y = f(x) y = x.^3-2*x-5; Bước2: Gọi hàm fminbnd [x,out] = fminbnd('f', 0, 2) Kết quả nhận được là: x = 0.8165, out=-6.0887 5. Tìm cực trị có điều kiện của hàm nhiều biến Bài toán tìm cực trị có điều kiện của hàm nhiều biến trong trường hợp tổng quát được xác định như sau: Trong đó:  x, b, beq, lb, và ub là các véc tơ,  A và Aeq là các ma trận,  c(x) và ceq(x) là các hàm nhiều biến đưa ra kết quả là các véc tơ,  f(x) là hàm nhiều biến đưa ra kết quả là giá trị vô hướng,  f(x), c(x), và ceq(x) có thể là các hàm phi tuyến. Matlab đưa ra hàm fmincon để giải quyết bài toán tìm cực trị có ràng buộc của hàm nhiều biến. Tuỳ theo yêu cầu về các điều kiện ràng buộc, chúng ta sử dụng các cú pháp khác nhau của lệnh fmincon. a. Các cú pháp của lệnh fmincon x = fmincon(fun,x0,A,b) x = fmincon(fun,x0,A,b,Aeq,beq) 206 x = fmincon(fun,x0,A,b,Aeq,beq,lb,ub) x = fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon) x = fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options) [x,fval] = fmincon(...) [x,fval,exitflag] = fmincon(...) [x,fval,exitflag,output] = fmincon(...) [x,fval,exitflag,output,lambda] = fmincon(...) [x,fval,exitflag,output,lambda,grad] = fmincon(...) [x,fval,exitflag,output,lambda,grad,hessian] = fmincon(...) Trong đó:  x = fmincon(fun,x0,A,b) : Xuất phát từ điểm tìm kiếm ban đầu x0 và tìm ra điểm cực tiểu của hàm số mô tả trong M-file 'fun' với ràng buộc biểu diễn bằng hệ bất phương trình: A*x <= b.  x = fmincon(fun,x0,A,b,Aeq,beq): tìm ra điểm cực tiểu của hàm số mô tả trong M-file 'fun' với ràng buộc biểu diễn bằng hệ bất phương trình: A*x <= b, hệ phương trình Aeq*x = beq. Nếu không có các ràng buuộc dạng hệ bất phương trình, chúng ta phải đưa vào các giá trị rỗng: A=[ ] và b=[ ].  x = fmincon(fun,x0,A,b,Aeq,beq,lb,ub) tìm ra điểm cực tiểu của hàm số mô tả trong M-file 'fun' với các biến thiết kế bị ràng buộc bởi vùng lb <= x <= ub. Nếu không có thêm các ràng buộc dạng hệ phương trình, hãy đưa vào các giá trị rỗng Aeq=[ ] và beq=[ ].  x=fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon):Tìm cực tiểu của hàm số có thêm các ràng buộc ở dạng các hệ phương trình phi tuyến ceq(x) = 0 và hệ các bất phương trình phi tuyến c(x) <= 0 được mô tả trong M-file 'nonlcon'. Nếu không có ràng buộc về điều kiện biên hãy gán lb=[ ] và (hoặc) ub=[ ].  x = fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options): cho phép đưa vào các thông số điều khiển quá trình tối ưu thông qua biến options được xác lập bằng lệnh optimset (xem cách sử dụng trong lệnh fsolve). Các thông số đầu ra của lệnh:  x, fval, exitflag, output: có ý nghĩa hoàn toàn giống như ở lệnh fsolve. 207  lambda: Là biến dạng cấu trúc chứa các nhân tử Lagrange tại điểm cực tiểu của hàm:  lambda.lower : Nhân tử đối với biên cận dưới lb.  lambda.upper : Nhân tử đối với biên cận trên ub.  lambda.ineqlin: Nhân tử đối với ràng buộc dạng hệ bất phương trình tuyến tính.  lambda.eqlin: Nhân tử đối với ràng buộc dạng hệ phương trình tuyến tính.  lambda.ineqnonlin: Nhân tử đối với ràng buộc dạng hệ bất phương trình phi tuyến.  lambda.eqnonlin: Nhân tử đối với ràng buộc dạng hệ phương trình phi tuyến.  grad : véc tơ giá trị độ dốc của hàm mục tiêu tại điểm cực trị tìm được.  hessian: véc tơ các giá trị đạo hàm riêng bậc hai của hàm mục tiêu tại điểm cực trị tìm được. b. Cách viết M-file mô tả hàm mục tiêu Hàm mục tiêu trong lệnh fmincon là một hàm nhiều biến đưa ra kết quả là giá trị vô hướng cần được tối thiểu hoá. Khi viết M-file mô tả hàm mục tiêu, chúng ta phải đưa vào đối số là một biến véc tơ có các phần tử là các biến thiết kế: Cấu trúc cơ bản của hàm mục tiêu: function f = myfun(x) f = ... % Tính giá trị của hàm tại điểm x Định nghĩa độ dốc trong hàm mục tiêu Lệnh fmincon cũng cho phép ta định nghĩa trước độ dốc của hàm mục tiêu. Độ dốc ở một điểm được hiểu là đạo hàm riêng của hàm mục tiêu theo các biến thiết kế tại điểm đó. Để có thể định nghĩa trước độ dốc, chúng ta phải xử lý hai vấn đề sau: * Bật chế độ định nghĩa trước độ dốc trong câu lệnh của fmincon: Sử dụng cú pháp: x = fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options) Trong đó options được xác lập bằng lệnh otimset như sau: 208 options = optimset('GradObj','on') * Thêm biến đầu ra của M-file định nghĩa hàm mục tiêu: function [f,g] = myfun(x) f = ... % tính giá trị hàm tại x if nargout > 1 % gọi hàm khi có 2 biến đầu ra g = ... % tính toán độ dốc tại x end; Định nghĩa trước đạo hàm bậc hai của hàm mục tiêu (ma trận hessian của hàm mục tiêu) Tương tự như việc định nghĩa trước độ dốc, để định nghĩa trước đạo hàm riêng phần bậc hai của hàm mục tiêu theo các biến thiết kế chúng ta cũng phải xử lý hai vấn đề: * Bật chế độ định nghĩa trước độ dốc trong câu lệnh của fmincon: Sử dụng cú pháp: x = fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options) Trong đó options được xác lập bằng lệnh otimset như sau: options = optimset('Hessian','on') * Thêm biến đầu ra của M-file định nghĩa hàm mục tiêu: function [f,g,H] = myfun(x) f = ... % Tính giá trị hàm mục tiêu tại x if nargout > 1 % hàm được gọi khi có 2 biến ra g = ... % Nhập độ dốc của hàm mục tiêu tại x if nargout > 2 H = ... % Nhập ma trận hessian tại x end; end; c. Cách viết M-file mô tả các ràng buộc phi tuyến Các ràng buộc phi tuyến ở dạng hệ phương trình hay bất phương trình cần được mô tả trong M-file (với tên gọi là 'nonlicon' trong lệnh fmincon). Trong 209 cấu trúc cơ bản của M-file này thông số đầu vào của hàm là một biến véc tơ có các phần tử là các biến thiết kế. Biến ra bao gồm 2 biến véc tơ, biến thứ nhất lấy ra các ràng buộc bất phương trình, biến thứ hai lấy ra các ràng buộc dạng phương trình: function [c,ceq] = mycon(x) c = ... % Tính các ràng buộc dạng hệ bất phương trình tại x ceq = ... % Tính các ràng buộc dạng hệ phương trình tại x Định nghĩa trước độ dốc của các ràng buộc Tương tự như việc định nghĩa trước mục tiêu, để định nghĩa trước độ dốc của hàm ràng buộc theo các biến thiết kế chúng ta cũng phải xử lý hai vấn đề: * Bật chế độ định nghĩa trước độ dốc trong câu lệnh của fmincon: Sử dụng cú pháp: x = fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options) Trong đó options được xác lập bằng lệnh otimset như sau: options = optimset('GradConstr','on') * Thêm biến đầu ra của M-file định nghĩa hàm ràng buộc: function [c,ceq,GC,GCeq] = mycon(x) c = ... % Tính các ràng buộc dạng hệ bất phương trình tại x ceq = ... % Tính các ràng buộc dạng hệ phương trình tại x if nargout > 2 % hàm được gọi khi có 4 biến ra GC = ... % gradients của các bất phương trình tại x GCeq = ... % gradients của các phương trình tại x end; d. Các ví dụ Ví dụ 1: Tìm điểm cực tiểu của hàm , xuất phát tại điểm x = [10; 10; 10] với ràng buộc : Bước 1: Viết M-file mô tả hàm mục tiêu: 210 function f = myfun(x) f = -x(1) * x(2) * x(3); Bước 2: Chuyển các ràng buộc tuyến tính về dạng chuẩn: Được chuyển về dạng ma trận sau: với Bước 3: gọi lệnh trong Matlab x0 = [10; 10; 10]; % điểm xuất phát A=[-1, -2, -2; 1, 2, 2]; b=[0; 72]; [x,fval] = fmincon('myfun',x0,A,b) Sau 66 lần tính toán Matlab sẽ đưa ra kết quả: x = 24.0000 12.0000 12.0000 Với giá trị cực tiểu của hàm là fval = -3.4560e+03 Ví dụ 2: Tính toán thông số tối ưu của nhíp lá Nhíp lá là bộ phận đàn hồi thường dùng cho các xe tải. Nó gồm nhiều lá thép có chiều dài thay đổi được bó với nhau thành một khối và được ngàm chặt vào cầu xe bằng các bu lông quang nhíp. Tải trọng tác động lên nhíp đặt tại các quang treo ở đầu mút của cụm nhíp. Công việc tối ưu nhíp bao gồm việc xác định chiều dài, chiều dầy các lá nhíp để đạt được mục đích tiết kiệm vật liệu nhất và vẫn đảm bảo các điều kiện sau:  ứng suất trong các lá nhíp không vượt quá giá trị cho phép  Độ cứng của nhíp không vượt quá giá trị yêu cầu (để đảm bảo độ êm dịu chuyển động) 211 Các thông số của xe là: TT Thông số Ký hiệu (đơn vị đo) giá trị 1 Trọng lượng toàn bộ của xe m (Kg) 7400 2 Chiều dài cơ sở l (m) 3.7 3 Toạ độ trọng tâm a, b (m) 2.795; 0.905 4 Bề rộng của lá nhíp b (Cm) 6.5 5 1/2 Chiều dài của lá đầu tiên l1 (Cm) 58.3 6 ứng xuất cho phép [σ] (Kg/cm2) 9400 Bước 1: Xác định các thông số thiết kế, hàm mục tiêu Các biến thiết kế:  Bề dầy các nhóm lá nhíp: h1, h2 ...  Số lá nhíp: n  Chiều dài từng lá nhíp: li với i=2→n Hàm mục tiêu: Trong trường hợp các lá nhíp có tiết diện chữ nhật, hàm mục tiêu sẽ là:     ρ=∑ = n 1k iii21 l.h.b.)l,h,h,n(FMin Giá trị ban đầu của biến thiết kế: được xác định bằng phương pháp thiết kế cổ điển của Parchilovskii: Chỉ số lá Bề dầy lá (cm) 1/2 Chiều dài lá nhíp Chỉ số lá Bề dầy lá (cm) 1/2 Chiều dài lá nhíp 1,2 0.7 58.3 8 0.7 23.03 3 0.7 46.58 9 0.7 18.32 4 0.7 41.87 10 0.7 13.61 5 0.7 37.16 11 0.7 10.25 6 0.7 32.45 12 0.7 8.9 7 0.7 27.74 13 0.7 4.2 Bước 2: Xác định các ràng buộc Ràng buộc về độ cứng: 212 Từ tính toán dao động chúng ta đã xác định được độ cứng của nhíp cần phải không được vượt quá giá trị C1=85 Kg/cm. Để xác định độ cứng theo biến thiết kế, chúng ta sử dụng công thức bán thực nghiệm của Parchilovskii: ∑ ++ − α= n k 1kk 3 1k )YY(a .E6 C Trong đó: + α là hệ số giảm độ cứng, thường có giá trị = 0.83 đến 0.87 + n là số lá nhíp, E là mô đun đàn hồi (Kg/cm2) + ai là độ lệch chiều dài các lá nhíp: a1=li-1-li ; i=2 →n + Yi=1/Ji ; Ji là mô men quán tính của mặt cắt ngang của phần thứ i của khối nhíp; nếu mặt cắt ngang của các lá nhíp có dạng hình chữ nhật, Ji được xác định bằng biểu thức : ∑ = = i 1k 3 i i 12 h.b J Ràng buộc về ứng suất: ứng suất trong các lá nhíp được xác định trên cơ sở tính toán các giá trị tải trọng đặt trên đầu mút của mỗi lá nhíp. Theo Parchilovskii, tải trọng này được xác định bằng hệ phương trình:        =+ =++ =++ =++ − 0XBXA .......................................... 0XCX.BX.A 0XCX.BX.A 0XCX.BP.A nn1nn 544434 433323 32222 Trong đó các hệ số Ak, Bk, Ck được xác định bằng các biểu thức: n2k l )ll( J J 1B;n2k;1 l l .3. J J 5.0A 3 k 3 1kk 1k k k k 1k 1k k k →=   −η++−=→=    −= + − − − 1n2k.1 l l 3. l l C 1k k 3 k 1k k −→=    −   = + + Trong biểu thức trên η là hệ số xét đến hình dạng đầu lá nhíp (tra trong bảng) [ 1 ]. Trong trường hợp đầu lá nhíp không bị vuốt ta có η=0. Từ giá trị của phản lực chúng ta sẽ xác định được biểu đồ mô men của mỗi lá, tìm ra mô men lớn nhất và từ đó xác định được giá trị ứng xuất lớn nhất trên mỗi lá. Các điều kiện ràng buộc khác: 213  Chiều dài của lá nhíp đầu tiên được xác định trên cơ sở bố trí chung trên xe: l1=58.3 cm  Chiều dài lá nhíp thứ k+1 phải ngắn hơn chiều dài lá nhíp thứ k: lk+1- lk≤0 với k=1→n-1  Các giá trị về chiều dài, bề dầy của các lá nhíp phải là số dương: h1>0, h2 >0, lk>0 với k=2→n Các điều kiện này có thể được viết dưới dạng chung:   = ≤ 0)x(Ceq 0)x(C , với x là véc tơ biến thiết kế Bước 3: Viết các M-file mô tả hàm mục tiêu và hàm ràng buộc M-file mô tả hàm mục tiêu: function y=nhip(x) n=13; h=x(3)*ones(1,n);h(1)=x(2); h(2)=x(2); b=x(1); l=x(4:end); y=sum(b*h.*l)*0.0087; M-file mô tả hàm ràng buộc: function [y,yy]=rangbuoc(x) nn=13; h=x(3)*ones(1,nn);h(1)=x(2); h(2)=x(2); b=x(1);al=0.85;E=2000000; l=x(4:end); l(nn+1)=0; %% Tính toán các đặc tính hình học của lá nhíp Jk=b*h.^3/12;Jk=cumsum(Jk); Yk=1./Jk;Yk(nn+1)=0;DYk=-diff(Yk); a=l(1)*ones(1,nn+1)-l;a=a(2:nn+1); %% Tính độ cứng của nhíp Ck=al*6*E/sum(a.^3.*DYk) yy(1)=Ck-85; %% ràng buộc về độ cứng %% Tính ứng suất lớn nhất ở các lá nhíp Q=1110*1.8;l=l(1:end-1); Jk=b*h.^3/12; Ak=0.5*Jk(2:end)./Jk(1:end-1).*(3*l(1:end-1)./l(2:end)-1); Bk=-(1+Jk(2:end)./Jk(1:end-1)); Ck=0.5*(l(3:end)./l(2:end-1)).^3.*(3*l(2:end-1)./l(3:end)-1); z=size(l);z=z(2)-1; V=zeros(z,1); V(1)=-Q/2*Ak(1); M=diag(Ak(2:end),-1)+diag(Bk)+diag(Ck,1); X=inv(M)*V; P=[Q/2 X'];lbc=l(1:end-1);lbd=l(2:end); 214 gt1=P(1:end-1).*(lbc-lbd); gt2=P(1:end-1).*lbc-P(2:end).*lbd; gt=[gt1;gt2];gtm=max(gt); %% ứng suất lớn nhất gtm=[gtm P(end)*l(end)]; us=gtm./Jk.*h/2 y(1:z+1)=us-9500; %% ràng buộc về ứng suất y=[y -x+0.004]; %% thêm các ràng buộc về chiều dài dương yy(2)=58.3-x(4); %% ràng buộc về chiều dài lá thứ nhất yy(3)=6.5-x(1); %% ràng buộc về chiều rộng các lá yyk=l(2:end)-l(1:end-1); %% ràng buộc lá sau ngắn hơn lá trước y=[y yyk]; %% ràng buộc dạng bất phương trình Bước 4: gọi lệnh fmincon từ Matlab x0=[6.5, 0.7, 0.7, 58.3, 46.58, 41.87, 37.16, 32.45, 27.74, ... 23.03, 18.32, 13.61, 10.25, 8.9, 4.2] [x, fval] = fmincon('nhip',x0,[ ], [ ], [ ], [ ], [ ],[ ],'rangbuoc') Chúng ta xác định được các giá trị tối ưu: Chỉ số lá Bề dầy lá (cm) 1/2 Chiều dài lá nhíp (Cm) 1 1.0911 58.30 2 1.0911 37.39 3 0.7530 17.94 4 0.7530 8.68 5 0.7530 2.79 Dưới đây là một số chỉ tiêu đánh giá chất lượng hai bộ nhíp dước tác dụng của tải trọng 1110 Kg (tải trọng lên một nhíp): TT Chỉ tiêu so sánh Nhíp sau khi tối ưu Nhíp cũ 1 Khối lượng (1/2 nhíp không kể phần ngàm) 7.157 (Kg) 15.456 (Kg) 2 Độ cứng 84.53 Kg/cm 87.62 Kg/cm 3 ứng suất lớn nhất 9140 Kg/cm2 10605 Kg/cm2 215 So sánh các chỉ tiêu đánh giá ta nhận thấy nhíp sau khi được tối ưu hoá có được các đặc tính vượt trội hơn hẳn nhíp tính sơ bộ như: Khối lượng giảm đi 50%, độ cứng giảm 4% (đảm bảo tốt hơn tính êm dịu chuyển động) và ứng suất ở trong các lá nhíp đã tối ưu giảm 14%. Điều này khẳng định được hiệu quả rõ rệt của phương pháp tối ưu đưa ra. 6.3. Giải các hệ phương trình vi phân thường 6.3.1. Giới thiệu về hệ phương trình vi phân thường a. Dạng chuẩn của hệ phương trình vi phân Phần lớn các hệ động lực đều được mô tả trong dạng toán học chuẩn có tên gọi là dạng Cauchy, đó là hệ n các phương trình vi phân bậc nhất sau:    =′ =′ )w,t(f)t(W ............ )w,t(f)t(W )n()n( )()( 11 (2-1) Trong đó :  t là biến độc lập (thường là thời gian),  w(i) (t) là thành phần thứ i của véc tơ W(t) có n phần tử và đạo hàm của nó theo biến t được ký hiệu bằng w'(t).  f(i)(t) là thành phần thứ i của véc tơ n hàm phụ thuộc vào t, w: F(t,w) Trong dạng véc tơ chúng ta có thể viết : W'(t)=F(t,w) (2-2) Các lệnh giải phương trình vi phân trong Matlab đều sử dụng dạng chuẩn này, chính vì vậy chúng ta phải chuyển các hệ phương trình vi phân về dạng Cauchy trước khi áp dụng các lệnh giải phương trình vi phân của Matlab. b. Dạng hệ phương trình vi phân có ma trận khối lượng Trong thực tế tính toán, không phải lúc nào chúng ta cũng có thể dễ dàng chuyển hệ phương trình vi phân mô tả hệ thống về dạng tường minh như trên. Chính vì vậy từ phiên bản 5.0, Matlab đã cho phép chúng ta giải các hệ phương trình vi phân ở dạng ẩn (có ma trận khối lượng) như sau: M.W'(t)=F(t,w) (2-3) Trong đó: 216 M là ma trận vuông n,n : các phần tử của ma trận này có thể là hằng số, có thể là hàm phụ thuộc vào t : M(t) hoặc phụ thuộc vào thời gian và biến trạng thái M(t,w). c. Các phương trình cứng (stiff equation) Các hệ phương trình vi phân có các giá trị riêng khác biệt rất lớn được gọi là những hệ phương trình cứng. Khác với các hệ phương trình vi phân thông thường, hệ phương trình vi phân cứng yêu cầu sự khác biệt rất lớn của giá trị bước thời gian khi sử dụng phương pháp số để giải chúng. Chính vì vậy một số phương pháp số dùng để giải hệ phương trình vi phân thông thường tỏ ra không hiệu quả khi dùng để giải các hệ phương trình vi phân cứng. Ví dụ về hệ phương trình vi phân cứng Chúng ta cần xác định nghiệm theo thời gian của hệ phương trình:   =−−=′ =+=′ 0)0(1999999 1)0(1998998 2212 1211 XXXX XXXX Trong dạng không gian trạng thái, phương trình được viết: x'= A.x +B.u Trong đó: 1999999 1998998 −−=A , 0 0=B Các giá trị riêng của ma trận A là λ1=-1, λ2=-1000 Nghiệm chính xác của phương trình có dạng: tt tt eCeCX eCeCX 1000 432 1000 211 −− −− += += Chúng ta nhận thấy rằng hệ thống có 2 thành phần có hằng số thời gian khác biệt rất lớn. Bước thời gian cần thiết khi tích phân số phương trình được điều khiển bởi thành phần hằng số thời gian có giá trị lớn hơn, để đảm bảo độ chính xác chúng ta phải chọn bước thời gian rất nhỏ: ví dụ như nếu chúng ta chọn độ lớn của bước tính bằng 1/10 giá trị ngịch đảo của λ2 thì giá trị bước thời gian phải là ∆t=0.00001 s. Tuy nhiên ảnh hưởng của thành phần thứ 2 này giảm rất nhanh khi thời gian tăng lên. Như vậy khi thời gian đủ lớn chúng ta có thể tăng giá trị bước thời gian lên nhiều lần để đẩy nhanh tốc độ tính toán. 217 6.3.2. Phương pháp giải hệ phương trình Cauchy trong MATLAB Để giải hệ phương trình vi phân trong Matlab chúng ta phải thực hiện 3 bước sau đây: a. Chuyển đổi hệ phương trình vi phân về dạng chuẩn. Dạng chuẩn ở đây được hiểu theo nghĩa rộng là dạng Cauchy ( 2-2) hoặc dạng mở rộng có ma trận khối lượng M, M(t) hay M(t,w). Khi chuyển đổi về dạng chuẩn trước hết chúng ta phải đổi biến để chuyển các đạo hàm bậc cao thành các đạo hàm bậc nhất, sau đó chúng ta chuyển hệ về dạng véc tơ (dạng 2-2) và chuyển điều kiện đầu về cũng về dạng véc tơ Chú ý: chúng ta phải chú ý đến trật tự đặt các biến trong véc tơ biến trạng thái để không nhầm lẫn chúng khi gán các điều kiện đầu cũng như khi xử lý kết quả. Ví dụ: có hệ phương trình vi phân mô tả chuyển động của cơ hệ sau: )t10(Sin.eYY.e.2Y).tcos1( t30t.152 =+′′+′′′+ (2-4) Với các điều kiện đầu là: Y(0)=0; Y'(0)=1; Y"(0)=2 Để chuyển về dạng chuẩn, chúng ta đổi các biến: Y=Y1; Y'=Y2; Y"=Y3 (2-5) Khi đó chúng ta sẽ có: Y'1=Y2; Y'2=Y3; (2-6) Thế các phương trình 2-6 vào phương trình 2-4 và chuyển vế sẽ được: tcos1 Y.e2Y)t10(Sin.eY 2 3 t15 1 t30 3 + −−=′ (2-7) Kết hợp các phương trình 2-6 và 2-7, chúng ta nhận được hệ ở dạng chuẩn:        + −−=′ =′ =′ tcos1 Y.e2Y)t10sin(.eY YY YY 2 3 t15 1 t30 3 32 21 (2-8) Véc tơ điều kiện đầu tương ứng sẽ là [0, 1, 2]T b. Viết M-file mô tả hệ phương trình vi phân (ODE File) ODE File là một tệp dạng m code do người dùng viết để định nghĩa các phương trình vi phân dùng cho các lệnh giải phương trình vi phân của Matlab 218 (lệnh ODE) tương ứng. M- file này sẽ được tham chiếu đến trong các câu lệnh ODE với tên 'ODE File' tuy nhiên chúng ta có thể đặt cho nó các tên bất kỳ. Mặc định các lệnh ODE sẽ giải các bài toán với điều kiện đầu ở dạng: dy/dt=F(t,y), ở đây t là một biến độc lập và Y là một véc tơ của các biến phụ thuộc. Để thực hiện điều này các lệnh ODE trong quá trình giải phương trình sẽ gọi ra hàm F = odefile(t,Y), trong đó t là dẫy số vô hướng, Y là véc tơ cột và kết quả đưa ra của hàm là một véc tơ cột. Chú ý rằng tệp ODE phải chấp nhận các đối số t và Y cho dù không dùng đến các đối số này. Trong trường hợp đơn giản tệp ODE được viết ở dạng: function F = odefile(t,y) F = ; Ví dụ: Để giải phương trình 2-4 đã được chuyển thành dạng chuẩn 2-8, chúng ta cần phải viết Odefile như sau: function F=Ptvp(t,Y) F(1,1)= Y(2); F(2,1)=Y(3); F(3,1)=(exp(30*t)*sin(10*t)-Y(1)-2*exp(15*t)*Y(3))/(1+cos(t)^2) Và ghi vào tệp có tên là Ptvp.m c. Gọi lệnh ODE Matlab cung cấp cho chúng ta 7 lệnh giải phương trình vi phân khác nhau. Đặc điểm và phạm vi sử dụng của các lệnh này được dẫn ra trong bảng 2-1: Lệnh Kiểu bài toán Bậc chính xác Phạm vi sử dụng ode45 Không cứng Trung bình Thường dùng nhất, hay dùng cho lần thử đầu tiên ode23 Không cứng Thấp Dùng để ước lượng thô sai số ode113 Không cứng thấp đến cao Dùng khi yêu cầu về sai số tương đối nghiêm khắc hoặc khi tính toán trên một tệp ODE cần nhiều thời gian tính toán.. 219 ode15s Cứng Thấp đến trung bình Dùng để giải hệ cứng hoặc khi có ma trận độ cứng ode23s Cứng thấp Dùng ước lượng sai số khi giải các hệ cứng hặc khi có ma trận khối lượng hằng số. ode23t Cứng trung bình Thấp Nếu bài toán chỉ hơi cứng và khi cần lời giải không có giảm chấn. ode23tb Cứng Thấp Sử dụng để ước lượng sai số khi giải các hệ cứng hoặc khi có ma trận khối lượng. Bảng 2-1: Các lệnh giải phương trình vi phân trong Matlab Tất cả các lệnh này đều có chung một cách sử dụng. Cú pháp đơn giản của chúng là: [t,Y] = solver('F',tspan,y0) Trong đó: Solver là tên lệnh ODE (ví dụ như ode45, ode15s ...) 'F' là tên của tệp m-file mô tả hệ phương trình vi phân y0 là véc tơ các điều kiện đầu (các điều kiện đầu phải đặt trong một véc tơ cột) tspan: xác định khoảng tích phân: tspan có thể là véc tơ có hai giá trị: giá trị đầu tiên quy định thời gian bắt đầu và giá trị thứ hai quy định thời gian kết thúc quá trình tính toán. Nếu tspan là một véc tơ nhiều phần tử, nó sẽ quy định các mốc thời gian cần xuất kết quả phân tích t là véc tơ các mốc thời gian xuất kết quả phân tích Y là một ma trận, trong đó mỗi cột của nó tương ứng là véc tơ các biến trạng thái, trật tự của các biến trạng thái được quy định trong m-file. Ví dụ: để giải hệ 2-4 trong khoảng thời gian từ 0 đến 10s, chúng ta cần gọi lệnh ODE như sau: [t,y]=ode45('Ptvp',[0, 10], [0; 1; 2]) Để vẽ ra đường y(t); y'(t); y"(t), chúng ta sẽ gọi các lệnh sau của matlab: plot(t, y(:,1)); plot(t,y(:,2); plot(t, y(:,3)) Chú ý: Nếu chúng ta gọi lệnh ODE không có vế trái, kết quả sẽ được xuất ra dưới dạng đồ thị các hàm y(t); y'(t); y"(t). 220 6.3.3. Giải hệ phương trình vi phân có ma trận khối lượng Việc giải hệ phương trình vi phân có ma trận khối lượng cũng tương tự như giải hệ ở dạng chuẩn có nghĩa là cũng cần thực hiện qua ba bước như trên. Các điểm khác biệt ở đây là:  Cách viết ODE file: Trong ODE file cần định nghĩa thêm ma trận khối lượng và bật chế độ giải hệ có ma trận khối lượng.  Lựa chọn lệnh ODE phù hợp (xem chi tiết ở bảng 2-1) Cách viết ODE file để giải hệ có ma trận khối lượng Để có thể định nghĩa được ma trận khối lượng cũng như bật chế độ giải hệ có ma trận khối lượng, cấu trúc của ODE file có dạng sau: function [out1, out2, out3] = odefile(t,y,flag) switch flag case '' % đưa ra hàm dy/dt = f(t,y). out1 = f(t,y); case 'init' % bật chế độ tính với ma trận khối lượng. out1=[ ]; out2=[ ]; out3=odeset('mass', 'giá trị'); case 'mass' % đưa ra ma trận khối lượng. out1= end; Chú ý: Khi flag='init' chúng ta sẽ gán thông số 'mass' của lệnh odeset một trong 3 giá trị sau:  'M' nếu ma trận khối lượng có các phần tử là hằng số,  'M(t)' nếu ma trận khối lượng có giá trị phụ thuộc vào thời gian  'M(t,y)' nếu ma trận khối lượng có giá trị phụ thuộc vào thời gian và biến trạng thái 6.3.4. Giải hệ phương trình vi phân có tham số Để đưa tham số vào phương trình vi phân chúng ta phải truyền tham số đó vào ODE file cũng như trong câu lệnh ODE chúng ta phải đưa giá trị tham số đó vào câu lệnh. Do trật tự các đối số trong câu lệnh ODE cũng như khi định nghĩa 221 ODE file rất chặt chẽ nên chúng ta phải đưa các tham số vào đúng vị trí của nó: a. Cấu trúc của tệp ODE function F=ODEFILE(t,y,flag, p1,p2,..) F=f(t,y,p1,p2..) % chèn hàm dy/dt với các thông số p1, p2.. vào đây b. Cú pháp câu lệnh ODE [t,Y] = solver('F',tspan,y0, [ ], val1, val2..) Trong đó:  val1, val2.. là giá trị tương ứng của thông số p1, p2..  Các đối số còn lại có ý nghĩa giống như phần trên.(chú ý đối số thứ 4 bắt buộc là giá trị rỗng [ ]) 6.3.5. Xác định các điểm cực trị, điểm uốn của đường cong nghiệm Tương tự như khi giải các hệ có ma trận khối lượng, để xác định các điểm đặc biệt như: cực trị, điểm uốn... của đường cong biến trạng thái chúng ta phải gọi lệnh Ode cũng như viết ode-file theo một cấu trúc khác: a. Viết Ode-file: Trong Ode-file chúng ta phải bật chế độ xác định các "sự kiện" và định nghĩa được các sự kiện cần xác định. Ode-file trong trường hợp này có cấu trúc như sau: function [out1, out2, out3] = odefile(t,y,flag) switch flag case '' % đưa ra các giá trị dy/dt = f(t,y). Out1 = f(t,y); case 'init' % bật chế độ xác định các sự kiện. Out1 = [ ]; Out2= [ ]; Out3=odeset('Events', 'on'); case 'events' % Định nghĩa các sự kiện. Out1 = Out2 = ; 222 Out3 = ; otherwise error(['Unknown flag ''' flag '''.']); end; Trong đó:  Véc tơ hàm sự kiện: là véc tơ cột các hàm phụ thuộc vào biến trạng thái và thời gian. Tuỳ thuộc vào yêu cầu xác định các sự kiện mà chúng ta chọn ra hàm sự kiện. Ví dụ như để xác định cực trị của hàm số, chúng ta cần xác định những điểm mà đạo hàm bậc nhất bằng 0 - có nghĩa là hàm y' được chọn là hàm sự kiện, còn nếu như xác định các điểm uốn thì hàm sự kiện sẽ là y".  Véc tơ lô gíc: là véc tơ cột các giá trị 0 và 1. Nếu giá trị lô gíc là 0, việc phân tích sẽ không dừng lại khi tìm thấy sự kiện, còn khi giá trị lôgic là 1, chương trình sẽ dừng lại ngay sau khi xác định được sự kiện.  Véc tơ hướng: Véc tơ hướng là một véc tơ cột các giá trị 0, -1, 1 dùng để xác định hướng phát triển của hàm số khi xác định sự kiện: giá trị 1, -1, 0 tương ứng chỉ ra rằng sự kiện sẽ được xác định khi hàm tăng, giảm, và không xét tới chiều biến thiên của hàm. Ví dụ như khi muốn xác định các điểm cực đại của hàm số chúng ta cần phải chọn hàm sự kiện là y', giá trị lô gíc =0, giá trị chỉ hướng bằng 1 b. Gọi lệnh Ode Trong trường hợp cần xác định các sự kiện, lệnh Ode có cú pháp sau: [t,Y,te,ye,ie] = solver('F',tspan,y0) Trong đó:  te là véc tơ các thời điểm xảy ra sự kiện.  ye là ma trận các biến trạng thái (mỗi cột của ma trậ là một biến) tương ứng với các thời điểm xảy ra sự kiện.  ie là véc tơ các chỉ số, chỉ ra sự kiện nào xảy ra ở thời điểm tương ứng. 6.3.6. Các thông số điều khiển chế độ tích phân Việc điều khiển hoạt động của các lệnh Ode được thực hiện thông qua lệnh odeset. cú pháp của lệnh này như sau: 223 Options = ODESET('Name1',Value1,'Name2',Value2,...) Trong đó:  'Name1', 'Name2'... Là tên các thuộc tính được quy định trong Matlab  Value1, Value2 ... Tương ứng là giá trị của các thuộc tính a. Các thuộc tính của ODESET Retol: Sai số tương đối, giá trị mặc định 1e-3. sai số trong mỗi bước tính cần thoả mãn điều kiện: e(i) <= max(RelTol*abs(y(i)),AbsTol(i)). AbsTol: Giá trị sai số tuyệt đối, giá trị mặc định 1e-6, giá trị này áp dụng cho tất cả các véc tơ tính toán. Refine: Hệ số độ mịn của dữ liệu xuất ra. Thuộc tính này sẽ tăng số lượng điểm của các dữ liệu ra. Giá trị mặc định bằng 1. OutputFcn: Tên của hàm dùng để xuất dữ liệu. Hàm xuất ra được gọi ra ở mỗi bước tích phân. Khi một lệnh ODE được gọi ra không có các đối số ra, máy sẽ dùng hàm quy định trong thuộc tính này, giá trị mặc định 'odeplot'. Giá trị của thuộc tính này là một chuỗi chỉ ra tên hàm cần gọi. OutputSel: Các chỉ số lựa chọn cho việc xuất kết quả. Véc tơ các chỉ số này sẽ chỉ ra hàm nào của véc tơ các hàm xuất ra sẽ được gửi tới hàm OutputFcn. Mặc định tất cả các hàm sẽ đưa ra. Giá trị của thuộc tính này là véc tơ các số nguyên. Stats: Hiển thị các thống kê của quá trình tính. Giá trị của nó là 'on' - bật chế độ hiển thị và 'off' - tắt chế độ hiển thị. Mặc định là 'off'. Jacobian - [on/off]: Hãy gán thuộc tính này là 'on' nếu tệp ODE được lập trình sao cho F(t,y,'Jacobian') đưa ra giá trị dF/dy. JConstant: Ma trận Jacobi hằng số [on/off]. Ấn định thuộc tính này là on nếu ma trận Jacobi là hằng số. JPattern - [on/off]: Hãy gán thuộc tính này on nếu tệp ODE được lập trình sao cho F([],[],'Jpattern') đưa ra ma trận thưa thớt. Vectorized - Véc tơ hoá tệp ODE [ on | {off} ]: Gán thuộc tính này là 'on' nếu tệp ODE file được lập trình sao cho F(t,[y1 y2 ...]) trả kết quả là véc tơ hàng [F(t,y1) F(t,y2) ...]. Events: Bật, tắt chế độ xác định các sự kiện. Gán thuộc tính này là 'on' nếu tệp ODE file được lập trình sao cho F(t,y,'events') trả các giá trị của các hàm event. Mặc định sẽ nhận giá trị 'off'. 224 Mass: Khai báo chế độ phân tích với ma trận khối lượng. giá trị của nó là [ {none} | M | M(t) | M(t,y) ]:  Chọn giá trị 'none' nếu tệp ODE được lập trình sao cho F(t,[],'mass') đưa ra một ma trận khối lượng.  'M' chỉ ra ma trận khối lượng là hằng số, 'M(t)' chỉ ra ma trận khối lượng phụ thuộc thời gian, 'M(t,y)' chỉ ra ma trận khối lượng phụ thuộc vào thời gian và các biến trạng thái. MassSingular: Ma trận khối lượng suy biến [ yes | no | {maybe} ]. Gán thuộc tính này là 'no' nếu ma trận khối lượng không suy biến. MaxStep: Là giá trị dương quy định cận trên của độ lớn bước tính. Mặc định giá trị bước tính bằng 1/10 bước trong TSPAN. InitialStep: Giá trị gợi ý của bước tính đầu tiên [ số dương ]. Mặc định máy sẽ định nghĩa chiều dài của bước tính đầu một cách tự động . MaxOrder: Bậc cao nhất của phép giải trong ODE15S [ 1 | 2 | 3 | 4 | {5} ] b. Phương thức sử dụng lệnh odeset Chúng ta có thể có hai cách khác nhau để sử dụng lệnh odeset: gọi lệnh odeset từ bên ngoài hay bên trong Ode-file: Gọi lệnh odeset từ môi trường làm việc của Matlab: Trường hợp này chúng ta sẽ gọi lệnh odeset trước khi gọi lệnh giải phương trình vi phân, sau đó đưa biến options (đã được xác định bằng lệnh odeset) vào trong lệnh Ode như sau: [t,Y] = solver('F',tspan,y0, options) Gọi lệnh odeset bên trong odefile Trong ode-file chúng ta gọi lệnh odeset khi cờ chỉ thị có giá trị là 'init' và chúng ta phải gán giá trị cho biến thứ 3 của hàm ode (xem các thí thí dụ ở trên). và khi chúng ta đã gọi được lệnh Odeset bên trong tệp Ode việc gọi lệnh Ode sẽ như trường hợp bình thường. Chú ý: Nếu chúng ta gọi lệnh odeset cả bên ngoài và bên trong, máy sẽ thực thi theo lệnh odeset gọi ở bên ngoài. 6.3.7. Các chú ý Ở phần trên chúng ta đã dẫn ra một số trường hợp điển hình của việc giải các hệ phương trình vi phân trong Matlab. Các trường hợp này có thể xuất hiện đơn lẻ, hoặc có thể phối hợp với nhau. Trong trường hợp phức tạp chúng ta phải 225 sử dụng cú pháp tổng quát để viết tệp Ode và gọi lệnh Ode như sau: a. Cú pháp tổng quát của lệnh Ode [T,Y,TE,YE,IE] = solver('F',tspan,y0,options) Các đối số F Tên của tệp ODE, đó là một hàm của MATLAB với các đối số là t và y, hàm này cần phải trả một véc tơ cột. Tất cả các lệnh ODE có thể giải phương trình vi phân ở dạng: Các hàm ode15s, ode23s, ode23t, and ode23tb có thể giải phương trình dạng : . Các hàm này ngoại trừ ode23s có thể giải phương trình dạng: tspan Là một véc tơ xác định khoảng tích phân [t0 tfinal]. Để nhận được kết quả tại các thời gian xác định (tăng dần hoặc giảm dần, hãy sử dụng tspan = [t0,t1, ..., tfinal]. y0 Véc tơ điều kiện đầu options Các đối số tích phân tuỳ chọn Optional (xem lệnh ODESET). p1,p2... Các thông số tuỳ chọn đưa vào hàm F t,Y Ma trận nghiệm Y, ở đây mỗi một hàng tương ứng với thời gian trong véc tơ t. b. Cấu trúc tổng quát của tệp Ode function varargout = odefile(t,y,flag,p1,p2) switch flag case '' % Return dy/dt = f(t,y). varargout{1} = f(t,y,p1,p2); case 'init' % Return default [tspan,y0,options]. [varargout{1:3}] = init(p1,p2); case 'jacobian' % Return Jacobian matrix df/dy. varargout{1} = jacobian(t,y,p1,p2); case 'jpattern' % Return sparsity pattern matrix S. varargout{1} = jpattern(t,y,p1,p2); 226 case 'mass ' % Return mass matrix. varargout{1} = mass(t,y,p1,p2); case 'events' % Return [value,isterminal,direction]. [varargout{1:3}] = events(t,y,p1,p2); otherwise error(['Unknown flag ''' flag '''.']); end; Các hàm được gọi trong ode- file được viết như sau: Hàm định nghĩa phương trình vi phân function dydt = f(t,y,p1,p2) dydt = Hàm xác định điều kiện đầu và các tuỳ chọn function [tspan,y0,options] = init(p1,p2) tspan = ; y0 = ; options = ; Hàm định nghĩa ma trận Jacobi function dfdy = jacobian(t,y,p1,p2) dfdy = ; Hàm định nghĩa ma trận khối lượng function M = mass(t,y,p1,p2) M = ; Hàm định nghĩa các sự kiện 227 function [value,isterminal,direction] = events(t,y,p1,p2) value = isterminal = ; direction = ; 6.3.8. Các ví dụ giải phương trình vi phân trong MATLAB Cho phương trình Van der Pol: Phương trình này có thể chuyển thành hệ phương trình: 1. Giải phương trình với : µ=1; t=[0 20]; y0=2; y'0=0 M-file function out1 = vdp1(t,y) out1 = [y(2); (1-y(1)^2)*y(2) - y(1)]; Câu lệnh gọi trong Matlab [t,y] = ode45('vdp1',[0 20],[2; 0]); plot(t,y(:,1),'-',t,y(:,2),'-.') 2. Giải phương trình trên với việc đưa điều kiện đầu trực tiếp vào ODE File M-file function [out1,out2,out3] = vdp1(t,y,flag) switch flag case '' out1 = [y(2); (1-y(1)^2)*y(2) - y(1)]; case 'init' 228 out1=[0 20]; out2=[2;0]; out3=[ ]; end; Câu lệnh gọi trong Matlab [t,y] = ode45('vdp1') plot(t,y(:,1),'-',t,y(:,2),'-.') 3. Giải phương trình vi phân trên với khả năng nhập giá trị µ từ dòng lệnh 3.1 Trường hợp lập trình đơn giản M-file function out1 = vdp1(t,y,flag,m) out1= [y(2); m*(1-y(1)^2)*y(2) - y(1)]; Câu lệnh gọi trong Matlab với µ=2 [t,y] = ode45('vdp1',[0 20],[2;0],[],2) plot(t,y(:,1),'-',t,y(:,2),'-.') Nếu gọi cách khác sẽ bị lỗi 3.2. Trường hợp lập trình mở rộng thêm tình huống M-file function out1 = vdp1(t,y,flag,m) if nargin<4 m=1; end; out1= [y(2); m*(1-y(1)^2)*y(2) - y(1)]; Câu lệnh gọi trong Matlab với µ=2 Nếu gọi trong Matlab bằng câu lệnh: 229 ode45('vdp1',[0, 20],[2;0],[ ], 2) Ta sẽ có nghiệm khi µ=2, Nếu gọi bằng câu lệnh: ode45('vdp1',[0, 20], [2;0]) Ta sẽ có nghiệm khi µ=1 (trường hợp mặc định) 3.3, Trường hợp mở rộng cho nhiều tình huống khác nhau: M-file function [out1,out2,out3] = vdp1(t,y,flag,m) if nargin<4 m=1; end; switch flag case '' out1 = [y(2); (1-y(1)^2)*y(2)*m - y(1)]; case 'init' out1=[0 20]; out2=[2;0]; out3=[ ]; end; Khi đó ta có thể gọi lệnh ODE bằng nhiều cách như: ode45('vdp1',[0 20],[2; 0]); ode45('vdp1'); ode45('vdp1',[0 20],[2;0],[],2) 4. Sử dụng ma trận khối lượng Phương trình có thể viết dưới dạng có ma trận khối lượng: 12 2 1 2 ' 2 ' 1 )1(10 01 yyy y y y x −−µ= M-file 230 function [out1,out2,out3] = vdp1(t,y,flag,m) if nargin<4 m=1; end; switch flag case '' out1 = [y(2); (1-y(1)^2)*y(2)*m - y(1)]; case 'init' out1=[0 20]; out2=[2;0]; out3=odeset('mass','M'); case 'mass' out1=[1 0; 0 1] end; Khi đó ta có thể gọi lệnh ODE bằng nhiều cách như: ode45('vdp1',[0 20],[2; 0]); ode45('vdp1'); ode45('vdp1',[0 20],[2;0],[],2) 5. Xác định các sự kiện đi qua 0 M-file function [out1,out2,out3] = vdp1(t,y,flag,m) if nargin<4 m=1; end; switch flag case '' out1 = [y(2); (1-y(1)^2)*y(2)*m - y(1)]; case 'init' out1=[0 20]; out2=[2;0]; out3=odeset('events','on'); 231 case 'events' out1=[ y(2); (1-y(1)^2)*y(2)*m - y(1)]; out2=[0;0];% không dừng lại khi gặp điểm 0 out3=[0;0]; % không xét đến hướng phát triển của hàm end; Câu lệnh gọi trong Matlab [t,y,te,ye,ie]= ode45('vdp1') Ở đây:  te là véc tơ các thời điểm mà các sự kiện(hàm đi qua điểm 0) xuất hiện;  ye là ma trận giá trị các hàm tại thời điểm sự kiện xuất hiện  ie là véc tơ các chỉ số xác định sự kiện nào sẽ xuất hiện (các giá trị của nó là 1 (có nghĩa là tại điểm đó hàm thứ nhất có giá trị 0) hoặc 2 đối với hàm thứ 2) Chính vì vậy để lấy các kết quả đưa ra ta dùng các dòng lệnh sau: e1=(ie==1); % trả giá trị đúng khi sự kiện thứ nhất xuất hiện e2=(ie==2); % trả giá trị đúng khi sự kiện thứ hai xuất hiện plot(t,y,te(e1),ye(e1,1),'ok',te(e2),ye(e2,2),'xk') xlabel('time'); title('Val Der pol Event Location') Kết quả thu được: 6. Thay đổi các tuỳ chọn để thay đổi các tuỳ chọn ta sẽ dùng câu lệnh odeset như sau: option=odeset('reltol',0.001,'abstol',0.0001) câu lệnh ode được viết : t,y] = ode45('vdp1',[0 20],[2; 0],option); plot(t,y(:,1),'-',t,y(:,2),'-.') 232 6.4. Bài tập ứng dụng BÀI 1: Đường đặc tính công suất của động cơ xăng được xác định theo công thức Lâydecman:        −   +   = 3 N e 2 N e N e xmaxee n n n n n n .NN Cho biết động cơ có số vòng quay cực tiểu là 1000 v/phút, tốc độ tại điểm có công suất lớn nhất là 3500 v/phút, hệ số λ=0.9, Công suất cực đại là 90KW. Dùng Matlab hãy xác định:  Véc tơ số vòng quay của động cơ ne (khoảng chia là 50 v/phút)  Véc tơ công suất ứng với số vòng quay  Véc tơ Mô men động cơ ứng với số vòng quay (Me=Ne/ωe) BÀI 2: Ô tô có đặc tính ngoài đưa ra bởi véc tơ công suất Ne và véc tơ số vòng quay ne. Hộp số có 5 tay số với các tỷ số truyền: i1 ..i5 = 5.7; 4.3; 3.5; 2.8; 1. tỷ số truyền của truyền lực chính io=5.3; hiệu suất truyền lực η=0.9; bán kính bánh xe r=0.48. Giả sử các véc tơ Ne, ne đã biết, sử dụng Matlab hãy xác định các véc tơ vận tốc V, Véc tơ lực kéo Pe ứng với từng tay số. Cho biết: công thức tính vận tốc của xe theo số vòng quay ở từng tay số là: ( ) 0hi e i.i.30 r.n iV π= Lực kéo được xác định bằng công thức: ηπ= .r.n. i.i.N30 Pe e hi0e )i( BÀI 3: Trên hình vẽ là mô hình của một cầu nhảy dưới dạng dầm bị ngàm một đầu, lực F thể hiện trọng lượng của người nhảy.  Hãy tạo ra một véc tơ F chứa các giá trị khác nhau của trọng lượng người nhảy (từ 300 N đến 750 N với gia số là 5N), một véc tơ X chứa khoảng cách a (từ 1000 đến 4000 mm, gia số là 75)  Tính ma trận giá trị mô men uốn tại ngàm 233  Tìm giá trị mô men uốn lớn nhất và bé nhất BÀI 4: Trên hình vẽ là mô hình dầm chịu lực:  Hãy tạo ra một véc tơ F chứa các giá trị khác nhau của lực (từ 200 N đến 750 N với gia số là 15N), một véc tơ X chứa khoảng cách A (từ 1000 đến 4000 mm, gia số là 30)  Tính ma trận giá trị mô men uốn tại điểm đặt lực  Tìm giá trị mô men uốn lớn nhất và bé nhất BÀI 5. TÍNH TOÁN ĐỘ CỨNG CỦA NHÍP: Độ cứng của nhíp được xác định từ biểu thức: ( )∑ = ++ − α= n 1k 1kk 3 1k YYa E6C Trong đó:  α là hệ số thực nghiệm: 0.85, E là mô đul đàn hồi của vật liệu  ak+1=l1- lk+1 , với lk là 1/2 chiều dài lá nhíp thứ k (tính từ quang treo đến đầu mút của lá)  Yk=1/Jk (Jk là mô men quán tính của mặt cắt ngang các lá nhíp) DỮ LIỆU ĐÃ BIẾT Chỉ số của lá Bề rộng b, Cm Bề dầy lá nhíp h, Cm Mô men quán tính J, Cm4 1-6 45 0,70 4,5.0,703/12=0,128 7-9 45 0,575 4,5.0,5753/12=0,0713 1/2 Chiều dài các lá nhíp (tính từ quang treo đến đầu mút của lá) l(1)=68.5;l(2)=67.2;l(3)=57.6;l(4)=50.4;l(5)=43.0;l(6)=35.0;l(7)=28.0;l(8)=20.6;l (9)=13.0 Dùng Matlab hãy tính độ cứng của nhíp với số liệu đã cho. BÀI 6. TÍNH TOÁN BỀN CÁC LÁ NHÍP 234 Phản lực tác động tại các đầu mút của lá nhíp được xác định từ hệ phương trình:    =+ =++ =++ − 0XBXA ............................... 0XCXBXA 0XCXBPA nn1nn 332313 22122 Trong đó Ak, Bk, Ck được xác định từ các biểu thức:     −= − − 1 l l 3 J J 5,0A k 1k 1k k k ,     +−= −1k k k J J 1B ,     −   = + + 1 l l .3. l l 5,0C 1k k 3 k 1k k ở đây:  lk là 1/2 chiều dài của lá nhíp thứ k (tính từ quang treo đến đầu mút của lá)  Jk là mô men quán tính mặt cắt ngang của lá nhíp thứ k : Hãy tính Tính ứng xuất tác dụng lên các lá nhíp có số liệu cho trong bảng dưới đây DỮ LIỆU CÁC LÁ NHÍP VÀ TẢI TRỌNG Số hiệu lá nhíp 1 2 3 4 5 6 Chiều dày lá nhíp 0.7 0.7 0.8 0.8 0.8 0.8 Bề rộng lá nhíp 6.5 6.5 6.5 6.5 6.5 6.5 Nửa chiều dài lá 50.0 48.5 37.5 28.0 20.0 13.0 Khoảng cách giữa bu lông ngàm nhíp: 12 cm Tải trọng đặt lên nhíp 900 Kg BÀI 7: Giả sử gia tốc của xe ở các tay số là các véc tơ J1, J2, J3, J4, J5 tương ứng với vận tốc là V1, V2, V3, V4, V5 (coi như đã biết), hãy lập chương trình tính và vẽ ra đồ thị thời gian và quãng đường tăng tốc phụ thuộc vào vận tốc của xe. biết rằng thời gian và quãng đường tăng tốc được xác định bằng công thức: dv J 1t v v0 ∫= ; ∫= t t 0 Vdts BÀI 8: Lập chương trình con tính ra giá trị của hai biểu thức phụ thuộc vào X:    >+ + ≤<−+ −≤++ = 2XKhi X1 1X 2X2Khi1X 2XKhi3X.2X Y 2 3 2 1 )X(cos1 )Xsin(1Y 22 + += 235 BÀI 9: Lập chương trình con tính giá trị của biểu thức phụ thuộc vào 2 biến X1, X2 như sau:        =+ >− + >− + = 12 2 1 12 12 2 1 21 21 2 1 XXKhiX1 XXKhi XX X1 XXKhi XX X1 Y  Tạo ra hai véc tơ X1 và X2 nằm trong khảng từ -10 đến 10, gia số bằng 0.5  Tính véc tơ Y theo 2 véc tơ X1 và X2 BÀI 10: Dưới đây là bảng biểu thị sự phụ thuộc của độ nhớt vào nhiệt độ, hãy đưa bảng số trên vào tệp có tên là kvisc.m và thực hiện các công việc sau: 1. Vẽ đồ thị biểu thị phụ thuộc độ nhớt vào nhiệt độ 2. Hãy xác định các giá trị của độ nhớt tại nhiệt độ 0,5 , 1.5, 2.5 ... 27.5oC 3. Tạo bảng với trật tự các biến ngược lại và tìm nhiệt độ tương ứng với độ nhớt 0,9 m/s2 Nhiệt độ (0C) Độ nhớt (m/s) Nhiệt độ (0C) Độ nhớt (m/s) 0 1.79 15 1.14 1 1.73 16 1.11 2 1.67 17 1.08 3 1.62 18 1.06 4 1.57 19 1.03 5 1.52 20 1.01 6 1.47 21 0.983 7 1.43 22 0.960 8 1.39 23 0.938 9 1.33 24 0.917 10 1.31 25 0.896 11 1.27 26 0.876 12 1.24 27 0.857 236 13 1.20 28 0.839 14 1.17 29 BÀI 11: HỆ SỐ TIẾT LƯU CỦA VAN XẢ Hệ số dòng chảy phụ thuộc chủ yếu vào đường kính lỗ xả d0 và chiều cao cột nước h. Các số liệu cho trên bảng.  Hãy tạo M-file để chứa dữ liệu nói trên  Sử dụng Matlab để xác định hệ số dòng chảy ứng với cột áp h=0.25 m và đường kính d=0.01, 0.02 ... 0.2 m Cột áp (h :m) Đường kính lỗ xả d0 (m) 0,39 0.18 0.06 0.03 0.015 0.006 0.21 0.590 0.594 0.601 0.611 0.622 0.651 0.24 0.591 0.594 0.601 0.610 0.620 0.648 0.27 0.591 0.595 0.601 0.609 0.618 0.646 0.30 0.591 0.595 0.600 0.608 0.617 0.644 0.40 0.593 0.596 0.600 0.605 0.613 0.638 BÀI 12: TÍNH TOÁN KÍCH THƯỚC KHỐI BÊ TÔNG Trên hình vẽ chỉ ra một khối hình trụ bao gồm một hình trụ bằng thép rỗng, hai đĩa thép và khối bê tông điền trong hình trụ. Khối lượng của toàn khối cần thiết kế là 45 kg, bề dày đĩa thép là t=5 mm, bề dày của hình trụ thép là h và đường kính trong là d (h=2d). hãy xác định đường kính d biết khối lượng riêng của thép là 7850 kg. m-3 và của bê tông là 2300 kg.m-3 237 BÀI 13: Xét khối bê tông- thép như hình vẽ trên. Hãy xác định đường kính d trong các trường hợp: bề dày =5, 5.5, 6...10 và khối lượng M=50, 75 ... 150 Đối với mỗi giá trị M hãy vẽ đường cong mô tả phụ thuộc của đường kính d vào bề dày t. BÀI 14: Dùng Matlab để giải các phương trình vi phân: a. (1+t2).y" + 2.t.Y' +3Y=0 với điều kiện đầu: t0=0; tf=5; Y(0)=0; Y'(0)=1 b. )( )( '" )( )(.'" tCosY tSin3 1YY t1 t2Cos5Y 2 =++++− với t0=0; t1=5 và điều kiện đầu: Y0=1; Y'0=0; Y"0=2 BÀI 15: Dùng Matlab để giải các phương trình vi phân và tìm điểm cực đại của hàm tìm được a. W' + (1.2 + sin10t)W=0; t0 = 0, tn=5, W(t0)=1 b. 3W' + W t1 1 2+ =cos t, t0 = 0, tn=5, W(t0)=1 c. (1 +t2)W" +2tW' +3W =2; t0=0, tn=5, W'(t0)=1, W(t0)=0 BÀI 16: Giải hệ phương trình:   =β+α =++− 1YX 0YCos1eXSin 3x1 .. )().( Với các bộ giá trị sau của α và β: (1, 1); (1, 2); (2, 1) BÀI 17: Tìm nghiệm của hệ hai phương trình 2 ẩn BÀI 18: Tìm cực tiểu của hàm số: Y=X12+ 2.X22 với các ràng buộc:    ≥ ≤− ≥+ 0X 1XX 2XX.2 1 21 21 BÀI 19: Tìm cực tiểu của hàm số: Y=3.X1 + 4.X2 với các ràng buộc ( ) ( )   ≥+− ≤++ 3X1X 6X1X 2 2 2 1 2 2 2 1 238 BÀI 20: Tìm cực tiểu của hàm số: f(x)=-X1.X2.X3, với các ràng buộc