Đề tài Matlab tìm hiểu và ứng dụng giải một số bài toán kĩ thuật

phép tính. Ví dụ: X.^Y, X.*Y, hay X.\Y. Nếu X và Y là các ma trận ( hay véc tơ) các phần tử của X sẽ đ−ợc nâng lên luỹ thừa hoặc đ−ợc nhân , chia bởi các phần tử t−ơng ứng của ma trận Y. Dẽ thấy là ma trận X và Y phải cùng cỡ. - Luỹ thừa các phần tử ma trận Ví dụ: >> x=[1 2 ; 3 4] x = 1 2 3 4 >> y=[ 3 4; 1 2] y = 3 4 1 2 >> x.^y ans = 1 16 3 16 x = 1 2 3 4 Nếu y không phải là ma trận mà là một số, ví dụ y=2, thì kết quả nh− sau: >> x.^2 ans = 1 4 9 16 - Nhân phần tử ma trận Ví dụ >> x.*y ans = 3 8 2 8 - Chia phần tử ma trận Ví dụ; >> x./y ans = 0.3333 0.5000 3.0000 2.0000 6. Các hàm thực hiện các phép tính với đa thức - Hàm poly : Xác định đa thức khi biết tr−ớc nghiệm Quy các viết p=poly(b) trong đó b là một ma trận hàng. Kết quả sẽ cho ra là một ma trận hàng mà mỗi phần tử của nó là một hệ số của một đa thức có nghiệm là các phần tử của ma trận b ( theo số mũ giảm dần). Ví dụ: b=[2 1 -4 3]; p=poly(b) p = 1 -2 -13 38 -24 Theo kết quả trên thì các số 2, 1, -4 và 3 là nghiệm của đa thức : x4- 2x 3 - 13x2 + 38x - 24 = 0 - Hàm roots : Xác định nghiệm của đa thức Quy cách viết : a=roots(b) trong đó b là ma trận hàng với các phần tử là các hệ số của đa thức (theo số mũ giảm dần). Kết quả cho ra là một ma trận cột mà các phần tử là nghiệm của đa thức. Ví dụ: Thử tìm lại nghiệm của đa thức trên. Ta viết các lện nh− sau: p=[ 1 -2 -13 38 -24 ]; r=roots(p) r = -4.0000 3.0000 2.0000 1.0000 Ta thấy kết quả hoàn toàn chính xác. - Hàm conv : Dùng nhân đa thức. Quy cách viết: a=conv(b,c) trong đó b,c là hai ma trận hàng có các phần tử là các hệ số của các đa thức cần nhân. Kết quả cho ra là ma trận a có các phần tử là hệ số của ma trận tích. Ví dụ : cần nhân hai đa thức x3+2x2+6 và 3x4-6x2+5x-10 ta làm nh− sau: b=[1 2 0 6]; c=[3 0 -6 5 -10]; a=conv(b,c) a = 3 6 -6 11 0 -56 30 -60 Vậy đa thức tích là: 3x7+ 6x6- 6x5 + 11x4 - 56x2 + 30x – 60 - Hàm deconv : Dùng chia hai đa thức Qui cách viết nh− sau: [ m , n] = deconv(p,q) với p và q là hai ma trận hàng có các phần tử là các hệ số của đa thức bị chia và đa thức chia, còn các phần tử của ma trận m,n là các hệ số của đa thức th−ơng và phần d−. Ví dụ : Ta thử chia ngay đa thức tích vừa có ở trên cho đa thức có các hệ số là các phần tử của ma trận c, tức đa thức: 3x4-6x2+5x-10. a=[3 6 -6 11 0 -56 30 -60]; c=[3 0 -6 5 -10]; [b,d]=deconv(a,c) b = 1 2 0 6 d = 0 0 0 0 0 0 0 0 Ta thấy kết quả là hoàn toàn đúng. 7. các Hàm dùng phân tích dữ liệu Các hàm tìm giá trị cực đại, cực tiểu và trung bình - Hàm max : Tìm giá trị lớn nhất Qui cách viết ln=max(a) với a là ma trận hàng. Kết quả cho ra là một phần tử có giá trị lớn nhất của a Ví dụ: a=[ 10 2 1 -30 23 8]; ln=max(a) ln = 23 - Hàm min : tìm giá trị cực tiểu Ví dụ: bn=min(a) bn = -30 bn=mi - Hàm mean : Tìm giá trị trung bình Ví dụ: tb=mean(a) tb = 2.3333 Nếu a là một ma trận có nhiều hàng nhiều cột thì các giá trị max, min hoặc trung bình sẽ là các giá trị t−ơng ứng đối với các cột của ma trận. Ví du : >>b=[1 2 3; 4 5 6 ; 7 8 9] b = 1 2 3 4 5 6 7 8 9 >>ln=max(b) ln = 7 8 9 >>bn=min(b) bn = 1 2 3 >>tb=mean(b) tb = 4 5 6 - Hàm sum : dùng tính tổng Ví dụ: Tính tổng các phần tử của ma trận a nh− sau T=sum(a) T = 26 - Hàm diff : Tính giá trị sai khác của hai số đứng liền nhau. Qui cách viết: s=diff(x) với x là một ma trận hàng hoặc cột. Ví dụ: >>x=[ 1.2 1.4 1.8 2.1 3 ]; >> s=diff(x) s = 0.2000 0.4000 0.3000 0.9000 Ta dễ dàng thấy rằng hàm diff này có thể dùng để tính gần đúng đạo hàm. Nếu x là một ma trận bình th−ờng thì quá trình tính sẽ đ−ợc thực hiện theo thứ tự các cột. - Hàm Interp1 : Dùng tìm các giá trị bị khuyết. Ví dụ: Đã biết giá trị của hàm y=x3-3x+4 tại các điểm có x=0,1,2,3,4 và 5. Hãy xác định giá trị của y tại các điểm có x= 0,3, 0,5, ..., 4.2, 4,8. Ta viết các lệnh nh− sau: x1=0:5; y1=x1.^3-3*x1+4; x2=[0 .3 .5 1.2 1.4 2.1 3.3 4.2 4.8 5]; y2=interp1(x1,y1,x2,'cubic') plot(x1,y1,'ro',x2,y2,'b+') Kết quả cho ở dạng số và đồ thị d−ới đây y2 = Columns 1 through 7 4.0000 2.7700 2.2500 2.2240 2.5920 7.0330 30.1210 Columns 8 through 10 65.6800 100.4800 114.0000 Đó là 10 giá trị của y ứng với 10 giá trị của x ( trong ch−ơng trình tính là x2). Còn đồ thị sẽ cho thấy sự ppù hợp của các kết quả này: Các điểm vẽ bằng dấu ‘+’ là biểu thị các điểm cần tìm, còn các điểm vẽ bằng dấu ‘o’ là t−ơng ứng các điểm đã cho. 0 1 2 3 4 5 0 2 0 4 0 6 0 8 0 1 0 0 1 2 0 8. Hàm của hàm Có nhiều hàm trong Matlab không chỉ làm việc với các ma trận số mà còn làm việc với các hàm toán. Các hàm của hàm này bao gồm các hàm dùng để tính tích phân, giải các ph−ơng trình phi tuyến và giải các ph−ơng trình vi phân. a. Hàm dùng để tích phân số: hàm quad hoặc quad8 Ví dụ: Cần tính tích phân của hàm f(x) = 3+sin2x/(4+2cosx) với cận từ 0 đến 4, ta tiến hành nh− sau: Viết một file hàm có tên tfan1.m nh− d−ới đây: function f=tfan1(x) f=3+sin(2*x)/(4+2*cos(x)); và một file chính có tên tfan.m: tf=quad8('tfan1',0,4) Chạy ch−ơng trình kết quả cho giá trị của tích phân: tf = 12.0517 b. Hàm để tìm nghiệm ph−ơng trình phi tuyến và các hàm tối −u - Hàm fzero : xác định nghiệm của hàm một biến. - Hàm fmin : tính giá trị cực tiểu của hàm một biến. - Hàm fmins : tính các giá trị cực tiểu của hàm nhiều biến Ví dụ: Dùng hàm fzero tìm nghiệm của ph−ơng trình sau tại vị trí gần x=1 : y=2- 6sinx / (1+x); Ta viết file hàm nghiem1.m nh− sau: function y=nghiem(x) y=2-6*sin(x)./(1+x); Sau đó viết file chính nghiem.m: x1=fzero('nghiem1',1) Chạy ch−ơng trình chính, kết quả cho trên màn hình nh− sau: Zero found in the interval: [0.36, 1.4525]. x1 = 0.5385 Để thử lại, ta thay x=x1 vào ph−ơng trình ban đầu để xem kết quả hai vế có bằng 0 hay không. y=2-6*sin(.5385)/(1+.5385) y = -6.0516e-005 Ta thấy y~ 0 và có thể nói kết quả đủ chính xác. c. Hàm để giải ph−ơng trình vi phân : hàm ode23 và ode34 Ví dụ: Giải ph−ơng trình vi phân Van de Pol có dạng nh− sau: x’’+( x2- 1 ) x’+ x=0 Ta viết lại ph−ơng trình trên d−ới dạng hệ ph−ơng trình vi phân bậc nhất: x1’= x1(1- x2 2 ) - x2 x2’= x1 Viết một file hàm mang tên ftvf.m mô tả hệ ph−ơng trình vi phân: function xc=ftvf(t,x) xc=zeros(2,1); xc(1)=x(1).*(1-x(2).^2)-x(2); xc(2)=x(1); Viết file chính có tên ftvf0.m nh− d−ới đây: t0=0; t1=20;% Thoi diem dau va cuoi x0=[0 .2]; % Dieu kien ban dau [t,x]=ode23('ftvf',t0,t1,x0); plot(t,x) Chạy ch−ơng trình kết quả cho d−ới dạng đồ thị nh− d−ới đây: 0 5 1 0 1 5 2 0 -3 -2 -1 0 1 2 3 9. Các toán tử điều khiển Cũng nh− các phần mềm lập trình khác, các toán tử điều khiển cũng có trong Matlab. Đó là toán tử vòng lặp for, while và toán tử lựa chọn if. - Toán tử for Dạng chung của toán tử này đ−ợc viết nh− sau: for biến = biểu thức các câu lệnh end Biểu thức ở đây th−ờng có dạng m:n hoặc m:i:n, trong đó m, n là giá trị đầu và cuối, còn i là gia số. Ví dụ: for k=1: n a(k)=sin(k*pi/5); b(k)=có(k*pi/5); end - Toán tử while Dạng chung của lệnh vòng lặp này có dạng sau: while biểu thức các câu lệnh end Biểu thức ở đây là biểu thức quan hệ. Ví dụ: d=1; while d>0.001 z1=z2-sin(z2)/(z2+2); d=abs(z2-z1); z2=z1; end - Toán tử điều kiện if Toán tử có dạng chung nh− sau: if biểu thức các câu lệnh elseif biểu thức các câu lệnh ... ... else các câu lệnh end Cũng nh− đối với toán tử while biểu thức ở đây cũng là biểu thức quan hệ. Ví dụ: for k=1: n for p=1: m if k= =p z(k,p)=1; elseif k<p z(k,p)=-1 ; else z(k,p)=0; end end end - Câu lệnh Break Lệnh Break cho phép thoát ra khỏi vòng lặp. Ví dụ: k=input(' k='); if k>0 break else a=5*k+4 end 10. các loại file dùng trong matlab Có hai loại file là M-File và file dữ liệu (data file) Các file có chứa các lệnh của Matlab đ−ợc gọi là M- file . Sở dĩ gọi là M- file vì phần mở rộng của các file này là .m. Có hai loại M- file là script file và function file (file hàm) Script File Khi Script đ−ợc kích hoạt Matlab đơn giản thực hiện các lệnh tìm thấy trong file. Các câu lệnh trong file này thì hoạt động trên toàn bộ các dữ liệu trong Workspace. Ví dụ : file lg.m d−ới đây là một script file: % File dung tinh ham luong giac va ve do thi f1=0:2*pi/60; a=3*sin(f1)+4.5*cos(f1+.6); plot(f1,a) Sau khi tính toán và vẽ đò thị xong các biến f1 và a vẫn còn l−u trong workspace Function File Là một M-File có chứa từ function tại vị trí đầu tiên trên dòng đầu của file. Với function file, các biến đ−ợc định nghĩa và hoạt động chỉ trong phạm vi file, chứ không có tính toàn cục nh− trong script file. Function file rất tiện ích trong việc mở rộng thêm khả năng của Matlab, cụ thể là tạo thêm đ−ợc các Matlab file mới. D−ới đây là một ví dụ về function file đ−ợc tạo để giải ph−ơng trình bậc hai, tên file là ftb2.m : function [x1 , x2] = ftb2(a,b,c) % File ham nay dung de giai Fuong trình bac hai : ax^2+bx+c=0 % voi a, b, c la cac hang so cho truoc. Cac nghiem là x1 va x2 delta=b^2-4*a*c; x1=(-b+sqrt(delta)/(2*a); x2=(-b-sqrt(delta)/(2*a); Hàm này dùng để giải ph−ơng trình bậc hai với lệnh gọi ra là: [nghiem1, nghiem2] = ftb2(10,-6.5,-3.5); File dữ liệu ( data file) Ta có thể ghi số liệu vào một file gọi là file số liệu với phần mở rộng là .dat. Khi cần file có thể đ−ợc gọi ra. D−ới đây là ví dụ về tạo và gọi file số liệu: Giả sử ta có một ma trận hàng a nh− sau: >>a=1:8 a = 1 2 3 4 5 6 7 8 Giờ ta muốn ghi các giá trị của ma trận a vào một file dữ liệu có tên là sl1.dat, ta dùng lệnh save với qui cách viết nh− sau: >> save sl1.dat a -ascii Khi đó ma trận a đã đ−ợc l−u trong fil sl1.dat. Nếu muốn gọi ra các số liệu này và gán cho cho nó một tên biến mới, ví dụ là x, ta dùng lệnh load với qui cách viết nh− sau: >>load sl1.dat >> x=sl1 x = 1 2 3 4 5 6 7 8 11. Xử lí tín hiệu Matlab có một số hàm dùng để xử lí tín hiệu số d−ới đây: abs - giá trị (mô đun) của số phức angle – góc pha conv - tích chập ( convolution) deconv – ( deconvotution) fft - biến đổi Fu ri ê nhanh ift - biến đổi ng−ợc Fu ri ê fftshift Ví dụ: cho x là một ma trận cột nh− sau: x=[ 4 3 7 -9 1 0 0 0 ]’ ; Biến đổi Fu ri ê nhanh của x đ−ợc viết nh− sau: y=fft(x) Kết quả cho nh− d−ới đây: y = 6.0000 11.4853 + 2.7574i -2.0000 +12.0000i -5.4853 -11.2426i 18.0000 -5.4853 +11.2426i -2.0000 -12.0000i 11.4853 - 2.7574i 12. vẽ đồ thị Matlab rất mạnh với việc vẽ đồ thị. D−ới đây trình bày cách vẽ đồ thị trong không gian hai chiều (2D) và ba chiều (3D). a.Vẽ đồ thị 2D Vẽ các đ−ờng cong thông th−ờng: Dùng lệnh plot với qui cách viết nh− sau: plot(x,y) trong đó x và y là các ma trận hàng (hoặc cột ) cùng cỡ. Matlab sẽ vẽ ra một đ−ờng cong của y theo x. Nếu bạn muốn chọn màu và kiểu nét vẽ cho đ−ơng cong có thể làm theo ví dụ sau: plot(x,y, 'r*') Khi chạy, Matab sẽ cho bạn đ−ờng cong tạo bởi các dấu *** và có màu đỏ (chữ r đứng cho từ red trong tiếng Anh, nghĩa là đỏ). D−ới đây là bảng kí hiệu màu và kiểu nét vẽ dùng khi vẽ đồ thị: Biểu t−ợng Màu sắc Biểu t−ợng Loại đ−ờng nét y m c r g b w k Vàng Đỏ t−ơi Xanh lá mạ Đỏ Xanh lá cây Xanh n−ớc biển Trắng Đen . ° ì + ∗ - : -- -. Điểm Vòng tròn Dấu nhân Dấu cộng Dấu sao Nét liền Nét chấm Nét đứt Nét chấm gạch Cũng có thể vẽ nhiều đ−ờng cong trên cùng một hệ trục toạ độ. Ví dụ: lệnh vẽ sau sẽ cho 3 đ−ờng cong khác nhau trên một hình: plot(x1,y1,’r *’,x2,y2,’g+’x3,y3,’w -‘) Nếu bạn muốn đặt tên cho đồ thị hãy dùng lệnh title. Ví dụ: title(' Do thi 1 ‘) Bạn cũng có thể đặt tên cho các trục toạ độ x và y bằng cách dùng lệnh xlabel và ylabel . Ví dụ: xlabel('x - m'); ylabel('y - m/s '); Còn để tạo l−ới trên hình vẽ dùng lệnh grid. Cụ thể ta làm một ví dụ d−ới đây: Vẽ trên cùng hệ trục toạ độ đồ thị hai hàm số sau y1=e-5tsin5t ; y2= e-5t cos5t; với biến t chạy từ - 4 đến +4: Ta viêt file p2d.m để vẽ đồ thị nh− sau: % Chuong trinh ve do thi 2D t=-4:.05:4; y1=exp(-.5*t).*sin(5*t); y2=exp(-.5*t).*cos(5*t); plot(t,y1,'*',t,y2,'-'); title('VE DO THI 2 CHIEU : y1( ***), y2 (---)') xlabel(' Gia tri t ') ylabel('gia tri y1, y2') grid Chạy ch−ơng trinh cho kết quả là đồ thị d−ới đây: -4 -3 -2 -1 0 1 2 3 4 -8 -6 -4 -2 0 2 4 6 8 V E D O T H I 2 C H IE U : y 1 ( * * * ) , y 2 ( - - - ) G i a tr i t gi a tri y 1, y 2 Có thể vẽ nhiều đồ thị trên một hình bằng cách dùng lệnh subplot. Ví dụ: ta vẽ hai đồ thị y1=y1(t) và y2=y2(t) vừa nêu trên cùng một hình nh−ng lại trên hai hệ trục riêng biệt, ta viết lại các các lệnh ( tính từ lệnh plot(t,y1,...) nh− sau: subplot(211) plot(t,y1,'*') title('VE DO THI 2 CHIEU : y1( ***), y2 (---)') xlabel('gia tri t') ylabel('gia tri y1') grid subplot(212) plot(t,y2,'-') xlabel('gia tri t') ylabel('gia tri y2') grid -4 -3 -2 -1 0 1 2 3 4 -1 0 -5 0 5 1 0 V E D O TH I 2 C H IE U : y1 ( * ** ), y2 (---) g ia tr i x gi a tri y 1 -4 -3 -2 -1 0 1 2 3 4 -5 0 5 1 0 g ia tr i x gi a tri y 2 Vẽ đồ thị độc cực Dùng lệnh polar Ví dụ: vẽ đồ thị hàm y=4e-4t với t=0:2pi/3 Ta viết ch−ơng trình vẽ nh− sau: t=0:pi/20:2*pi/3; r=4*exp(-4*t); polar(t,r) Kết quả cho ra nh− d−ới đây: 1 2 3 4 30 210 60 240 90 270 120 300 150 330 180 0 Vẽ Đồ thị dạng thanh (Bar). Đồ thị bar biểu diễn giá trị của một véc tơ hay một ma trận d−ới dạng các thanh thẳng đứng hay nằm ngang. Bar(y) biểu diễn các phần tử của véc tơ y nếu y là một véc tơ. Nếu y là một ma trận, một nhóm thanh sẽ biểu diễn các phần tử trong một hàng, trục x sẽ biểu diễn số hàng. Ví dụ: Vẽ đồ thị biểu thị số sinh viên tốt nghiêp của một tr−ờng Đại học từ năm 1995 đến 2001, biết số sinh viên tốt nghiệp t−ơng ứng là 1000, 1200, 1800, 1850, 1750, 2100. Lập ch−ơng trình vẽ nh− sau: x=[1000, 1200, 1800, 1850, 1750, 2100] bar(x); title('SO LUONG SINH VIEN TOT NGHIEP HANG NAM') xlabel('Nam') ylabel('So luong sinh vien') set(gca,'XTickLabel',{'1996';'1997';'1998';... '1999';'2000';'2001'}) Chạy ch−ơng trình cho kết quả là đồ thị d−ới đây: 1 9 9 6 1 9 9 7 1 9 9 8 1 9 9 9 2 0 0 0 2 0 0 1 0 5 0 0 1 0 0 0 1 5 0 0 2 0 0 0 2 5 0 0 S O L U O N G S IN H V IE N T O T N G H IE P H A N G N A M N a m S o lu on g si nh v ie n Đồ thị pie. Lệnh Pie(x) biểu thị thành phần phần trăm theo dữ liệu trong x. Mỗi phần tử trong x đ−ợc đại diện bằng một phần của đồ thị Pie. Pie(x,explode) đặt các phần tử cho đồ thị Pie. Explode là một véc tơ hay một ma trận của số 0 và 1 t−ơng ứng với x. Các phần tử là 0 sẽ t−ơng ứng với các phần tạo thành một khối liền trong đồ thị pie. Nếu là giá trị 1 thì phần tử t−ơng ứng với mẩu trên đồ thị sẽ đ−ợc tách rời khỏi khối liền trong đồ thị Pie. Explode phải có cỡ giống với x. Ví dụ : Biểu diễn thành phần phần trăm của các loại vật liệu có trong một hỗn hợp. x=[7 15 56 12 10]; eplode= [1 0 0 0 0]; pie(x,eplode) colormap jet Ta đ−ợc đồ thị : 7 % 1 5 % 5 6 % 1 2 % 1 0 % b. Vẽ đồ thị 3D : dùng lệnh surf, mesh, contour3 Ví dụ: Vẽ đồ thị z=-20x2+x-15y2+5y , với -4<x<4 và -4<y<4. Viêt file p3d.m để vẽ đồ thị ( dùng lệnh surf )nh− sau: [x,y]=meshgrid(-4:.2:4,-4:.2:4); z=-20*x.^2+x-15*y.^2+5*y; surf(x,y,z) title(' Do thi 3D') xlabel(' truc x') ylabel(' truc y') zlabel(' truc z') pause Kết quả cho nh− hình d−ới đây: Còn hai đồ thị d−ới đây đ−ợc vẽ bằng lện mesh và lệnh contour3: - Đồ thị 3D dùng lệnh mesh - Đồ thị 3D dùng lệnh contour 1 0 2 0 3 0 4 0 1 0 2 0 3 0 4 0 -5 0 0 -4 0 0 -3 0 0 -2 0 0 -1 0 0 0 truc x truc y tr uc z D o th i 3 D d ung le nh c o nto ur Phần 2 ứng dụng Matlab giải một số bài toấn kĩ thuật 1. Bài toán về mạch điện Ví dụ 1: Tính c−ờng độ dòng điện của mạch điện có các trở kháng phức d−ới đây, biết hiệu điện thế nguồn U=5V, các trở kháng: z1=1+0,5i, z2=2+0,5i, z3=4+i, z4=z1 (Ω ). z1 z1 z1 z1 aaa U=5v fff z4 z4 z4 z4 Giải: Viết các ph−ơng trình định luật Kiếc sốp cho mạch nh− sau: (z1+z2+z4)I1- z2I2 = U (z1+2z2+z4)I2 - z2I1- z2I3= 0 (z1+2z2+z4)I3- z2I2- z2I4= 0 (z1+z2+z3+z4)I4- z2I3 = 0 Viết lại hệ ph−ơng trình trên d−ới dạng ma trận nh− sau: (z1+z2+z4) - z2 0 0 I1 = U - z2 (z1+2z2+z4)I2 - z2 0 I2 = 0 0 - z2 (z1+2z2+z4) - z2 I3 = 0 0 0 - z2 (z1+z2+z3+z4) I4 = 0 hay Z.I = U với Z là ma trân trở kháng, I là ma trận cột dòng điện, còn U là ma trận cột thế hiệu. Nghiệm I đ−ợc tính nh− sau: I=Z-1.U D−ới đây là file ml12.m để giải ph−ơng trình ma trận trên % File ml12.m % Chuong trinh giai mach dien tro khang phuc z1=1+.5i; z2=2+.5i; z3=4+i; z4=z1; % Cac tro khang U=5; % Hieu dien the nguon % Ma tran tro khang Z=[z1+z2+z4 -z2 0 0;... -z2 z1+2*z2+z4 -z2 0;... 0 -z2 z1+2*z2+z4 -z2;... 0 0 -z2 z1+z2+z3+z4 ]; F=[ V 0 0 0 ]'; % Ma tran cot ve phai I=inv(Z)*F % Tinh dong dien Chạy ch−ơng trình cho kết quả d−ới đây: I = 1.3008 - 0.5560i 0.4560 - 0.2504i 0.1530 - 0.1026i 0.0361 - 0.0274i Ví dụ 2: Tính dòng điện cho mạch có R, L và C nh− hình vẽ, cho E=60v, R=50 , L=1H, Ω C=10-4F, tại thời điểm ban đầu i1(0)=0, i2(0)=0. i1 L i3 E oi R C Viết hệ ph−ơng trình vi phân cho mạch nh− sau: Ldi1/dt + Ri2 = E RCdi2/dt+ i2 - i1 = 0 Chuyển vế, viết lại hệ ph−ơng trình trên nh− d−ới đây: i1’ = - Ri2/L + E/L (1) i2’ = i1/RC - i2 /RC Dùng hàm ode45, lập ch−ơng trình để xác định i1, i2 bằng các file ml13.m và ml13a.m d−ới đây: % File ml13.m CHUONG TRINHF TINH MACH DIEN % CO R-L-C global E R L C E=60; L=1; R=50;C=.005 % Nhap cac so lieu i0=[0 0]'; % Dieu kien ban dau t=5; % thoi gian [t,i]=ode45('ml13a',t,i0); subplot(211) plot(t,i(:,1)) title('GIAI MACH DIEN CO R- L- C') ylabel('D.dien trong mach chinh I1 - A') grid subplot(212) plot(t,i(:,2)) ylabel('D.dien trong mach re I2 - A') xlabel(' Thoi gian') grid % File ml13a.m function ic=ml13a(t,i) global E R L C ic=zeros(2,1); ic(1)=-R/L*i(2)+E/L; ic(2)=i(1)/(R*C)-i(2)/(R*C); Chạy ch−ơng trình cho kết quả là trị số c−ờng độ dòng điện trong mạch chíhh i1 và mạch rẽ i2 nh− hình d−ới đây: 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 0 0.5 1 1.5 G IA I M A CH D IE N CO R- L- C D .d ie n t ro ng m ac h ch in h I1 - A 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 0 0.5 1 1.5 D .d ie n t ro ng m ac h re I2 - A Tho i g ian Để kiểm tra kết quả trên, ta so sánh nó với kết quả chính xác của hệ ph−ơng trình vi phân (1) là: i1= 6/5- 6e -100t/5- 60te-100t i2= 6/5- 6e -100t/5- 120te-100t Đồ thị chính xác của i1 và i2 có dạng nh− hình d−ới đây: 0 0 . 5 1 1 . 5 2 2 . 5 3 3 . 5 4 4 . 5 5 0 0 . 2 0 . 4 0 . 6 0 . 8 1 1 . 2 1 . 4 D o n g d i e n I 1 = 6 / 5 - 6 / 5 e ( - 1 0 0 t ) - 6 0 t e ( - 1 0 0 t ) 0 0 . 5 1 1 . 5 2 2 . 5 3 3 . 5 4 4 . 5 5 0 0 . 2 0 . 4 0 . 6 0 . 8 1 1 . 2 1 . 4 D o n g d i e n I 2 = 6 / 5 - 6 / 5 e ( - 1 0 0 t ) - 1 2 0 t e ( - 1 0 0 t ) Ta thấy kết quả của hai cách tính là nh− nhau. Để thấy sự dao động của mạch điện một cách rõ nét ta tăng điện dung của tụ điên: lấy C=0,005F rồi chạy lại ch−ơng trình ta đ−ợc kết quả nh− hình d−ới đây: 0 0 .5 1 1 . 5 2 2 . 5 3 3 . 5 4 4 . 5 5 -2 0 2 4 6 G IA I M A C H D IE N C O R - L - C D on g di en t ro ng m ac h ch in h - A 0 0 . 5 1 1 . 5 2 2 . 5 3 3 . 5 4 4 . 5 5 0 0 .5 1 1 .5 2 D on g di en t ro ng m ac h re R - A Th o i g ia n 2. Giải bài toán động học cơ cấu phẳng Ví dụ : Giải bài toán động học cơ cấu cu lit Hãy xác định chuyển vị, vận tốc góc, gia tốc góc của thanh cu lit 3 của cơ cấu cu lit, biết tr−ớc kích th−ớc các khâu: l0=0,5m, l1=0,2m, tay quay 1 quay đều với w1=8rad/s. l1 1ϕ l0 l3 3ϕ Ta có quan hệ véc tơ sau: l3 = l0 + l1 Viết l0 , l1 d−ới dạng số phức ta đ−ợc: l3 = l0e i 2/π + l1e i 1ϕ (1) Dùng hàm abs ta tính đ−ợc trị số của véc tơ l3 , dùng hàm angle để tính góc chỉ ph−ơng của véc tơ này cũng là góc 3ϕ của cu lit. l3a = abs( l3 ) (2) 3ϕ = angle( l3 ) (3) Để tính đ−ợc vận tốc góc và gia tốc góc của cu lít ta thay l3 = l3aei 3ϕ vào (1) rồi đạo hàm theo t, ta đ−ợc: l3ae i 3ϕ = l0ei 2/π + l1ei 1ϕ l3a’e i 3ϕ + i 3ω l3aei 3ϕ = i 1ω l1ei 1ϕ Chia cả hai vế cho ei 3ϕ : l3a’ + i 3ω l3a = i 1ω l1ei )( 31 ϕϕ − Tách phần thực và ảo, ta rút ra: l3a’ = - 1ω l1sin( )( 31 ϕϕ − (4) 3ω = 1ω l1cos )( 31 ϕϕ − / l3a (5) Để có gia tốc góc, ta đạo hàm (5) : α 3 = ( - 1ω l1sin )( 31 ϕϕ − ( 1ω - 3ω ) l3a - l3a’ 3ω )/ l3a (6) Căn cứ vào các biểu thức (1), (2), (3), (4), (5) và (6) lập ch−ơng trình tính các chuyển vị vận tốc và gia tốc nh− đ−ợc viết trong file cl3.m d−ới đây. Kết quả về chuyển vị, vận tốc và gia tốc góc cu lít đ−ợc cho d−ới dạng đồ thị ở phía d−ới ch−ơng trình. % CHUONG TRINH TINH CHUYEN DONG CUA CO CAU CU LIT l1=.2; l0=.5; % Kich thuoc tay quay va gia w1=8; % Van toc goc khau dan - rad/s n=51; % So vi tri tay quay df=2*pi/n; f1=0:df:2*pi; f1d=f1*180/pi; % Goc tay quay – do l3=l0*exp(i*pi/2)+l1*exp(j*f1); f3=angle(l3); l3a=abs(l3); l3ac=-l1*w1*sin(f1-f3); w3=l1*w1*cos(f1-f3)./l3a; alpha3=-(v32.*w3+l1*w1*(w3-w1).*sin(f3- f1))./l3a; subplot(311) plot(f1d,f3*180/pi) title(' CHUYEN DONG CUA CU LIT ') ylabel(' C.vi goc Culit - do') grid subplot(312) plot(f1d,w3) ylabel(' Van toc goc Culit - rad/s') grid subplot(313) plot(f1d,alpha3) ylabel(' Gia toc goc Culit - rad/s^2') xlabel(' Goc tay quay') grid 0 50 100 150 200 250 300 350 400 60 80 100 120 CHUYEN DONG CUA CU LIT C .v i g oc C ul it - do 0 50 100 150 200 250 300 350 400 -10 -5 0 5 V an t oc g oc C ul it - ra d/ s 0 50 100 150 200 250 300 350 400 -50 0 50 G ia t oc g oc C ul it - ra d/ s2 Goc tay quay 3. Giải bài toán cân bằng máy Ví dụ 1 : Giải bài toán cân bằng tĩnh Có trục quay trên có gắn các khối l−ợng lệch tâm đặt trong cùng một mặt phẳng nh− hình vẽ. Trục hiện đang mất cân bằng tĩnh.Ta tiến hành cân bằng trục quay bằng ph−ơng pháp số phức nh− sau: Từ điều kiện cân bằng tĩnh ta có: m4 m1 m1r1 + m2r2 + m3r3 + m4r4 + mcbrcb = 0 r1 r4 r2 mcb rcb r3 Biểu diễn các véc tơ d−ới dạng số phức ta đ−ợc: mcbrcb = - (m1r1e i 1ϕ + m2r2ei 2ϕ + m3r3ei 3ϕ + m4r4ei 4ϕ ) Từ đó tính đ−ợc trị số của tích mcbrcb và góc chỉ ph−ơng ϕ của bán kính véc tơ rcb bằng cách dùng hàm abs và hàm angle. mcbrcb = abs( mcbrcb); ϕ = angle(mcbrcb ) Ch−ơng trình tính nh− sau: % File canbt.m % TINH CAN BANG TINH m1=12; m2=15; m3=27;m4=30; % Cac khoi luong r1=.6; r2=.8; r3=.5; r4=.4; % Tri so cac ban kinh vec to f1=pi/2; f2=pi/12; % Cac goc chi phuong f3=-pi/6;f4=2*pi/3; % cua ban kinh ri mrcb=-m1*r1*exp(i*f1)-m2*r2*exp(i*f2)- m3*r1*exp(i*f3)-... m4*r4*exp(i*f4); mr=abs(mrcb) % Tri so mcg*rcb f=angle(mrcb)*180/pi % Goc chi phuong cua ban % kinh vecto rcb (do) Chạy ch−ơng trình trên kết quả nhận đ−ợc là tích mcbrcb ( kí hiệu là mr) và góc chỉ ph−ơng (kí hiệu là f) của bán kính véc tơ của khối l−ợng cân bằng (tính ra độ): mr = 23.3171 f = -147.2962 Nếu chọn tr−ớc r = rcb= 0,6m thì khối l−ợng cân bằng cần đặt vào là: m = mcb = mr/rcb =38,86 kg. Ví dụ 2: Giải bài toán cân bằng động vật quay Bài toán: Hãy cân bằng động trục quay có gắn các khối l−ợng đặt lệch tâm trên những mặt phẳng khác nhau nh− hình vẽ d−ới đây. Biết hai mặt phẳng T và P là hai mặt phẳng cân bằng. m1 T m1 P m2 m2 r1 r1 r2 r2 m4 r4 A B m4 rP mP l2 mP l1 r3 l4 r3 m3 l3 m3 rT mT mT l Các khối l−ợng, độ dài các bán kính véc tơ t−ơng ứng và các kích th−ớc khác đ−ợc cho nh− sau: m1=20 kg, m2=15 kg, m3=17 kg, m4=26 kg, r1=0,6m, r2= 0,3m, r3=0,4m, r4=0,4m, l1=0,3m, l2=0,5m, l3=0,8m, l4=1,2m, l=1,5m . Các bán kính véc tơ r1, r4 t−ơng ứng có ph−ơng thẳng đứng và nằm ngang, còn véc tơ r2, r3 làm với ph−ơng ngang góc và . 02 30=ϕ 03 60−=ϕ Ta giải nh− sau: Giả sử hai đối trọng khối l−ợng mT và mP đã đ−ợc đặt vào hai mặt phẳng cân bằng T và P, các bán kính véc tơ t−ơng ứng là rT và rP , và hiện giờ trục đã đ−ợc cân bằng tĩnh và động. Từ ph−ơng trình cân bằng mô men đối với điểm A của các lực quán tính và ph−ơng trình tổng véc tơ các lực quán tính bằng 0