Giáo trình Cơ sở MATLAB và SIMULINK - Nguyễn Quang Hoàng

đẳng kỹ thuật. Matlab cĩ thể trợ giúp đắc lực các sinh viên và kỹ sư trong việc giải quyết các vấn đề tính tốn số các bài tốn kỹ thuật. Đặc biệt đối với sinh viên và kỹ sư ngành Cơ khí, điện, điện tử, cơ điện tử Matlab là một cơng cụ khơng thể thiếu. Mục đích c a cuốn sách là cung cấp cơ s Matlab và Simulink cho các sinh viên kỹ thuật từ n m thứ hai sau khi đ cĩ các kiến thức cơ bản về tốn, vật lý, cơ học kỹ thuật cũng như kỹ thuật lập trình. Ngồi ra nếu cĩ được thêm các kiến thức về kỹ thuật điều khiển, xử lý dữ liệu số người đọc cĩ thể m rộng thêm được các ứng dụng c a Matlab và Simulink. Nội dung c a cuốn sách này được phân bố trong chín chương. Các chương từ một đến bảy trình bày việc sử dụng các lệnh c a Matlab cho các bài tốn cơ bản như tính tốn trên v ctơ, ma trận, đ họa, giải phương trình vi phân thường, và một số ph p biến đổi tích phân như ouri r, Laplace. Chương giới thiệu về ph n imulink – một cơng cụ sử dụng các khối hàm để mơ ph ng hệ. Trong mỗi chương, sau ph n giới thiệu cách sử dụng các lệnh c a Matlab đều cĩ những ví dụ cụ thể và ph n bài tập thực hành để người học cĩ thể tự thực hành. Trong chương một số bài tốn kỹ thuật thường gặp trong l nh vực cơ học và kỹ thuật được trình bày. Cuốn sách này được biên soạn trên cơ s bài giảng c a tác giả cho sinh viên ngành Cơ điện tử Trường Đại học Kinh doanh và Cơng nghệ Hà nội. Tuy nhiên, cuốn sách này khơng chỉ là tài liệu học tập cho sinh viên các trường đại học định hướng ứng dụng mà cịn là tài liệu học tập cho sinh viên các trường đại học và cao đẳng kỹ thuật. Cuốn sách cũng là tài liệu bổ ích cho các kỹ sư trong cơng việc chuyên mơn c a họ. Trong quá trình biên soạn chúng tơi đ nhận được sự trợ giúp quý báu c a nhiều đ ng nghiệp. Chúng tơi xin cảm ơn G .TSKH. Nguyễn V n Khang và G .T . Nguyễn hong Điền đ xem giúp bản thảo và cĩ nhiều đề nghị cải tiến quý báu. iii Mặc dù đ cố gắng nhiều, nhưng chắc chắn khơng tránh kh i các sai sĩt, tác giả mong muốn nhận được sự gĩp ý c a các bạn đ ng nghiệp và c a các m sinh viên để cĩ điều kiện sửa chữa, hồn thiện hơn trong các l n tái bản sau. Mọi ý kiến đĩng gĩp xin gửi về địa chỉ c a tác giả: T . Nguyễn uang Hồng, Bộ mơn Cơ học ứng dụng, Trường Đại học Bách khoa Hà Nội. E-mail: hoangnq-dam@mail.hut.edu.vn, hoặc T l. 04.3 6 046 . 10 ăm 2010 Tác giả iv M Lời nĩi đ u iii 1.1 Matlab là gì? 1 1.2 Giao diện người sử dụng 2 1.3 Các ph p tính số học cơ bản 3 1.4 Tốn tử gán 4 1.5 Các định ngh a tốn học cơ bản 8 1.6 ố phức 11 1.7 Xử lý lỗi khi gõ lệnh 12 1.8 Kết thúc phiên làm việc với Matlab 12 1.9 Bài tập thực hành 13 , 2.1 V ctơ và các ph p tính trên v ctơ 15 Nhập v ctơ 15 Các ph p tính cộng và trừ hai v cto cùng cỡ 16 Tạo v ctơ cỡ lớn từ các biến sẵn cĩ 17 Tạo v ctơ cĩ các ph n tử cách đều 18 Các đặc trưng c a v ctơ 19 Tích vơ hướng và tích cĩ hướng hai v ctơ 21 Tham chiếu đến các ph n tử v ctơ 22 Một số hàm định sẵn cho v ctơ 23 2.2 Biểu diễn một đa thức b ng v ctơ 23 Nhập đa thức 23 Các ph p tính trên đa thức 24 h p nhân đa thức 24 Phép chia đa thức 25 h p cộng và trừ đa thức 25 Khơng điểm hay nghiệm c a phương trình đa thức 26 Xây dựng đa thức từ các khơng điểm cho trước 26 Giá trị c a đa thức tại một điểm 26 Đạo hàm đa thức 27 2.3 Ma trận và các ph p tính cơ bản trên ma trận 28 Nhập ma trận 28 v Nhân ma trận với một số, ph p cộng và trừ hai ma trận 28 Chuyển vị ma trận 29 h p nhân ma trận 29 Một số ph p tính cơ bản khác 30 Các ma trận đặc biệt 31 Tham chiếu đến các ph n tử c a ma trận 33 Ma trận khối 35 h p nhân mảng hai ma trận cùng cỡ 36 Tính định thức và giải hệ phương trình đại số tuyến tính 37 Tìm hạng c a ma trận 39 Tìm ma trận nghịch đảo và tựa nghịch đảo 41 Trị riêng và v ctơ riêng c a ma trận vuơng 46 hân tích ma trận vuơng A thành tích các ma trận 47 2.4 Bài tập thực hành 51 3.1 Các kiểu dữ liệu 53 Kiểu v ctơ và ma trận 53 Chuỗi ký tự (ký tự, xâu - string) 54 Kiểu ơ, mảng (C ll-Array) 54 Kiểu cấu trúc 54 3.2 oạn thảo Script file trong Matlab 56 Scripts 57 Các hàm, MATLAB - Function 60 3.3 Các vịng lặp và rẽ nhánh 62 Vịng lặp OR 62 Vịng lặp WHILE 63 Lệnh if, cấu trúc if - else - end 65 Cấu trúc switch-case 66 3.4 Các Mat-File 71 3.5 Bài tập thực hành 71 4 Đồ ọ 4.1 Đ họa 2D 73 Đặt màu và kiểu đường cho đ thị 76 Một số tùy chọn khi vẽ đ thị 2D 77 Vẽ nhiều đ thị trên cùng một hệ trục 79 vi Thêm các chú thích 81 Các lệnh axis 82 Đặt giới hạn miền vẽ với lệnh axis 83 Lệnh ubplots 84 Vẽ các đ thị xếp ch ng và lệnh linspace 87 Đ thị th o tọa độ cực và Đ thị với thang chia Logarith 89 Đ thị c a dữ liệu rời rạc 92 Vẽ biểu đ với lệnh contour – vẽ đường đ ng mức 96 Thêm chú thích trên đ thị 101 Vẽ đ thị các hàm cĩ điểm khơng xác định 101 4.2 Các lệnh vẽ trong khơng gian – 3D 103 4.3 Bài tập thực hành 113 5 ị ĩ đ 5.1 Mịn hĩa b ng đa thức 115 5.2 Mịn hĩa b ng hàm e mũ 120 5.3 Nội suy đa thức 122 5.3 Bài tập thực hành 125 6 P â 6.1 Giải phương trình vi phân bậc nhất với od 23 và od 45 127 6.2 Giải các phương trình vi phân bậc hai 133 6.3 Bài tập thực hành 139 7 á ế đổ í â 7.1 h p biến đối Laplac 141 7.2 h p biến đổi Laplac ngược 143 7.3 p dụng Laplac giải phương trình vi phân 145 7.4 h p biến đổi ouri r 148 7.5 h p biến đổi ouri r ngược 151 7.6 h p biến đổi ouri r một tín hiệu rời rạc 151 7.7 h p biến đổi ouri r nhanh (FFT) 152 7.8 Bài tập thực hành 154 8 G ớ ệu ề S u k 8.1 Khái niệm về simulink 157 vii 8.2 Nguyên lý hoạt động và việc thực hành trong simulink 159 Kh i động simulink 159 Các khối trong simulink 160 8.3 Một số ví dụ đơn giản 164 Mơ hình hĩa một phương trình b ng sơ đ khối 164 Mơ ph ng một quá trình động học 166 Mơ hình hĩa một hệ động lực liên tục đơn giản 168 Mơ tả hệ dao động một bậc tự do 169 Mơ ph ng số tay máy một bậc tự do 171 8.4 Đơn giản sơ đ simulink 174 Đơn giản sơ đ simulink b ng khối cn 174 Đơn giản sơ đ simulink b ng khối subsyst m 175 Kết hợp simulink và script fil (m-file) 177 8.5 Xử lý kết quả mơ ph ng 184 8.6 Bài tập thực hành 185 g 9. G á k u 9.1 Bài tốn hệ thanh 189 Hệ thanh t nh định 189 Hệ thanh siêu t nh 191 9.2 Bài tốn uốn phẳng c a d m 193 9.3 Bài tốn qu đạo chuyển động c a viên đạn 199 9.4 Bài tốn bắn trúng đích 203 9.5 Bài tốn dao động 206 Con lắc đơn 206 Con lắc đơn dây tr o đàn h i 208 Con lắc k p 211 Dao động nh c a con lắc lliptic 213 Hệ dao động ba bậc tự do 214 9.6 hân tích động học cơ cấu 216 hân tích động học cơ cấu bốn khâu bản lề 218 9.7 Bài tốn động học ngược rơbốt 221 Tài liệu tham khảo 225 viii Chương 1. Mơi trưng Matlab 1.1 Matlab là gì? MATLAB đưc vit tt t “MATrix LABoratory”, là mt cơng c tính tốn s và mơ phng s, cơng c này đưc phát trin da trên các thư vin hàm tính tốn s vit bng ngơn ng lp trình FORTRAN. Matlab cĩ giao din thân thin vi ngưi s dng. Mt đnh nghĩa khác: Matlab là mt phn mm đa năng tính tốn s, th hin các s liu bng hình nh trc quan và là mt ngơn ng đa năng cung cp mt mơi trưng linh hot cho vic tính tốn k thut. Matlab bao gm nhiu cơng c đ: - thu thp d liu (Data acquisition) - phân tích và x lý d liu (Data analysis and exploration) - hin th hình nh trc quan và x lý hình nh (Visualization and image processing) - to mu và phát trin thut tốn (Algorithm prototyping and development) - mơ hình hĩa và mơ phng (Modeling and simulation) - lp trình và phát trin ng dng (Programming and application development). Matlab cĩ th chy trên hu ht các h máy tính: máy tính sách tay - Labtop, máy tính cá nhân PC, đn các h thng siêu máy tính (super computer). Matlab đưc điu khin bi các tp lnh, tác đng qua bàn phím trên ca s điu khin. Nĩ cũng cho phép mt kh năng lp trình vi cú pháp thơng dch lnh – cịn gi m-file. Các lnh hay b lnh ca Matlab lên đn con s hàng trăm và ngày càng đưc m rng nh các thư vin tr giúp hay do ngưi s dng to ra. Đ khi đng Matlab trong mơi trưng Windows bn ch cn nháy đúp chut vào biu tưng Matlab cĩ trên màn hình - Desktop. Các lnh ca Matlab rt mnh và hiu qu, nĩ cho phép gii các loi bài tốn khác nhau và đc bit khi x lý các d liu cĩ cu trúc kiu véctơ và ma trn. Trong các phiên bn mi, ngưi ta đã đưa vic tính tốn ký t vào phn mm Matlab. 1 1.2 Giao din ngưi s dng Giao din ca Matlab sau khi kích hot đưc th hin như trên hình 1-1 hoc mt dng tương t tùy theo la chn ca ngưi s dng (bng các la chn khác nhau trong mc Desktop cho ta các dng giao din tương ng). 4 5 3 1 2 6 Hình 1-1. Cửa số giao diện của Matlab Các phn t chính ca giao din này gm: 1. Ca s lnh (Command Window). 2. Ca s ghi li các lnh đã thc hin (Command History). 3. Ca s cho bit thư mc hin thi (Current Directory) và danh sách các bin đang s dng (Workspace). 4. Thanh ch đưng dn ca thư mc hin thi, cho phép chn thư mc hin thi. 5. Thanh Shortcut. 6. Nút Start. 2 1.3 Các phép tính s hc cơ bn Đi vi ngưi mi bt đu s dng Matlab, trưc ht ta cn làm quen vi các phép tính s hc (cng, tr, nhân và chia: +, − , *, / ). Vi hai s a và b trong Matlab ta cĩ th thc hin đưc các phép tính s hc như lit kê trong bng dưi đây. Phép tính Cách vit tốn hc cách vit trong Matlab phép cng c= a + b c= a + b phép tr c= a − b c= a − b phép nhân c= a ⋅ b c= a* b phép chia c= a: b c= a/ b phép chia phi c= a: b c= a/ b phép chia trái c= ba/ c= a\ b lũy tha c= ab c = a^b Ví d >> a=3.5; b=7.5; % gan gia tri so cho hai bien a >> a+b % phep cong hai so ans = 11 >> a-b % phep tru hai so ans = -4 >> a*b % phep nhan hai so ans = 26.2500 >> a/b % phep chia hai so, a:b (= phep chia phai) ans = 0.4667 >> a\b % phep chia trai (hieu la b chia cho a, b/a) ans = 2.1429 >> a^b % phep luy thua, a mu b ans = 1.2037e+004 Chú ý rng, mi dịng ch sau du ‘%’ đu đưc MATLAB b qua và đưc s dng đ din gii, bình lun. Dịng lnh kt thúc vi du chm phy ‘;’ s khơng 3 xut kt qu ra, cịn dịng lnh khơng cĩ du ‘;’ kt thúc (đ trng) s đưa kt qu ra khi dịng lnh đưc thc hin. Mt dịng cĩ th đưc kéo dài bng vic đánh ‘’ vào cui dịng và tip tc câu lnh (phép tính) dịng k tip. Th t ưu tiên các phép tốn Khi tính tốn mt biu thc gm nhiu s hng, nhiu phép tính thì th t ưu tiên các tốn t rt quan trng. Th t ưu tiên Tốn t 1 Ngoc đơn 2 Lũy tha 3 Nhân và chia, t trái qua phi 4 Cng và tr, t trái qua phi Ví d Cn tính giá tr hàm xxx2( 4− 5 3 + 1.5) + 5( x + 10) f( x ) = xx2( 2 + 2.5 x + 5) + 1.5( x + 10) = = ti x x0 nào đĩ, chng hn x 5 . >> x=5; >> tuso=x^2*(x^4-5*x^3+1.5)+5*(x+10); >> mauso=x^2*(x^2+2.5*x+5)+1.5*(x+10); >> f=tuso/mauso f = 0.1037 1.4 Phép gán Trong Matlab du bng “=” đưc s dng cho phép gán. Mc dù trong cách vit thơng thưng ta hiu du bng th hin mt phương trình, nhưng trong Matlab du bng đưc đnh nghĩa là phép gán. Đ phân bit gia hai cách th hin ta xét ví d sau. Nu trong ca s lnh ta vit >> x+18=120 Matlab s báo li vi dịng hin th mu đ: ??? x+18=120 Error: The expression to the left of the equals sign is not a valid target for an assignment. 4 Matlab khơng hiu cách bn vit mt phương trình như trên giy, mà nĩ ch cĩ th thc hin đưc các phép tính và gán giá tr tính đưc cho mt bin nào đĩ. Chng hn ta cĩ th gán giá tr ca phép tính 120-18 cho bin x bng cách vit >> x=120-18 Mt cách vit khác ta cĩ th s dng phép gán cho mt phép tính lp trong chương trình tính. Tc là ta cĩ th vit >> x=x+12 Phép gán này ch hat đng đưc nu như trưc đĩ ta đã cĩ giá tr ca bin x. Ví d, chui các phép tính sau đây đưc thc hin >> x=32^3 x = 32768 >> x=x+124 x = 32892 Đ s dng mt bin trong v phi ca phép gán, chúng ta phi gán mt giá tr cho bin đĩ trưc khi s dng. Khi vit các dịng lnh sau Matlab s báo li >> x= 124 x = 124 >> t=x+y ??? Undefined function or variable 'y'. Lý do ca li trên là do bin y chưa đưc gán giá tr trưc khi thc hin cng vi bin x. Trong khi đĩ các dịng lnh sau đưc thc hin đúng. >> x=124 x = 124 >> y=126 y = 126 >> t=x+y t = 250 Trong nhiu phép gán (cĩ th là phép tính trung gian) nu khơng mun kt qu hin ra dưi phép gán thì ta s dng du chm phy (;) vào cui biu thc. Như th trong ca s lnh ta khơng b lãng phí khơng gian. Ví d >> x=124; >> y=126; >> t=x+y t = 250 Chúng ta cũng cĩ th gom nhiu phép gán trên cùng mt dịng lnh. Chng hn >> x=124; y=126; t=x+y t = 250 5 Lưu ý rng, các du “;” trong dịng lnh trên báo cho Matlab bit là ta khơng mun các giá tr ca x và y hin ra dịng dưi. Nu mun các giá tr này đưc hin th, ta thay chúng bng các du phy “,”. Nu ta khơng s dng du “;” hoc du “,” gia các phép gán Matlab s báo li. >> x=124, y=126, t=x+y >> x=124 y=126 t=x+y x = 124 ??? x=124 y=126 t=x+y y = 126 | t = 250 Error: Unexpected MATLAB expression. Khi thc hin nhiu phép tính, nu mun xĩa bt hot xĩa tt c các bin đã s dng đ gim bt khơng gian s dng ca b nh, ta s dng lnh clear tên-bin hoc clear all Trưc khi làm vic đĩ ta nên xem li danh sách các bin đã s dng bng lnh who. Khi thc hin lnh who, Matlab s lit kê cho bn bit tt c các bin đã đưc s dng. Chng hn, vi các lnh đã thc hin trên ta cĩ >> r=10; a=5^5; x=182; y=235; z=x+y z = 417 >> V=4/3*pi*r^3 V = 4.1888e+003 >> t=x+a t = 3307 >> who Your variables are: V a r t x y z Bng lnh whos, ta s nhn đưc nhiu thơng tin hơn v các bin như: tên bin s dng trong b nh, c, s lưng bytes đã s dng, kiu ca bin. >> whos Name Size Bytes Class V 1x1 8 double array a 1x1 8 double array r 1x1 8 double array t 1x1 8 double array x 1x1 8 double array y 1x1 8 double array z 1x1 8 double array Grand total is 7 elements using 56 bytes Vi lnh clear hoc clear all tt c khơng gian làm vic s b xĩa. Khi khơng cn s dng đn mt hay nhiu bin nào đĩ ta cĩ th xĩa chúng bng lnh clear var_name. Chng hn đ xĩa bin V, a và t ta thc hin 6 >> clear V a t Nu mun s dng li dịng lnh đã nhp vào cho các cơng vic tip theo (gi nguyên hoc đ sa đi thành dịng lnh mi), ta s dng hai mũi tên lên xung trên bàn phím (↑, ↓). Các dịng lnh đã nhp vào s xut hin và bn cĩ th sa đi nu cn. Các lnh này cĩ th đưc sa đi ngay ti dịng lnh hin thi. Ta cũng cĩ th copy và paste các dịng lnh ngay trên Command-Window. Đi vi các dịng lnh dài phi vit xung dịng thì trưc khi xung dịng, bn phi kt thúc dịng th nht bng du chm lng – du ba chm (). Ví d >> FirstClassHolders=109; >> SecondClass=100; >> Coach=121; >> Crew=8; >> TotalPeopleonPlane=FirstClassHolders + SecondClass + Coach... + Crew TotalPeopleonPlane = 338 Du ba chm “” phía sau Coach trên dịng đĩ cho bit dịng lnh này chưa kt thúc. Sau du ba chm này bn đánh ENTER, Matlab s chuyn xung dịng dưi và ch các đu vào tip theo. Cho đn đây chúng ta đã bit cách đưa vào và s dng các bin cho các phép tính gán. Các kt qu tính tốn đưa ra màn hình cĩ bn s sau du chm nu s đĩ cĩ phn l. >> co = cos(0.2), si = sin(0.2), tn = tan(0.2) co = 0.9801 si = 0.1987 tn = 0.2027 Đĩ là đnh dng short trong Matlab. Dng này đã đưc mc đnh khi bn khi đng và s dng Matlab. Nu các kt qu vi bn s sau du chm khơng đt đưc đ chính xác yêu cu ca bn, hãy thc hin lnh >> format long Matlab s hin th kt qu tính vi 14 s sau du chm. Hãy quan sát ví d sau >> format long >> co = cos(0.2), si = sin(0.2), tn = tan(0.2) co = 0.98006657784124 si = 0.19866933079506 tn = 0.20271003550867 >> format short 7 >> co = cos(0.2), si = sin(0.2), tn = tan(0.2) co = 0.9801 si = 0.1987 tn = 0.2027 Hãy so sánh các kt qu trên, chú ý rng vi format short thì s hng th tư sau du chm đã đưc làm trịn. Nu bn tính tốn vi ngành tài chính, liên quan đn đơn v đo là tin thì ch cn hai s thp phân sau du chm là đ, khi đĩ bn s dng format bank. Sau đĩ các kt qu tính tốn s đưc làm trịn vi hai s sau du chm. >> format bank >> luonggio=35.55 luonggio = 35.55 >> luongtuan=luonggio*40 luongtuan = 1422.00 Vi các s ln, Matlab s hin th bng ký hiu e mũ. S 4.1264× 105 s đưc vit dng e mũ là 4.1264e + 005 . Nu mun tt c các s đưc hin th dng e mũ thì bn s dng lnh format short e hoc format long e. C th như trong ví d sau: >> format short e >> x=5.125*3.16 x = 1.6195e+001 >> format long e >> x=5.125*3.16 x = 1.619500000000000e+001 Nu bn gõ vào lnh format rat, Matlab s hin th kt qu bng mt phân s gn nht vi kt qu ca bn, ví d >> format rat >> x=5.125*3.16 x = 3239/200 1.5 Các đnh nghĩa tốn hc cơ bn Đ thun tin cho vic tính tốn, trong Matlab ngưi ta đã đnh nghĩa sn rt nhiu đi lưng tốn hc và các hàm cơ bn. Ví d như s π đã đưc đnh nghĩa sn vi = 4 π 3 tên gi là pi, và khi tính th tích hình cu bán kính R theo cơng thc V3 R , trong Matlab ta thc hin như sau >> R = 2; 8 >> V = 4/3*pi*R^3 V = 33.5103 Mt s quen thuc khác trong tốn hc đưc s dng trong nhiu ng dng là s e vi giá tr xp x bng 2.718. Trong Matlab hàm mũ cơ s e , ex đưc vit là exp(x ). Dưi đây là mt s ví d >> exp(1) ans = 2.7183 >> exp(2) ans = 7.3891 Đ tính căn ca mt s x ta s dng hàm sqrt(x ). Ví d >> x = sqrt(9) x = 3 >> y = sqrt(11) y = 3.3166 Đ tính lơgarít t nhiên ca s x , ta s dng hàm log(x ). Ví d >> log(3.2) ans = 1.1632 >> x = 5; log(x) ans = 1.6094 ði vi lơgarít cơ s 10 ta s dng hàm log10(x ) >> x = 3; log10(x) ans = 0.4771 Dưi đây là các hàm tốn hc cơ bn Các hàm tốn hc sqrt(x) Tính căn ca bin x exp(x) Hàm e mũ, Exponential expm1(x) Tính giá tr exp(x)-1 log(x) Logarithm cơ s t nhiên, cơ s e log1p(x) Tính giá tr log(1+x) log10(x) Logarithm cơ s 10 log2(x) Logarithm cơ s 2 pow2(x) Hàm mũ cơ s 2, pow2(x) = 2^x 9 realpow(x,y) Tính x^y, vi x, y là thc reallog(x) Logarithm cơ s t nhiên ca s thc realsqrt(x) Căn bc hai ca s ln hơn hoc bng 0 nthroot(x, n) Căn bc n ca s thc nextpow2(n) Tr li s p đu tiên sao cho 2^p >= abs(n) Các hàm làm trịn và ly phn dư (Rounding and remainder) fix(x) Làm trịn s x bng mt s nguyên gn 0 floor(x) Làm trịn s x bng mt s nguyên bé hơn gn nht ceil Làm trịn s x bng mt s nguyên ln hơn gn nht round Làm trịn s x bng mt s nguyên gn nht mod Ly phn nguyên ca phép chia rem Ly phn dư ca phép chia sign Hàm du Các hàm lưng giác cơ bn như sin, cos, tan, cot đu đưc đnh nghĩa trong Matlab vi đi s đưc mc đnh cho bng radian. Các hàm lưng giác ngưc đưc mc đnh là tr v giá tr radian. Các hàm này đưc đnh nghĩa trùng vi tên hàm thưng dùng, đưc vit bng ch nh. >> cos(pi/4) ans = 0.7071 Đ s dng các hàm lưng giác ngưc như arcsin, arccos, arctan ta ch vic thêm a vào phía trưc tên ca hàm lưng giác, ví d như asin(x), acos(x), atan(x) >> format rat >> atan(pi/3) ans = 1110/1373 Bng dưi đây lit kê mt s hàm lưng giác và các hàm lưng giác ngưc: Các hàm lưng giác - Trigonometric functions sin Hàm sin sind Hàm sin vi đi s là đ sinh Hàm sin Hyperbolic asin Hàm sin ngưc, tc arcsin 10 asind Hàm sin ngưc cho kt qu là đ asinh Hàm sin Hyperbolic ngưc cos Hàm cos cosd Hàm cos vi đi s là đ cosh Hàm cos Hyperbolic acos Hàm cos ngưc, tc arccos acosd Hàm cos ngưc cho kt qu là đ acosh Hàm cos Hyperbolic ngưc tan Hàm tang, Tangent. tand Hàm tang vi đi s là đ tanh Hàm tang Hyperbolic atan Hàm tang ngưc, arctang atand Hàm tang ngưc cho kt qu là đ atan2 Hàm tang ngưc cho kt qu t –pi đn pi atanh Hàm tang Hyperbolic ngưc cot Hàm cotang cotd Hàm cotang vi đi s là đ coth Hàm cotang Hyperbolic acot Hàm cotang ngưc, arccotang acotd Hàm cotang ngưc cho kt qu là đ acoth Hàm cotang hyperbolic ngưc 1.6 S phc Chúng ta cũng cĩ th đưa s phc vào trong Matlab. Nh li rng đơn v phc đưc đnh nghĩa là căn ca −1 . Trong Matlab s phc ký hiu là i hoc j i = −1 Mt s phc cĩ th đưc vit dng z= x + iy , trong đĩ x là phn thc và y là phn o ca z . Vic nhp mt s phc vào Matlab tht đơn gin, vi i đưc mc đinh là đơn v o. Các phép tính trên s phc cũng đưc thc hin mt các d hiu. Ví d vi hai s phc a = 2 + 3i b = 1 - i ⇒ a + b = 3 + 2i Phép tính này đưc thc hin trong Matlab như sau >> format short >> a = 2 + 3i; 11 >> b = 1 - i; >> c = a + b c = 3.0000 + 2.0000i Mt s hàm liên quan đn s phc đưc lit kê trong bng dưi đây Các hàm vi s phc abs Cho đ ln hay modul véctơ phc angle Cho gĩc pha hay argument ca s phc, radian complex To lp d liu phc t các phn thc vào o conj Cho s phc liên hp imag Cho phn o real Cho phn thc isreal Tr li giá tr đúng, 1, nu đi s là thc cplxpair Sp xp các s phc thành cp liên hp Ví d vi hai s phc z1 và z2 ta cĩ: >> format short >> z1=2.5+6.5i >> z1/z2 z1 = 2.5000 + 6.5000i ans = 1.3659 - 0.7073i >> z2=-0.5+4.5i >> real(z1) z2 = -0.5000 + 4.5000i ans = 2.5000 >> z1+z2 >> imag(z1) ans = 2.0000 +11.0000i ans = 6.5000 >> z1-z2 >> abs(z1) ans = 3.0000 + 2.0000i ans = 6.9642 >> z1*z2 >> angle(z1) ans = -30.5000 + 8.0000i ans = 1.2036 1.7 X lý li khi gõ lnh Như đã thy trên, khi thc hin các dịng lnh bn cĩ th làm khơng đúng và Matlab đã báo li cho bn. Error !. Nu khi bn kt thúc dịng lnh bng gõ phím ENTER và nhn ra rng dịng lnh trên cĩ li, bn khơng nht thit phi gõ li dịng lnh đĩ, mà đơn gin bn ch cn s dng các phím mũi tên lên hoc xung, ↑, ↓ , đ hin th li dịng lnh cn sa. Sau khi sa li, và đánh phím ENTER Matlab s đưa ra kt qu. 12 1.8 Kt thúc phiên làm vic vi Matlab Chúng ta va bt đu vi mt s lnh cơ bn ca Matlab, và bây gi cĩ th bn mun ghi li các cơng vic va thc hin và thốt khi Matlab. Làm th nào đ kt thúc Matlab trên màn hình ca bn? Tht đơn gin bn vào exit trong menu File, như khi s dng các chương trình khác trong Windows. Mt cách khác là bn gõ lnh quit vào du nhc trong ca s lnh và Matlab s đĩng li. Kt thúc phiên làm vic. 1.9 Bài tp thc hành S dng Matlab đ thc hin tính các đi lưng sau: 13 1. x = 12 ; 17 13 2. y =16 + 47 17 3. z = 91.25 4. Đúng hai sai. Khi bin y chưa đưc gán giá tr, Matlab cho phép bn thc hin phép gán x= y2 vào trong b nh đ s dng cho vic tip theo. 5. Cho bit th tích ca hình tr cĩ chiu cao h và bán kính đáy r là V= π rh2 , hãy s dng Matlab đ tính th tích hình tr cĩ chiu cao 12 cm và đưng kính đáy là 4 cm. 6. S dng Matlab tính giá tr sin ca π / 3 và biu din kt qu dưi dng mt s hu t. 7. Hãy to mt m-file đ tính kt qu ca sin(π /4),sin( π /3) và sin(π /2) ; biu din chúng bng các s hu t. 13 Chương 2. Véctơ, ma trn và Matlab Các phép tính s hc thơng thưng đưc thc hin trên các s riêng l đưc gi là đi lưng vơ hưng (scalar). Tuy nhiên trong rt nhiu bài tốn ta cn thc hin các phép tính lp li nhiu ln vi nhng b s liu. Cơng vic này đưc thc hin mt cách d dàng trong Matlab khi các b s liu đĩ đưc lưu tr dưi dng mng. Mng trong Matlab cĩ th là mt, hai hay nhiu chiu. Mng mt chiu cịn gi là véctơ, mng hai chiu là ma trn. Véctơ và ma trn là hai đi tưng này đưc Matlab coi như kiu d liu chính ca mình. Và do đĩ trong Matlab cĩ rt nhiu cơng c hu ích x lý các mng s kiu véctơ và ma trn. Chú ý rng trong Matlab véctơ đưc biu din là mt ma trn ct hoc ma trn hàng. 2.1 Véctơ và các phép tính trên véctơ Nhp véctơ Véctơ là mt mng mt chiu cha các s. Matlab cho phép bn to ra các véctơ ct hoc véctơ hàng. Mt véctơ ct cĩ th đưc to ra trong Matlab bng cách lit kê các phn t trong du ngoc vuơng [ ] và gia các phn t là du chm phy “;”. Các véctơ cĩ s phn t tùy ý. Chng hn đ to ra véctơ ct vi ba phn t, ta vit >> a = [2; 1; 4] a = 2 1 4 Các phép tính cơ bn trên véctơ ct cĩ th đưc thc hin vi tên bin đã s dng đ to ra nĩ. Nu mun nhân mt s vi mt véctơ, hay tích vi mt vơ hưng. Gi s mun to ra mt véctơ cĩ các phn t bng ba ln các phn t ca véctơ va đưc to ra, tc là nhân véctơ đã cĩ vi 3. Tt nhiên ta cĩ th gán s 3 này vào mt bin, c, chng hn. Ta thc hin >> c = 3; >> b = c*a 15 b = 6 3 12 Đ to ra mt véctơ hàng, chúng ta lit kê các phn t trong du ngoc vuơng [], và s dng du cách hoc du phy đ nm gia các phn t. Ví d >> v = [2 0 4] v = 2 0 4 Hoc s dng du phy >> w = [1,1,9] w = 1 1 9 Véctơ ct cũng cĩ th chuyn thành véctơ hàng và ngưc li nh phép chuyn v. Gi s cĩ véctơ ct vi n phn t ký hiu bi: v1    v2 v =   ...  v  n  chuyn v ca véctơ này là T = v [v1 v 2 ... vn ]. Trong Matlab, phép chuyn v đưc th hin bng du nháy đơn (’). Chuyn v véctơ hàng cho ta véctơ ct như trong ví d >> a = [2; 1; 4]; >> y = a' y = 2 1 4 Chuyn v véctơ ct cho ta véctơ hàng >> Q = [2 1 3] Q = 2 1 3 >> R = Q' R = 2 1 3 16 Các phép tính cng và tr hai véctơ cùng c Các phép tính cng và tr hai véctơ to ra cho ta mt véctơ mi. Đ thc hin đưc các phép tính cơng hoc tr, hai véctơ phi cùng dng (cùng là ct hoc cùng là hàng) và cùng cĩ s phn t. Các phép tính đưc thc hin vi tên bin ca chúng. Ví d, thc hin phép cng hai véctơ >> A = [1; 4; 5]; >> B = [2; 3; 3]; >> C = A + B C = 3 7 8 Và đây là phép tr hai véctơ hàng >> W = [3,0,3]; >> X = [2,1,1]; >> Y = W – X Y = 1 –1 2 To mt véctơ t các véctơ khác Matlab cho phép bn ni các véctơ li vi nhau đ to ra mt véctơ cĩ nhiu phn t. Gi u và v là hai véctơ ct vi s phn t tương ng là m và n , mà đã đưc to ra trong Matlab. Chúng ta cĩ th to ra véctơ w vi m+ n phn t, trong đĩ m phn t đu là ca véctơ u và n phn t cui là ca véctơ v . Vic này đưc thc hin bng cách vit w= [u; v]. Ví d >> A = [1; 4; 5]; >> B = [2; 3; 3]; >> D = [A;B] D = 1 4 5 2 3 3 Vic này cũng cĩ th thc hin đi vi các véctơ hàng. Đ to véctơ hàng w vi m+ n phn t t hai véctơ hàng u cĩ m phn t và v cĩ n phn t, ta vit w = [u, v]. Ví d 17 >> u = [12, 11, 9] >> v = [1, 4]; >> w = [u, v] w = 12 11 9 1 4 To véctơ hàng cĩ các phn t cách đu Matlab cung cp cơng c đ to ra véctơ cĩ các phn t cách đu vi bưc h , xut phát t giá tr a1 và kt thúc ti giá tr an . Phn t th k ca véctơ này cĩ giá tr là = + − ⋅ aakk 1 ( 1) h . Tốn t hai chm (:) cho ta thc hin vic này vi cú pháp như sau = a [a1 :: han ] Ví d đ to ra mt véctơ cha danh sách các s chn t 0 đn 10 ta vit >> x = [0:2:10] x = 0 2 4 6 8 10 Nh cĩ k thut này mà vic v đ th hàm s y= fx( ) trên đon x= [, a b ] đưc thc hin mt cách d dàng. Đon [,a b ] s đưc chia đu vi bưc h đ nh thành mt véctơ x , sau đĩ tính giá tr ca hàm s ti các đim chia và ta nhn đưc mt = véctơ y cĩ s phn t bng s phn t ca véctơ x , yi fx( i ). Ví d >> x = [0:0.1:1] x = 0 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 1.00 Trưng hp f( x ) = ex ta cĩ >> y = exp(x) y = Columns 1 through 9 1.0000 1.1052 1.2214 1.3499 1.4918 1.6487 1.8221 2.0138 2.2255 Columns 10 through 11 2.4596 2.7183 hoc vi fx( ) = x 2 >> y = x.^2 y = 0 0.01 0.04 0.09 0.16 0.25 0.36 0.49 0.64 0.81 1.00 Chú ý rng phép tính lũy ta đưc ký hiu bi du mũ (^), và phép tính lũy tha trên đưc vit là .^ ; nu trong ví d trên bn gõ vào >> y = x^2 thì Matlab s báo li ! 18 ??? Error using ==> mpower Matrix must be square. Các phép tính .* ./ .^ s đưc nĩi k hơn phn sau. Trong quá trình to véctơ vi các phn t cách đu, bn cĩ th s dng bưc h âm, tc là véctơ bn to ra s gim đu t giá tr đu v giá tr cui. Ví d đ to ra dãy s gim t 100 v 80 vi bưc là 5, ta vit >> u = [100:–5:80] u = 100 95 90 85 80 Mt cách khác đ to ra véctơ hàng vi các phn t cách đu là s dng lnh linspace. Lnh x = linspace(a,b) cho ta mt véctơ hàng x gm 100 phn t...c gii h phương trình đi s tuyn tính. Các ví d dưi đây s ch ra mt cách đ nhn đưc nghim ca phương trình Ax= b . Xét h phương trình 3x+ 2 y = 5 6x+ 2 y = 2 T đây ta bit đưc ma trn h s A và véctơ v phi b : 3− 2  5  = = A− , b   6 2  2  Ta s nhp chúng vào Matlab: >> A = [3 –2; 6 –2]; >> b = [5; 2]; Trưc ht kim tra xem ma trn h s A cĩ phi là ma trn chính qui: >> det(A) ans = 6 Do đĩ tn ti ma trn nghch đo, và tính đưc nghim theo: >> x=inv(A)*b x = -1.0000 -4.0000 43 Cách mà ta va s dng đ nhn đưc nghim ca h phương trình đi s tuyn tính, ch áp dng đưc khi ma trn h s là vuơng và khơng suy bin. Trong h này ta thy s phương trình bng s n cn tìm. Trưng hp s phương trình ít hơn s n ta gi là h ht (thiu) phương trình (underdetermined). Khi đĩ h phương trình cĩ th cĩ vơ s nghim tha mãn. Bi vì ta ch cĩ th xác đnh đưc mt s bin trong các n đĩ, s cịn li cĩ th nhn các giá tr tùy ý. Ví d như h hai phương trình dưi đây vi ba n x+2 y − z = 3 5y+ z = 0 Sau khi bin đi ta nhn đưc x=3 − 7 y z= −5 y Giá tr ca hai bin x và z ph thuc vào bin y . Vi mi giá tr ca y ta s nhn đưc các giá tr tương ng ca x và z . H cĩ vơ s nghim. Nu vit h trên dng ma trn Ax= b vi 1 2− 1  3    A=0 5 1  , b =  0    0 0 0   0 Ta thy ngay det(A )= 0 . Đ bit đưc nghim ca h này cĩ tn ti khơng ta cn so sánh hng ca ma trn A và ma trn m rng C = [A b], nu chúng bng nhau thì tn ti nghim, và cĩ th s dng phép chia trái. Ví d cn gii h phương trình 3x+ 2 y − z = 7 4y+ z = 2 ta thc hin: >> A= [3 2 -1; 0 4 1]; b = [7;2]; >> C = [A b] C = 3 2 -1 7 0 4 1 2 >> rank(A) ans = 2 >> rank(C) ans = 2 >> x=A\b x = 2.0000 0.5000 0 44 Trong trưng hp này Matlab đã cho z = 0 , thc t z cĩ th nhn giá tr bt kỳ. Vi mt h phương trình cĩ vơ s nghim, thì ta cĩ th tìm đưc mt nghim cĩ đ dài nh nht chng hn. Cơng thc ca nghim cĩ đ dài nh nht như sau x= A† b = AT[ AA T ]− 1 b . Ma trn A†= AAAT[ T ]− 1 đưc gi là ma trn ta nghch đo trái ca ma trn A . Trong Matlab đã sn cĩ mt hàm đ tính ma trn ta nghch đo này, đĩ là pinv(A). Như th nghim cĩ chun nh nht (ti ưu) ca h trên đưc tính là >> A= [3 2 -1; 0 4 1]; b = [7; 2]; >> x=pinv(A)*b x = 1.6667 0.6667 -0.6667 Hãy th li vi >> x=A'*inv(A*A')*b x = 1.6667 0.6667 -0.6667 Khi tính tốn vi các đi lưng, ta cn phi lưu ý đn th t ưu tiên các phép tốn, ví d khi gii h phương trình đi s tuyn tính Kx = 0.5 b⇒ x = 0.5 K−1 b Sau đây s tính trong Matlab vi các kiu vit khác nhau >> K=[2 2 3; 4 5 6; 7 8 9]; b=[2 3 4]'; >> x=0.5*K\b % cho ket qua sai x = 0 -2.0000 2.6667 >> x=0.5*(K\b) % cho ket qua dung x = 0 -0.5000 0.6667 >> x=0.5*inv(K)*b x = 0.0000 -0.5000 0.6667 45 >> x=K\b*0.5 x = 0 -0.5000 0.6667 Tr riêng và véctơ riêng ca ma trn vuơng Cho A là mt ma trn vuơng cp n . S λ đưc gi là tr riêng và véctơ khác khơng x là véctơ riêng ca A nu chúng tho mãn điu kin : Ax= λ x hay (A−λ E ) x = 0 vi E là ma trn đơn v. Vi mi giá tr λi , ta s cĩ vơ s các véctơ xi tho mãn. Các véctơ riêng cùng tương ng vi mt λi rõ ràng là ph thuc tuyn tính và ch khác nhau mt hng s α. Do đĩ ta cĩ th chn mt véctơ duy nht làm cơ s (ví d đưc chun hố theo chun Euclide đ cĩ chun bng 1...). Tp hp n véctơ riêng ng vi n tr riêng khác nhau to thành mt h véctơ đc lp tuyn tính. Ma trn gm các ct là các véctơ riêng ca ma trn A gi là ma trn dng riêng ca A (Modal matrix). Nu cĩ hai ma trn vuơng A và B cùng cp n . S λ và véctơ khác khơng x tho mãn điu kin : Ax= λ Bx hay (A− λ Bx ) = 0 thì s λ gi là tr riêng suy rộng ca hai ma trn A và B, véctơ x là véctơ riêng tương ng. Vic tìm tr riêng và véctơ riêng trong Matlab đưc thc hin bng lnh eig(): Lnh eig tính tốn tr riêng và véctơ riêng ca ma trn vuơng cho mt véctơ l cha các tr riêng ca l = eig(A); ma trn vuơng A. cho ta ma trn chéo D cha các tr riêng [V, D] = eig(A); ca A, ma trn V cĩ các ct là các véctơ riêng ca A tha mãn AV = VD. cho ta véctơ cha các tr riêng suy rng d = eig(A,B); ca các ma trn vuơng A và B. cho ta ma trn chéo D cha các tr riêng suy rng ca hai ma trn A và B, ma trn [V, D] =eig(A,B); V cĩ các ct là các véctơ riêng tha mãn AV = BVD. tương t như eig(A,B) đi vi ma trn đi xng A và ma trn đi xng xác đnh eig(A,B,'chol') dương B. Phương pháp tính đây da trên khai trin Cholesky ma trn B. Các ví d sau minh ha vic tìm tr riêng và véctơ riêng trong Matlab 46 >> C=[2 1 0; >> [V, D]=eig(C) 1 4 -1; V = 0 -1 7]; 0.9119 -0.4064 -0.0571 -0.4037 -0.8630 -0.3037 -0.0742 -0.3000 0.9510 >> eig(C) D = ans = 1.5573 0 0 1.5573 0 4.1233 0 4.1233 0 0 7.3194 7.3194 Phân tích ma trn vuơng A thành tích các ma trn Phân tích LU Vi mi ma trn vuơng A khơng suy bin, det(A )≠ 0 , luơn phân tích đưc thành tích ca hai ma trn tam giác A= LU , vi L là ma trn tam giác dưi cĩ các phn t trên đưng chéo chính bng 1, và U ma trn tam giác trên. Trong trưng hp cĩ s dng s hốn v các hàng cho nhau trong quá trình phân tích thì, các phép hốn v s đưc th hin trong ma trn P , tha mãn PA= LU Ma trn hốn v cĩ tính cht sau PPEPPT= ⇒ T = −1 Ví d >> A = [-1 2 0; 4 1 8; 2 7 1] A = -1 2 0 4 1 8 2 7 1 >> [P, L, U]=lu(A) P = 1.0000 0 0 0.5000 1.0000 0 -0.2500 0.3462 1.0000 L = 4.0000 1.0000 8.0000 0 6.5000 -3.0000 0 0 3.0385 47 U = 0 1 0 0 0 1 1 0 0 Ý nghĩa ca phân tích LU cĩ th thy trong vic gii h Ax= b . Sau phi phân tích LU ta nhn đưc Ax= LUx = b ðt y= Ux , suy ra Ly= b . Như th thay vì gii h Ax= b , ta gii hai h Ly= b và Ux= y Do L và U là các ma trn tam giác nên vic gii hai h này khá đơn gin. Vi L là ma trn tam giác dưi nên d dàng nhn đưc nghim y . Sau đĩ gii tìm x . Ví d cn gii h phương trình 3x+ 2 y − 9 z = − 65 −9x + 5 y + 2 z = 16 6x+ 7 y + 3 z = 5 Nhp h vào Matlab, phân tích LU và gii tìm nghim: >> A = [3 2 -9; -9 -5 2; 6 7 3]; b = [-65; 16; 5]; >> [L, U] = lu(A) L = -0.3333 0.0909 1.0000 1.0000 0 0 -0.6667 1.0000 0 U = -9.0000 -5.0000 2.0000 0 3.6667 4.3333 0 0 -8.7273 >> x = U\(L\b) x = 2.0000 -4.0000 7.0000 Phân tích Cholesky Ma trn vuơng A đi xng xác đnh dương luơn phân tích đưc thành tích ALL= T , vi L là ma trn tam giác dưi. Ma trn vuơng A xác đnh dương nu tha mãn xAxT >0, ∀ x ≠ 0 . 48 Ví d >> M = [2 3 4;3 5 6;4 6 9]; >> L=chol(M) L = 1.4142 2.1213 2.8284 0 0.7071 0.0000 0 0 1.0000 >> L'*L-M ans = 1.0e-015 * 0.4441 0 0 0 0 0 0 0 0 Phân tích QR Phân tích QR là phép phân tích mt ma trn vuơng hoc ch nht thành tích ca ma trn trc giao và mt ma trn tam giác trên: A= QR vi QQT = E , R là ma trn tam giác trên. Ví d >> A = [2 3 4 5; 3 5 6 8; 4 6 9 10]; >> [Q, R] = qr(A) Q = -0.3714 0.2491 -0.8944 -0.5571 -0.8305 -0.0000 -0.7428 0.4983 0.4472 R = -5.3852 -8.3563 -11.5131 -13.7415 0 -0.4152 0.4983 -0.4152 0 0 0.4472 -0.0000 >> Q*R - A ans = 1.0e-014 * 0 -0.1332 -0.1776 0 0.0444 -0.0888 -0.0888 0 0 -0.0888 0 0.1776 49 >> Q*Q' ans = 1.0000 -0.0000 -0.0000 -0.0000 1.0000 -0.0000 -0.0000 -0.0000 1.0000 >> inv(Q)-Q' ans = 1.0e-015 * 0 0 0 -0.0833 0.1110 -0.0555 0.1110 -0.1066 -0.2220 Phân tích SVD (Singular Value Decomposition) Phân tích SVD là biu din ma trn A c m× n dng tích các ma trn A= USVT vi U là ma trn trc giao c m× m , V là ma trn trc giao c n× n , và S là ma trn thc c m× n cĩ dng đưng chéo cha các giá tr kỳ d (singular values) ca A theo trt t t ln đn bé. Đĩ là các giá tr khai căn ca các tr riêng ca ma trn AT A . Ví d >> A = [1 2;2 3;3 5]; >> [U,S,V] = svd(A) U = -0.3092 0.7557 -0.5774 -0.4998 -0.6456 -0.5774 -0.8090 0.1100 0.5774 S = 7.2071 0 0 0.2403 0 0 V = -0.5184 -0.8552 -0.8552 0.5184 >> err = U*S*V'-A err = 1.0e-015 * 0 0.4441 0.8882 0.4441 0.8882 0 50 2.4 Bài tp thc hành 1. Tìm đ dài ca véctơ a = [−1 7 3 2]. 2. Tìm đ dài ca véctơ phc b = [−1 + i 7i 3 −2−2i]. 3. Nhp các s 1, 2, 3 đ to ra mt véctơ hàng và mt véctơ ct. 4. Cho hai véctơ ct a = [1; 2; 3] và b = [4; 5; 6]. Hãy tìm tích phn t hai véctơ. 5. Lnh nào cĩ th to ra mt ma trn c 5 × 5 cĩ các phn t trên đưng chéo chính bng 1 cịn li là 0? 6. Cho hai ma trn 8 7 11  2 1 2   AB=6 5 1 , =  1 6 4   0 2 8  2 2 2 Hãy tính tích phn t (tích mng) và tích hai ma trn A và B. 7. Cho ma trn 8 7 11  3 2 8      A = 6 5 1  , t A hãy to ra ma trn B = 3 2 8  .     3 2 8  6 5 1  8. Tìm nghim ca h phương trình : x+2 y + 3 z = 12 −4x + y + 2 z = 13 9y− 8 z = − 1 Đnh thc ca ma trn h s bng bao nhiêu? 9. Hãy xem xét h sau cĩ nghim duy nht khơng? ti sao x−2 y + 3 z = 1 x+4 y + 3 z = 2 2x+ 8 y + z = 3 10. S dng phân tích LU tìm nghim ca h sau: x+7 y − 9 z = 12 2x− y + 4 z = 16 x+ y −7 z = 16 51 Chương 3. Lp trình trong Matlab Trong các chương trưc đã trình bày vic s dng các hàm, lnh đnh nghĩa sn ca Matlab đ gii các bài tốn đơn gin. Ngồi các lnh, hàm đnh sn cĩ này, Matlab cịn là mơi trưng cho phép lp trình gii các bài tốn phc tp, hoc to ra các hàm tin li cho ngưi s dng. Trong chương này s trình bày kh năng lp trình trong Matlab. 3.1 Các kiu d liu Trong Matlab cĩ nhiu kiu d liu, sau đây ch trình bày mt s kiu liên quan đn dng th hin và lưu tr. Đ cĩ th nhn đưc s tr giúp trc tip t Matlab, bn cn tn dng câu lnh help. ==> help class, help strfun, help struct Kiu véctơ và ma trn Biu din s các ma trn hai hay nhiu chiu cũng như các véctor trong chương 2 là mt kiu d liu đc bit. Mi phn t ca véctơ hay ma trn theo chun cn ơ nh 8 Byte (class double) hoc 4 Byte (class single). Các đi lưng phc cn mt ơ nh gp đơi, phn thc và o đưc ghi riêng r. Ví d đi vi mng 3 chiu (3D), vi cu trúc K(hàng, ct, chiu sâu) >> K(1,1,1)=2; K(2,2,2)=4; K(3,3,3)=5; >> K(:,:,3) % lop 3 ans = 0 0 0 0 0 0 0 0 5 K(:,:,2) % lop 2 ans = 0 0 0 53 0 4 0 0 0 0 K(:,:,1) % lop 1 ans = 2 0 0 0 0 0 0 0 0 Chui ký t (ký t, xâu - string) Chui ký t đưc đt trong cp du nháy đơn trên >> ’Day la chuoi ky tu’ Mi phn t ca chui ký t (Ma trn ký t) cn ơ nh 2 Byte. Kiu ơ, mng (Cell-Array) Các d liu ca các kiu khác nhau, ví d chui ký t, các ma trn cĩ c khác nhau, cĩ th đưc s dng như các ơ ca mt mng. Các phn t ơ đưc gi đn thơng qua ch s ca nĩ. Các phn t ca mt ơ đưc đt trong du ngoc nhn {}. Ví d >> A(1,1)={[1 2 3; 4 2 6; 1 7 9]}; >> A(1,2)={'Ma tran thu'}; >> A(2,1)={3+8i}; >> A(2,2)={0:pi/20:2*pi}; %hoac >> A{1,1}=[1 2 3; 4 2 6; 1 7 9]; >> % goi den cac phan tu >> A(1,1) ans = [3x3 double] >> A{1,1} ans = 1 2 3 4 2 6 1 7 9 >> A{1,2} % o (1,2) ans = Ma tran thu >> A(1,2) ans = 'Ma tran thu' Kiu cu trúc Cu trúc đưc s dng đ làm vic vi các kiu d liu khác nhau. Tên ca mt cu trúc gm hai phn, tên cu trúc trưc du chm ( . ) và tên trưng trong cu 54 trúc sau du chm, (struct_name.field_name). Các phn t ca cu trúc đưc gi đn qua tên và ch s. Cú pháp: % tao mot cau truc structur=struct(’name_1’, value1, ’name_2’, value2,.. ) % tro den phan tu cua cau truc structur.name hoc structur.name_1 = value_1; structur.name_2 = value_2; Ví d >> A=[1 2 3; 4 2 6; 1 7 9]; >> my_structur=struct('data', A, 'dimension',[3 3]) my_structur = data: [3x3 double] dimension: [3 3] >> my_structur.dimension ans = 3 3 >> my_structur.data ans = 1 2 3 4 2 6 1 7 9 >> my_structur(2).data=inv(A) % mo rong truong data my_structur = 1x2 struct array with fields: data dimension >> my_structur.data % in ra ma tran A va inv(A) ans = 1 2 3 4 2 6 1 7 9 ans = 4.0000 -0.5000 -1.0000 5.0000 -1.0000 -1.0000 -4.3333 0.8333 1.0000 >> % goi den cac phan tu >> my_structur(1).data(1,2) % phan tu A(1,2) ans = 55 2 >> my_structur(2).data(1,2) % phan tu (1,2) cua inv(A) ans = -0.5000 >> my_structur(2).data(3,2) % phan tu (3,2) cua inv(A) ans = 0.8333 3.2 Son tho Script file trong Matlab Matlab là mt ngơn ng tương tác, nghĩa là các lnh đưc đánh ti du nhc ‘>>’ đưc thc hin ngay sau khi kt thúc dịng lnh bng phím ‘Enter’. Tuy nhiên tht đơn điu khi đánh (gõ) mt chui lnh dài mi khi s dng Matlab thc hin cơng vic. Đ tránh vic này trong Matlab cĩ hai bin pháp m rng cơng năng: scripts và functions (dùng tp văn bn và các hàm). C hai bin pháp này s dng loi tp cĩ tên là m-file (file cĩ phn m rng là .m) đưc to ra nh trình son tho. Ưu đim ca m-file là cĩ th ghi li các dịng lnh và cĩ th sa đi mt cách d dàng mà khơng phi đánh li tồn b danh sách các lnh. Script và Matlab-function function [out_1, . . .] = name(in_1, . . .) % Matlab–function / Hàm trong matlab @ Tham chiu đn hàm function handle nargin S tham s đu vào narout S tham s đu ra global Đnh nghĩa, khai báo bin tng th / tồn cc persistent Đnh nghĩa, khai báo bin persistent %, %{ . . . %} Các dịng chú thích, din gii Xung dịng khi dịng lnh quá dài eval(str) Xác đnh giá tr mt xâu feval(F, in_1, ,in_n) Xác đnh giá tr mt hàm inline(’function’,’t’,) Hàm ngay trong mt script clear function_name Xĩa các hàm == > help function, help function_handle, help funfun 56 Các th tc vào / ra disp(string) Hin th mt chui ký t, mt xâu Hin th d liu khơng đưc đnh disp(variable) dng Chuyn đi mt bin thc thành xâu num2str(variable) ký t Chuyn đi mt bin nguyên thành int2str(variable) xâu ký t variable = input(string) Đưa vào mt bin string = input(string,’s’) Đưa vào mt xâu ký t Dng cho đn khi gõ vào phím bt pause kỳ Dng li mt khong thi gian n pause(n) giây In ra mt d liu cĩ đnh dng fprintf(format, variable) sprintf(format, variable) In ra mt ký t, mt xâu → tic operations toc Cho ta thi gian tính tốn tic toc == > help function, help function_handle, help funfun Scripts biu tưng “t giy trng” Hình 3-1. 57 Matlab Script files là mt chui các lnh đưc gõ vào trong ca s son tho và đưc ghi vào tp cĩ phn m rng là m (filename.m, hay đưc gi là m-file). Đ to m-file, trong menu chính bn chn New ==> M-file hoc bn nháy chut vào biu tưng “t giy trng” trên thanh cơng c (hình 3-1). Vic chy m-file tương đương vi vic đánh tồn b các dịng lnh trên ca s lnh ti du nhc ‘>>’ ca Matlab. Các bin s dng trong m-file đưc đt vào trong khơng gian làm vic ca Matlab. Khơng gian làm vic này là trng khi khi đng Matlab, và nĩ s cha tt c các bin đưc đnh nghĩa trong phiên làm vic. Mun xĩa tt c các bin ta s dng lnh. >> clear all Mt script-file cĩ th s dng các bin trong mơi trưng làm vic và đưa ra mơi trưng làm vic (workspace) các d liu to ra. Các d liu này tn ti trong quá trình chy chương trình và nĩ cĩ th đưc s dng cho các cơng vic x lý tip theo và cĩ th hin th chúng bng đ th chng hn. Script-file khơng cha các phn khai báo và khơng b gii hn bi begin / end. Các chú thích và din gii trong file đưc vit tng dịng mt sau du %. Đi vi MATLAB 7 ta cĩ th son mt khi gm nhiu dịng bt đu bng %{ và kt thúc bi %}. Đ ngt dịng trong mt biu thc hay mt dịng lnh ta s dng ba chm liên tip cui dịng, . Matlab s hiu dịng dưi là ni tip ca dịng trên. Hình 3-2. Ca s son tho script – file, MATLAB 7 58 Vic gi hay chy mt chương trình trong m-file đã đưc son tho trưc cĩ th đưc thc hin bng cách: trong ca s lnh sau du nhc >> ta gõ tên m-file đĩ, cn nh là khơng cĩ phn m rng (.m). Hoc ta cĩ th chy chương trình ngay trên ca s son tho bng cách n F5, hoc dùng chut chn Debug -- > Run trên thanh Menu. Dưi đây là mt Script file thc hin vic tính và đưa th tích ca mt hình cu, mà tham s bán kính đưc ngưi s dng nhp vào t bàn phím function volume r = input('Hay cho ban kinh cua hinh cau') vol = (4/3)*pi*r^3; disp('The tich hinh cau nay la:') disp(vol) Dưi đây là mt s dịng lnh minh ha cho eval, feval, inline, hàm n tên và fprintf, spritf, disp. • eval >> x='1/y*sin(y)'; % Bieu dien ham nhu mot xau ky tu >> y=0.725; % gan cho bien y mot gia tri >> a=eval(x) % a = 0.9147 Nu cĩ mt hàm foo.m đã đưc son trong mt m-file thì sau dịng lnh >> fooHandle=@foo ba lnh sau đây s cho kt qu như nhau. >> feval(fooHandle,8.75) >> fooHandle(8.75) >> foo(8.75) • inline >> g=inline('5*a+6.5*sin(b)') % hoac >> g=inline('5*a+6.5*sin(b)','a','b') % a, b la hai bien g = Inline function: g(a,b) = 5*a+6.5*sin(b) >> c=0:2 c = 0 1 2 >> b=pi/4 b = 59 0.78539816339745 >> z=g(c,b) z = 4.59619407771256 9.59619407771256 14.59619407771256 Hàm n tên, cú pháp fhandle = @(argumlist) . . . . . >> om=5; % gan gia tri 5 cho bien om >> f = @(y) cos(om*y); % viet bieu thuc ham >> t=linspace(0,2*pi/om); % chuoi thoi gian >> f(t) Đnh dng hin th >> x='1/y*sin(y)'; % Bieu dien ham nhu mot xau ky tu >> y=0.725; % gan cho bien y mot gia tri >> a=eval(x) % a = 0.91466957608738 >> fprintf('%3.3f',a) 0.915 >> fprintf('% s %3.3f','Gan cho bien a gia tri', a) Gan cho bien a gia tri 0.915 >> sprintf('%3.3f',a) ans = 0.915 >> sprintf('%3.3e',a) ans = 9.147e-001 >> disp(['a = ', num2str(a, '%4.3f')]) a = 0.915 >> disp(['a = ', sprintf('%4.3f',a)]) a = 0.915 Các hàm, MATLAB - Function Ngồi các hàm sn cĩ (built-in functions), Matlab cho phép to ra các hàm ca riêng bn, đĩ là các hàm dng m-file. Các hàm m-file này cĩ phn m rng là .m tương t như các Script file. Tuy nhiên, điu khác vi các script file là các hàm m- file đưc s dng vi các tham s vào và tham s ra. Vic s dng các hàm m-file là rt cn thit, chúng cho phép chia nh mt bài tốn ln thành nhiu mơđun, mi 60 mơđun s đưc gii bng mt hàm. Các hàm m-file dng này đưc gi là Matlab- function vi cu trúc tng quát như sau: function [output 1, output 2] = function_name(input1,input2) % % Some comments that explain what the function does go here. % MATLAB command 1; MATLAB command 2; MATLAB command 3; Tên ca m-file lưu tr hàm này phi ly đúng tên hàm đã đnh nghĩa trong chương trình function_name và nĩ đưc gi t dịng lnh ca MATLAB hoc t mt m-file khác như lnh sau: >> [output1, output2] = function_name(input1, input2) Khi làm vic vi các hàm m-file ta cn phân bit hai loi bin đưc s dng: Bin tồn cc (global variable) và bin cc b hay bin đa phương (local variable). Bin tồn cc đưc khai báo theo cú pháp global var_1 var_2 Các bin tồn cc đưc xĩa khi mơi trưng tính tốn bng lnh clear [tên bin]. Bin cc b ch đưc s dng trong mt hàm, bin này khơng đưc lit kê trong ca s workspace, và do đĩ khơng th tác đng đn đưc trong mơi trưng làm vic. Ngồi ra, các bin persistent trong mt hàm cĩ th đưc thng nht vi cách khai báo persistent var_1 var_2 Trái vi các đi lưng đã đưc khai báo bng global, các đi lưng vi khai báo persistent ch đưc bit đn trong hàm mà nĩ đưc khai báo. Do đĩ các hàm khác khơng th tác đng đưc đn bin này. Các bin persistent ch đưc xĩa khi hàm cha nĩ thốt ra khi b nh (clear function_name) hoc khi hàm đưc thay đi sau đĩ đưc ghi li. Sau đây s trình bày hai ví d đ làm sáng t vic thao tác đi vi hàm. Ví d đu tiên là vic son mt hàm đ v các đưng trịn cĩ các bán kính khác nhau vi tên hàm (file_name=fcircle.m, function_name=fcircle). Thơng s đu vào là bán kính r , đu ra là ta đ x, y ca các đim trên đưng trịn. Trong hàm này cĩ s dng lnh plot đ v đưng trịn. Hàm này đưc vit như sau: function [x,y]=fcircle(r) % chia nho goc 2pi thanh 100 goc nho bang nhau goc=linspace(0,2*pi,100); x=r*cos(goc); % toa do x y=r*sin(goc); % toa do y plot(x,y), grid on; % ve duong tron 61 axis('equal'); % hai truc cung ty le axis([-1.1*r 1.1*r -1.1*r 1.1*r]); title(['Duong tron co ban kinh r = ',num2str(r)]); Trong mơi trưng làm vic (Workspace) hàm này cĩ hai cách gi vi sau >> fcircle(3) % ve duong tron ban kinh r = 3 >> [x,y]=fcircle(4)% ve duong tron ban kinh r = 4, va dua % ra man hinh hai vector chua toa do cac diem Hãy thc hin và xem kt qu ! 3.3 Các vịng lp và r nhánh Khi lp trình tính tốn trên d liu kiu ma trn, vectơ các vịng lp điu khin rt hu ích và đưc s dng rng rãi. Matlab đưa ra các dng vịng lp và r nhánh: vịng lp for, vịng lp while, cu trúc if-else-end và cu trúc switch-case. Các vịng lp và r nhánh for var_=dieukien, Lnh for các lnh end while dieukien, Lnh while các lnh end if dieukien Lnh if (nu .. thì ..) ... end switch .. Lnh switch (chuyn) case .. end break Lnh ngt vịng lp trong for và while Vịng lp FOR Vịng lp for cho phép thc hin mt nhĩm lnh lp li vi s ln xác đnh vi cú pháp như sau 62 for index = start: increment : finish statement_1, ..., statement_n end Ví d cn xây dng mt hàm s dng vịng lp for đ tính tng các phn t ca mt véctơ v. Trên thanh cơng c, ta la chn New –> m-file đ m mt ca s son tho, trong ca s này ta son ni dung sau function tong = tongvector(v) % ham tinh tong cac phan tu cua mot vector n=length(v); % tra lai so phan tu, chieu dai cua vector s = 0; % khoi gan cho tong for i = 1:n % hoac for i = 1 : 1 : n s = s + v(i); end; tong = s; % save with file name tongvector.m Vi hàm đnh nghĩa như vy, ta cĩ th s dng ngay trên du nhc trong ca s lnh, ví d >> x=[1:2:50]; >> s=tongvector(x) s = 625 Đ kim tra vi lnh sn cĩ ca Matlab ta gi lnh sum >> sum(x) ans = 625 Mt s ví d khác v s dng lnh for: n=36; n=5; for i = 1:n for i = 1:n a(i) = sin(i*pi/n); for j = 1:n end a(i,j) = 1/(i+j-1); end end Vịng lp WHILE While là lnh lp vi s ln lp khơng xác đnh. Cú pháp WHILE expression statements END 63 Ví d cn tính tng sau 1 1 1100 1 S =1 + + + ... + = ∑ 2 3 100 i=1 i % Tinh tong S = 1+1/2+1/3+ ... +1/100 % su dung vong lap while % su dung vong lap for function tongS function tong_Sn n=100; n=100; S=0; i=1; S=0; while i<=n for i=1:n S = S+1/i; S=S+1/i; i = i+1; end end tong=S; tong=S; disp('Voi n = 100, tong la :') disp('Voi n = 100, tong la :') disp(tong) disp(tong) % Ket qua la % Ket qua la >> tong_S >> tong_Sn Voi n = 100, tong la : Voi n = 100, tong la : 5.1874 5.1874 Đ vit thành mt hàm cĩ giao din vi ngưi s dng ta vit mt hàm như sau function S = tong_Sn(n) % su dung vong lap for T=1; while T disp('Vao so duong n: ') n = input('n = ') if n>1 T = 0; end end tong=0; for i=1:n tong=tong+1/i; end S=tong; disp(['Voi n = ', sprintf('%4d',n),' Tong S(n) la']) disp(S) % save with file name tong_Sn.m 64 Dưi đây là mt s ví d v vic s dng lnh while: while 1 disp('Vao so n. Neu n <= 0, thi thoat chuong trinh.') n = input(n = ') if n <= 0, break, end r = rank(magic(n)) end disp('Nhu vay do !.') E = 0*A; F = E + eye(size(E)); N = 1; while norm(E+F-E,1) > 0, E = E + F; F = A*F/N; N = N + 1; end Ví d tìm s n ln nht sao cho tng sau nh hơn 7 1 1 1n 1 S =1 + + + ... + =∑ < 7 2 3 ni=1 i function n = timn(S) % tim so n lon nhat thoa man % tongS(n) = 1+1/2+1/3+...+1/n < S % gia tri S > 1 vao tu ban phim T=1; while T disp('Vao so duong S') S = input('S = ') if (S>1 & S<10) T = 0; end end format long tong=0; i=0; while(tong<S) i = i+1; tong = tong+1/i; end n=i-1; tong=tong-1/i; disp('So n tim duoc la: ') disp(n) disp('Gia tri cua tong S(n) la: ') disp(tong) 65 Lnh if, cu trúc if - else - end Lnh điu kin nu - thì. Cú pháp IF expression statements ELSEIF expression statements ELSE statements END Ví d if I == J A(I,J) = 2; elseif abs(I-J) == 1 A(I,J) = -1; else A(I,J) = 0; end k=5; % ket qua chuong trinh for i=1:k for j=1:k >> A if i==j A = A(i,j)=4; 4 0 1 0 0 elseif(abs(j-i)==2) 0 4 0 1 0 A(i,j)=1; 1 0 4 0 1 else 0 1 0 4 0 A(i,j)=0; 0 0 1 0 4 end end end Ma trn A trên cĩ th nhn đưc bng các lnh sau >> k=5; >> d1=4*ones(k,1); d2=ones(k-2,1); >> a=diag(d1)+diag(d2,-2) +diag(d2, 2) Cu trúc switch-case Lnh chuyn trong nhiu trưng hp da trên biu thc và cĩ dng tng quát như sau 66 SWITCH switch_expr CASE case_expr, statement, ..., statement CASE {case_expr1, case_expr2, case_expr3,...} statement, ..., statement ... OTHERWISE, statement, ..., statement END Ví d method = 'Bilinear'; switch lower(method) case {'linear','bilinear'} disp('Method is linear') case 'cubic' disp('Method is cubic') case 'nearest' disp('Method is nearest') otherwise disp('Unknown method.') end % ket qua la Method is linear Ví d diem = 1; switch diem case 1 disp('Diem gioi') case 2 disp('Diem tot') case 3 disp('Diem kha') case 4 disp('Diem trung binh') otherwise disp('Diem kem.') end 67 Đ thốt ra khi vịng lp for hoc while, ta s dng lnh break. Khi chương trình chy đn dịng lnh này thì nĩ t đng nhy ra khi vịng lp cĩ cha break. Ví d for p = 7:8 for q = 3:5 for r = 1:2 fprintf('\n %3.0f, %3.0f, %3.0fn' , p, q, r) end if q == 4, break; end end end fprintf('\n Xuong dong \n') % Ket qua la 7, 3, 1n 7, 3, 2n 7, 4, 1n 7, 4, 2n 8, 3, 1n 8, 3, 2n 8, 4, 1n 8, 4, 2n Xuong dong Ví d sau đây ch ra s tương tác gia các hàm khi lp trình. Xét phương trình vi phân = = yɺ f(, t y ) vi điu kin đu y( t0 ) y 0 . = Ta cn tìm nghim y( t ) trong khong thi gian t[, t0 tf ]. Đ gii bài tốn bng phương pháp Euler, ta s chia khong thi gian trên thành N = − = đon, vi đ dài hay cịn gi là bưc thi gian là h()/ tf t0 N . Gi zi zt( i ) = + = là tr gn đúng ca hàm y( t ) ti thi đim ttihii 0 , ( 1,2,..., N ) , theo phương pháp Euler dãy các giá tr zi đưc tính theo cơng thc sau = = + ⋅ z0 y(), t 0 zi z i− 1 h f (,) t i − 1 z i − 1 Khi tin hành thc hin trong Matlab ta s son ba m-file: mt mơ t phương trình vi phân, mt mơ t phương pháp gn đúng Euler, và mt chương trình chính. C th đ gii phương trình vi phân =1 − = = yɺ m ( F0 byyt ), ( 0 ) 0 , t0 0, tf =10 . Ta bit gnhim chính xác ca phương trình này là − =F0 − (/b mt ) y() tb [1 e ] 68 Dưi đây là ba m-file đưc son. % file_name = ptvp_vantoc.m function dydt=ptvp_vantoc(t,y) global m b F0 % khai bao cac bien tong the dydt=1/m*(F0-b*y); = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = % file_name = Euler_method.m function [t, z]=Euler_method(ptvp, t0, tf, y0, N) % Cac bien: % ptvp - ten file chua ve phai cua phuong trinh vi phan can giai % t0 - thoi diem dau % tf - thoi diem cuoi (t0 < tf) % y0 - gia tri ham tai thoi diem dau % N - so doan ma khoang thoi gian (tf-t0) duoc chia ra global m b F0 % khai bao cac bien tong h = (tf - t0)/N; % buoc thoi gia t = t0+[0:N]'*h; % vector chua cac diem c z(1,:) = y0; for k = 1:N z(k + 1,:) = z(k,:)+h*ptvp(t(k),z(k,:)); end = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = % file_name = tuy_chon__.m % chuong trinh chinh global m b F0 % khai bao cac bien tong the m=1; b=0.2; F0=4; % Gan gia tri cho cac tham so t0=0; tf=10; N=20; y0=0; % nghiem gan dung z [t, z]=Euler_method(@ptvp_vantoc, t0, tf, y0, N); %nghiem chinh xac y=F0/b*(1-exp(-b/m.*t)); figure(1) subplot(211) plot(t,y,'k',t,z,'-.k','Linewidth',1); % so sanh ket qua tinh voi nghiem chinh xac legend('y','z') ylabel('y, z [m/s]'); subplot(212) plot(t,abs(z-y),'k','Linewi... khơng khí: cl=1, c n = 0 Qui dao chuyen dong khi bo qua luc can khong khi 2500 2000 1500 y y [m] 1000 500 0 0 500 1000 1500 2000 2500 3000 3500 4000 x [m] Hình 915. Quĩ đo chuyn đng ca viên đn khi cĩ cn tuyn tính vi các gĩc bn khác nhau 10, 20, . , 90o T các mơ phng trên, ta cĩ th tng hp các kt qu trong bng 93. Bng 92. Bng 93. Kt qu mơ phng khi khơng cĩ cn Kt qu mơ phng khi cĩ cn t l Thi Thi Gĩc bn Tm xa Tm cao Gĩc bn Tm xa Tm cao gian bay gian bay (đ) (m) (m) (đ) (m) (m) (s) (s) 10 10.60 3146.461 138.319 10 9.95 2435.036 121.402 20 20.90 5905.968 536.595 20 18.60 3705.257 422.214 30 30.55 7950.113 1146.789 30 26.10 4213.168 823.918 40 39.30 9043.155 1895.302 40 32.45 4179.547 1263.938 50 46.85 9044.022 2691.851 50 37.70 3755.920 1691.532 60 52.95 7950.000 3440.366 60 41.80 3046.891 2066.010 70 57.45 5899.847 4050.560 70 44.75 2137.729 2355.927 80 60.20 3138.691 4448.835 80 46.55 1100.429 2538.861 90 61.15 0000.000 4587.154 90 47.15 0000.000 2601.335 c) Trưng hp lc cn khơng khí t l bình phương vn tc: cl=0, c n = 0.1 202 Qui dao chuyen dong khi bo qua luc can khong khi 400 300 y y [m] 200 100 0 0 100 200 300 400 500 600 x [m] Hình 916. Quĩ đo chuyn đng ca viên đn khi cn t l bình phương vn tc vi các gĩc bn khác nhau 10, 20, . , 90o T các mơ phng trên, ta cĩ th tng hp các kt qu trong bng 94 sau Bng 94. Kt qu mơ phng khi cĩ cn t l bc nht và bc 2 Gĩc bn (đ) Thi gian bay (s) Tm xa (m) Tm cao (m) 10 5.40 498.857 39.726 20 8.70 578.956 103.061 30 11.35 592.698 173.318 40 13.65 566.398 244.065 50 15.65 507.712 310.553 60 17.35 419.968 368.575 70 18.60 304.410 414.172 80 19.45 162.797 443.653 90 19.70 000.000 453.958 9.4 Bài tốn bn trúng đích Xét bài tốn: Mt qu đn khi lưng m bay t do trong khơng khí chu lc cn 2 khí đng hc FD = cv , vi v là vn tc và c là h s cn. Phương trình vi phân phân chuyn đng ca viên đn là 203 mxɺɺ= − cvx ɺ, myɺɺ=− mg − cvy ɺ, v= xɺ2 + y ɺ 2 Bit rng vn tc ca viên đn khi thốt khi ming là nịng v0 . Hãy xác đnh gĩc nghiêng θ đ viên đn đt mc tiêu khong cách 8 km. S dng các s liu −4 2 m = 20 kg, c =3.2 × 10 kg/m ,v0 = 500 m/s , và g = 9.80665 m/s . Phương pháp gii bài tốn như sau, trưc ht ta bn th vi gĩc bn θ nào đĩ gii phương trình vi phân chuyn đng ta nhn đưc đim tip đt vi tm xa L= L(θ ) , so sánh đim rơi này mc tiêu ta đưc mt sai lch r()θ= L () θ − Lmt T phương trình sai s này cho ta xác đnh đưc gĩc bn cn thit. Trong matlab ta cn xây dng các mfile sau: function F = dEqs(t,y) % First-order differential equations. m = 20; c = 3.2e-4; g = 9.80665; x_ = y(1); vx = y(2); y_ = y(3); vy = y(4); v = sqrt(vx^2+vy^2); xdot = vx; vxdot = 1/m*(-c*v*vx); ydot = vy; vydot = 1/m*(-m*g-c*v*vy); F = [xdot, vxdot, ydot, vydot]'; function y = inCond(v0, theta) y = [0 v0*cosd(theta) 0 v0*sind(theta)]'; function r = residual(theta) global tSTART tSTOP h k v0 k=k+1 [tSol,ySol] = ode45(@dEqs,[tSTART:h:tSTOP],inCond(v0, theta)); % tim thoi gian chuyen dong for j = 50:length(tSol) if ySol(j,3)>=0 && ySol(j+1,3)<0 t_ground(j) = tSol(j); ki = j+1; end end % Tam xa L(theta) tamxa = ySol(ki, 1); r = tamxa-8000; function root = newtonRaphson2(func,x,tol) if nargin == 2; tol = 1.0e-8; end if size(x,1) == 1; x = x'; end % x must be column vector for i = 1:30 [jac,f0] = jacobian(func,x); if sqrt(dot(f0,f0)/length(x)) < tol root = x; return end 204 dx = jac\(-f0); x = x + dx; if sqrt(dot(dx,dx)/length(x)) < tol*max(abs(x),1.0) root = x; return end end error('Too many iterations') function [jac,f0] = jacobian(func,x) % Returns the Jacobian matrix and f(x). h = 1.0e-4; n = length(x); jac = zeros(n); f0 = feval(func,x); for i =1:n temp = x(i); x(i) = temp + h; f1 = feval(func,x); x(i) = temp; jac(:,i) = (f1 - f0)/h; end function main global tSTART tSTOP h k v0 k=0; %so lan ban thu tSTART = 0; tSTOP = 150; h = 0.01; v0 = 500; theta = 40; % thu ban voi theta = 40 do theta = newtonRaphson2(@residual,theta); [tSol,ySol] = ode45(@dEqs,[tSTART:h:tSTOP],inCond(v0, theta)); % tim thoi gian chuyen dong for j = 50:length(tSol) if ySol(j,3)>=0 && ySol(j+1,3)<0 t_ground(j) = tSol(j); ki = j+1; end end % Tam xa L(theta) tSTOP = tSol(ki) [t,y] = ode45(@dEqs,[tSTART:h:tSTOP],inCond(v0, theta)); plot(y(:,1), y(:,3), 'k-', 'Linewidth', 2), grid on xlabel('x'), ylabel('y') disp('Goc ban can tim (do)') theta Kt qu thu đưc gĩc bn cn thit đ đt mc tiêu là θ = 77.8777o . Thi gian bay ca viên đn là 91.35 giây. Quĩ đo chuyn đng ca viên đn như trên hình 917. 205 15000 10000 y 5000 0 -5000 0 1000 2000 3000 4000 5000 6000 7000 8000 9000 x Hình 917. Quĩ đo chuyn đng ca viên đn 9.5 Bài tốn dao đng Con lc đơn Xét con lc tốn là mt qu cu nh, khi lưng m treo vào dây mm khơng khi lưng, khơng giãn. Phương trình vi phân chuyn đng ca con lc trong trưng hp khơng cn: ml2ϕɺɺ + mgl sin ϕ = 0 ⇒ϕɺɺ =−(g / l )sin ϕ Nu k đn lc cn t l vi vn tc h s k , lc cn l F= kv = klϕɺ , ta cĩ c ϕ ml2ϕɺɺ=− mglsin ϕ − kl 2 ϕ ɺ m ðt y1=ϕ, y 2 = ϕɺ , ta cĩ Hình 918 yɺ1= y 2 1 yɺ =−( mgl sin ykly +=−2 ) (/)sin gl y − (/) kmy 2ml 2 1 2 1 2 Trin khai trong Matlab function ydot = conlac_don(t,y) global m l k g = 9.81; % gia toc trong truong v = l*y(2); ydot = zeros(2,1); ydot(1) = y(2); ydot(2) = -g/l*sin(y(1)) - k/m*y(2); 206 Trong phn sau ta s kho sát chuyn đng ca con lc cĩ khi lưng m = 0.25 kg, chiu dài dây l = 0.5 m, vi các gĩc lch ban đu khác nhau t 10 đn 110 đ. Các kt qu mơ phng đưc đưa ra như trên hình 919, 920 và 921. global m l k m = 0.25; l = 0.5; k = 0.0; % khong can % k = 0.1; % co can t_start = 0; t_end = 6.0; phi = [10 : 20 : 110]; for i = 1:length(phi) y1 = phi(i)*pi/180; y2 = 0; y0 = [y1 y2]'; [t,y] = ode45('conlac_don',[0: 0.02: t_end],y0); plot(t, y(:,1),'k-', 'linewidth',1), grid on xlabel('t [s]'), ylabel('\phi [rad]'), hold on end 2 1 0 [rad] φ -1 -2 0 1 2 3 4 5 6 t [s] Hình 919. ð th theo thi gian ca gĩc lc khi khơng cĩ cn, k = 0 2 1 0 [rad] φ -1 -2 0 1 2 3 4 5 6 t [s] Hình 920. ð th theo thi gian ca gĩc lc khi cĩ cn, k = 0.2 207 0.4 0 0.1 0.2 0.2 0 [rad] y [m] φ 0.3 -0.2 0.4 0.5 -0.4 0 2 4 6 -0.2 0 0.2 t [s] x [m] Hình 921. Chuyn đng ca con lc, vi k = 0.20 Xét trưng hp dao đng nh, ta cĩ th tuyn tính hĩa quanh v trí cân bng và nhn đưc phương trình vi phân: yɺ1= y 2, y ɺ 2 =− (/)(/) g l y 1 − k m y 2 Kt qu mơ phng trong trưng hp khơng cn vi các điu kin đu khác nhau đưc đưa ra như trên hình 922. T trên hình v ta thy rng con lc dao đng vi chu kỳ khơng đi, khơng ph thuc vào điu kin đu. 0.5 0 [rad] φ -0.5 0 1 2 3 4 5 6 t [s] Hình 922. Dao đng nh, khơng cn k = 0 Con lc đơn dây treo đàn hi Xét con lc là mt qu cu nh đưc treo vào lị xo cĩ đ cng c , chiu dài t nhiên là l (hình 923). H hai bc t do vi các ta đ suy rng s,ϕ . ðng năng và th năng ca h: 1 2 2 21 2 Tmls=2 [() ++ϕɺ sɺ ], Π=−+2 csmgls ()cosϕ S dng phương trình Lagrange loi 2 ta nhn đưc phương trình vi phân chuyn đng ca con lc như sau: 208 mls(+++++ )2ϕɺɺ 2( mlss )ɺ ϕ ɺ mgls ( )sin ϕ = 0 msɺɺ −+− m( l s )ϕɺ2 mg cos ϕ += cs 0 Trin khai trong Matlab c l+ s function ydot = conlac_danhoi(t,y) ϕ % ptvp cd cua con lac dan hoi (m, c, L) % vector y = [q1, q2, v1, v2]' global m c L g = 9.81; m q1 = y(1); q2 = y(2); qd1 = y(3); qd2 = y(4); Hình 923 qd = [qd1; qd2]; M_q = [m, 0; 0, m*(L+q1)^2]; C_q = [ 0, -m*(L+q1)*qd2; m*(L+q1)*qd2, -m*(L+q1)*qd1]; g_q = [c*q1-m*g*cos(q2); m*g*(L+q1)*sin(q2)]; ydot = zeros(4,1); qdot = qd; qd_dot = inv(M_q)*(-C_q*qd - g_q); ydot = [qdot; qd_dot]; % save conlac_danhoi.m = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = % main program global m L c m = 0.5; L = 0.5; c = 50; % gan cac thong so he t_start = 0; t_end = 3.0; q10 = 0.05; q20 = 30.0; y1 = q10; y2 = q20*pi/180; y3 = 0; y4 = 0; y0 = [y1 y2 y3 y4]'; [t,y] = ode45('conlac_danhoi',[0: 0.01: t_end],y0); figure(1) subplot(2,1,1), plot(t, y(:,1),'k-','linewidth',1), grid on ylabel('q_1 [m]') subplot(2,1,2), plot(t, y(:,2),'k-','linewidth',1), grid on ylabel('q_2 [rad]'), xlabel('t [s]'), figure(2) x0 = 0; y0 = 0; hold on, plot(x0,y0,'o-','linewidth',1) r = 0.02; for i = 1:5:length(t)/2 x1(i) =(L+y(i,1))*sin(y(i,3)); y1(i)=-(L+y(i,1))*cos(y(i,3)); x_day = [x0, x1(i)]; y_day = [y0, y1(i)]; 209 plot(x_day, y_day, 'k-','linewidth',1), grid on, axis equal tam1 = [x1(i), y1(i)]; vehinhtron(tam1,r); end xlabel('x [m]'), ylabel('y [m]') axis([-0.8*L 0.8*L -1.5*L 0.05*L]); figure(3) for i = 1:1:length(t) xC(i) = (L+y(i,1))*sin(y(i,3)); yC(i) = -(L+y(i,1))*cos(y(i,3)); end plot(xC,yC,'k-','linewidth',1), grid on = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = Dưi đây là kt qu mơ phng vi các s liu khác nhau ca lị xo, và các điu kin đu khác nhau. 0.2 0.1 [m] 0.1 [m] 0.05 1 1 q q 0 0 0 1 2 3 0 1 2 3 1 1 0 0 [rad] [rad] 2 2 q q -1 -1 0 1 2 3 0 1 2 3 t [s] t [s] -0.45 -0.45 -0.5 -0.5 -0.55 -0.55 y [m] y y [m] y -0.6 -0.6 -0.65 -0.7 -0.65 -0.4 -0.2 0 0.2 0.4 -0.4 -0.2 0 0.2 0.4 x [m] x [m] 210 0 0 -0.1 -0.1 -0.2 -0.2 -0.3 -0.3 y [m]y y y [m] -0.4 -0.4 -0.5 -0.5 -0.6 -0.2 0 0.2 -0.2 0 0.2 0.4 x [m] x [m] m = 0.5; L = 0.5; c = 50; m = 0.5; L = 0.5; c = 150; q10 = 0.05; (mét) q10 = 0.05; (mét) q20 = 30.0; (đ) q20 = 30.0; (đ) Hình 924. Kt qu mơ phng con lc đàn hi Con lc kép Phương trình vi phân chuyn đng ca con lc kép nhn đưc nh phương trình Lagrănge 2 đưc vit gn li dng ma trn như sau: l1 Mqq( )ɺɺ+ Cqqq ( , ɺ ) ɺ += gq ( ) 0 ϕ 1 m vi 1 2 2 l2 m1 l 1+ m 2 l 2 m 2 l 1 l 2cos( q 1 − q 2 )  M( q ) =   mllcos( qq− ) ml 2 ϕ2 m2 2 1 2 1 2 2 2   0m2 l 1 l 2 sin( q 1− q 2 ) qɺ 2  Cq(, qɺ ) =   Hình 925. Con lc kép −m l lsin( q − q ) qɺ 0  2 1 2 1 2 1  T g( q )= (m l + m l ) g sin q m gl sin q  . 1 1 2 2 2 2 2 2  Trin khai trong Matlab function ydot=conlactoan_kep(t,y) % vector y = [q1, q2, v1, v2]' global m1 m2 l1 l2 g = 9.81; q1 = y(1); q2 = y(2); qd1 = y(3); qd2 = y(4); qd = [qd1; qd2]; 211 M_q = [m1*l1^2+m2*l1^2, m2*l1*l2*cos(q1-q2); m2*l1*l2*cos(q1-q2), m2*l2^2 ]; C_q = [ 0, m2*l1*l2*sin(q1-q2)*qd2; -m2*l1*l2*sin(q1-q2)*qd1, 0]; g_q = [g*l1*sin(q1)*(m1+m2); g*m2*l2*sin(q2)]; ydot = zeros(4,1); qdot = qd; qd_dot = inv(M_q)*(-C_q*qd - g_q); ydot = [qdot; qd_dot]; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % main program global m1 m2 l1 l2 m1 = 0.25; l1 = 0.5; m2 = 0.25; l2 = 0.5; t_start = 0; t_end = 6.0; q10 = 10; q20 = 12; y1 = q10*pi/180; y2 = q20*pi/180; y3 = 0; y4 = 0; y0 = [y1 y2 y3 y4]'; [t,y] = ode45('conlactoan_kep',[0: 0.01: t_end],y0); plot(t, y(:,1),'k-', t, y(:,2),'k--', 'linewidth',1), grid on legend('q_1','q_2'), xlabel('t [s]') axis([0 t_end -0.3 0.3]); Chy chương trình vi các điu kin đu khác nhau cho ta các kt qu như trên hình 926 và 927. 0.3 0 0.2 -0.2 0.1 -0.4 0 -0.1 y y [m] -0.6 -0.2 -0.8 0 2 4 6 t [s] -1 q q -0.4 -0.2 0 0.2 1 2 x [m] Hình 926. Kt qu mơ phng vi điu kin đu q10 = 10; q20 = 12; 212 0.3 0 0.2 -0.2 0.1 -0.4 0 y [m] y -0.6 -0.1 -0.2 -0.8 0 1 2 3 -1 t [s] q q -0.4 -0.2 0 0.2 1 2 x [m] Hình 927. Kt qu mơ phng vi điu kin đu q10 = 10; q20 = 12; Dao đng nh ca con lc elliptic Xét con lc elliptic như hình 928. Xe A cĩ khi x lưng m1 chuyn đng trên đưng ngang, lị xo cĩ đ m1 cng c . Dây AB cĩ chiu dài l, khi lưng khơng c đáng k và luơn căng. Ti trng đưc coi như cht A đim cĩ khi lưng m2. Bit khi x = 0 lị xo khơng b bin dng. l S dng phương trình Lagrange loi 2, ta nhn đưc ϕ m phương trình vi phân chuyn đng: 2 B 2 (mmxml1++ 2 )ɺɺ 2ϕɺɺ cos ϕ − ml 2 ϕ ɺ sin ϕ += cx 0, Hình 928 2 ml2ϕɺɺ + mlx 2ɺɺcos ϕ + mgl 2 sin ϕ = 0 Trưng hp xét xét dao đng nh, coi sinϕ≈ ϕ , cos ϕ ≈ 1 và b qua s hng phi tuyn ϕɺ2 , ta nhn đưc phương trình vi phân dao đng nh. 2 (m1+++= mxml 2 )ɺɺ 2ϕɺɺ cx 0, mlxml 2 ɺɺ ++ 2 ϕ ɺɺ mgl 2 ϕ = 0 hay dng ma trn m+ m ml ɺɺ 1 2 2  x c 0  x  0    +  =   . ml ml 2 ϕɺɺ 0 m gl ϕ  0 2 2  2      T đây ta s s dng Matlab tìm các tn s dao đng riêng và các dng dao đng riêng. 213 function Daodongnho_enliptic m1 = 50; m2 = 40; % kg l = 5; % m g = 9.81; % m/s^2 c = 2000; % N/m M = [m1 + m2, m2*l; m2*l m2*l^2]; C = [c, 0; 0, m2*g*l]; n = size(M,1); [V, D] = eig(C,M); for i=1:n ome(i)=sqrt(D(i,i)); end ome, v1 = V(:,1), v2 = V(:,2) x = [0:1:n]; subplot(2,1,1) plot(x, [0; V(:,1)],'k-','linewidth',1), grid on, ylabel('v_1') subplot(2,1,2) plot(x,[0; V(:,2)],'k-','linewidth',1), grid on, ylabel('v_2') Các tn s riêng và các dng 0 dao đng riêng 1 -0.02 ome = v 1.3727 -0.04 6.4535 0 0.5 1 1.5 2 0.2 v1 = -0.0062 2 0 -0.0303 v v2 = -0.2 0 0.5 1 1.5 2 -0.1413 0.0297 Hình 929. Các dng dao đng riêng H dao đng ba bc t do Cho h dao đng ba bc t do như trên hình 930. Các lị xo tuyn tính cĩ đ cng ci , khi lưng các vt nng là mi , (i = 1,2, 3 ). Gi xi là dch chuyn ca các khi lưng k t v trí mà các lị xo khơng b bin dng. Áp dng phương pháp tách vt (hoc phương trình Lagrange loi 2) nhn đưc phương trình vi phân chuyn đng dng ma trn như sau: T Mxɺɺ + Cx = p , x = [x1 x 2 x 3 ] 214 vi   m1 0 0 m g  c c   1  1 2 M = 0m 0  , p = m g 2 2  m x1     g 1 0 0 m m3 g 3    c3 c+++ c c c − c − c  x c5 1 2 3 5 3 5  2 m2 C =−c3 c 3 +− c 4 c 4  , c4   −c5 − c 4 c 4 + c 5   m3 x3 Cho bit các thơng s ca h: Hình 930 c1=== c 2 c 310 N/mm, c 4 == c 5 5 N/mm , m1 g== m 3 g98.1 N, m 2 g = m 3 g . 1. Tìm v trí cân bng tĩnh ca h, bng cách gii h phương trình đi s tuyn tính đ tìm bin dng ca các lị xo Cx= p . 2. Tìm các tn s dao đng riêng và các dng dao đng riêng ca h t phương trình dao đng t do Mxɺɺ + Cx = 0 . Bài tốn đưc gii trong Matlab bng mt mfile như sau: function Daodon_3DOF m1 = 10; m2 = 5; m3 = m1; % kg c1 = 10000; % N/m c2 = c1; c3 = c1; c4 = 5000; c5 = c4; g = 9.81; % m/s^2 M = diag([m1, m2, m3]); C = [ c1+c2+c3+c5, -c3, -c5; -c3, c3+c4, -c4; -c5, -c4, c4+c5]; p = g*[m1; m2; m3]; x0 = inv(C)*p % tinh do dan tinh cua cac lo xo n = size(M,1); [V, D] = eig(C,M); for i=1:n ome(i)=sqrt(D(i,i)); end ome' v1 = V(:,1) v2 = V(:,2) v3 = V(:,3) 215 x = [0:1:n]; subplot(3,1,1) plot(x,[0; V(:,1)],'k-','linewidth',1), grid on, ylabel('v_1') subplot(3,1,2) plot(x,[0; V(:,2)],'k-','linewidth',1), grid on, ylabel('v_2') subplot(3,1,3) plot(x,[0; V(:,3)],'k-','linewidth',1), grid on, ylabel('v_3') ð dãn tĩnh ca lị xo x0 = 0.0123 0.0201 0.0260 Các tn s dao đng riêng ome = 0.4 21.2572 1 0.2 48.5851 v 68.4662 0 Các véctơ dng riêng 0 1 2 3 v1 = 0.5 0.1053 2 0 0.1875 v 0.2671 -0.5 v2 = 0 1 2 3 0.1756 0.5 0.2847 3 0 -0.1692 v v3 = -0.5 -0.2410 0 1 2 3 0.2895 -0.0066 Hình 931. Các dng dao đng riêng 9.6 Phân tích đng hc cơ cu Bài tốn phân tích đng hc cơ cu da trên các phương trình liên kt thưng dn đn vic gii h phương trình đi s phi tuyn. Sau đây trình bày vic gii h phương trình đi s phi tuyn bng phương pháp NewtonRaphson. Bài tốn: cn tìm nghim x tha mãn h phương trình đi s phi tuyn fx0(),= f ∈RRn , x ∈ n hay vit dưi dng h n phương trình fxxxk(1 , 2 ,..., n )= 0, k = 1,2,..., n . 216 ð thit lp các bưc ca bài tính, trưc ht ta khai trin Taylor hàm f( x ) ti lân cn x0 ∂f fx(+=+δ xfx )() () xx δ +O () δ x 2 0 0∂x 0 2 =fx()()0 + Jxx 0 δ +O () δ x vi J( x0 ) là ma trn Jacobi đưc xác đnh ti đim x0 ∂f1 ∂ f 1 ∂ f 1  ..  ∂xx1 ∂ 2 ∂ xn ∂f ∂ f ∂ f  ∂f 2 2.. 2  Jxx()()0= 0 = ∂xx1 ∂ 2 ∂ xn  () x0 ∂x .. .. .. ..  ∂f ∂ f ∂ f  nn.. n  ∂xx1 ∂ 2 ∂ xn  Gi s x0 là mt xp x ban đu ca nghim cn tìm và đ (x0 + δ x ) là mt giá tr gn nghim hơn (chính xác hơn), ta cho f( x0 +δ x ) = 0 . T đĩ nhn đưc mt h phương trình đi s tuyn tính đi vi δx −1 Jx()0δ x= − fx () 0 suy ra δx= − Jx()0 ⋅ fx () 0 . Ly giá tr x1= x 0 + δ x là xp x gn đúng mi ca nghim. Lp li nhiu ln cơng vic trên, nu phương pháp hi t s cho ta nghim ca h phương trình đi s phi tuyn. Thut gii Bưc 1. Khi gán k = 0 , chn xp x ban đu x(0) . Bưc 2. Tính f( x(k ) ) . Nu ||f ( x(k ) )||< ε thì dng, nu khơng thì tip tc t 3. Bưc 3. Tính ma trn Jacobi ti x(k ) , tc là J( x(k ) ). Bưc 4. Gii Jx()()kδ x= − fx () () k đ tìm δx(k ) . Bưc 5. Ly x(1)k+ = x () k + δ x . Bưc 6. Tăng k , k= k + 1 . Nu k > M vi M là s bưc lp đã chn trưc, thì dng. Nu khơng, tip tc t Bưc 2. Trong thut gii trên ta đã thốt ra khi vịng lp trong hai tình hung. Khi khơng đt đưc đ chính xác cn thit, ta buc phi kt thúc sau mt s bưc lp nht đnh (bưc 6). Quá trình lp s kt thúc tt nu như ta đt đưc đ chính xác chn trưc (bưc 2). Chú ý rng đây ta chn điu kin dng bng cách xét chun Euclide ca vectơ f. Cũng cĩ th thay điu kin này bng các điu kin khác như 217 ||δx ||< ε . S thành cơng ca phương pháp ph thuc vào vic chn xp x ban đu. Nu chn đưc đim xut phát tt, phương pháp hi t rt nhanh. Ngưc li, thì khĩ cĩ th đốn đưc và yêu cu thut tốn phi dng li sau mt s bưc lp nht đnh. Khi áp dng phương pháp NewtonRaphson đ gii bài tốn, ngồi vic nhp các hàm f( x ) vào, cịn yêu cu phi cĩ ma trn Jacobi J( x ) . Trong trưng hp nu vic đo hàm ∂fi/ ∂ x j gp khĩ khăn, hoc khĩ thc hin đưc ta cĩ th s dng xp x sau đây đ tính ma trn jacobi ∂fi fi(xe+ j h ) − f i () x ≈ ∂xj h vi h là bưc nh và ej là véctơ đơn v theo hưng ca trc x j . Phân tích đng hc cơ cu bn khâu bn l. Cho cơ cu 4 khâu bn l OABC như hình 9 y 32. Các kích thưc ca cơ cu OA = l2, AB = B l3, BC = l4, dx và dy. Hãy vit các phương trình liên kt cho cơ cu, t đĩ gii bài tốn v trí ca cơ cu bng phương pháp Newton q A 2 q3 Raphson. Cho bit qq1= 10 + t , hãy xác đnh các gĩc ϕ3(t ) và ϕ4(t ) ti các thi đim q1 dy C tk , k = 1,2,3,... ng vi các gĩc quay ca O x thanh OA, qq1= 10 + tk . Cho bit các dx thơng s ca cơ cu l1=0.20, l 2 = 0.60, Hình 932 l3 =0.40, dx = 0.45, d y = 0.12 m. T hình v ta vit đưc các phương trình liên kt fqqq1(,,) 1 2 3= l 1 cos+cos ql 1 2 qdl 2 −−=x 3 cos q 3 0 fqqq2(,,) 1 2 3= l 1 sin ql 1 + 2 sin qdl 2 −−y 3 sin q 3 = 0 Tính đưc các ma trn jacobi −l2sin q 2 l 3 sin q 3  −l1sin q 1  Jq =   , Jq =   23 lcos q− l cos q 1 lcos q 2 2 3 3  1 1  ðo hàm các phương trình liên kt theo thi gian cho ta phương trình đ gii bài tốn vn tc và gia tc: J qɺ+ J q ɺ = 0 ⇒qɺ =− J−1 J q ɺ q2323 q 1 1 23q23 q 1 1 J qɺɺ+ Jɺ q ɺ + Jq ɺɺ +J ɺ q ɺ =⇒0 qɺɺ=− J−1 ( Jɺ q ɺ + Jq ɺɺ + Jɺ q ɺ ) q2323 q 23 23 q 1 1 q 1 1 23q23 q 23 23 q 1 1 q 1 1 218 vi −qɺ lcos q q ɺ l cos q −lcos q ɺ 2 2 2 3 3 3  ɺ 1 1  Jq =   , Jq =  qɺ1 23 −qɺ lsin q q ɺ l sin q 1 −lsin q 2 2 2 3 3 3  1 1  Trin khai trong Matlab vi các mfile 1. Các phương trình liên kt function f = bkbl_ ptlk(q1, q2, q3) global l1 l2 l3 dx dy f1 = l1*cos(q1) + l2*cos(q2) - l3*cos(q3) - dx; f2 = l1*sin(q1) + l2*sin(q2) - l3*sin(q3) - dy; f=[f1; f2]; 2. Các ma trn Jacobi function [J1, J23] = bkbl_Jacob(q1,q2,q3) global l1 l2 l3 dx dy J1 = [-l1*sin(q1); l1*cos(q1)]; J23 = [-l2*sin(q2), l3*sin(q3); l2*cos(q2), -l3*cos(q3)]; 3. ðo hàm ca ma trn Jacobi function [dJ1, dJ23] = bkbl_Jacdot(q1,q2,q3, dq1, dq2,dq3) global l1 l2 l3 dx dy dJ1 = -l1*dq1*[cos(q1); sin(q1)]; dJ23 = [-dq2*l2*cos(q2), dq3*l3*cos(q3); -dq3*l2*sin(q2), dq3*l3*sin(q3)]; 4. Chương trình chính function bkbl_main global l1 l2 l3 dx dy l1 = 0.12; l2 = 0.70; l3 = 0.40; dx = 0.50; dy = 0.10; %cho goc cua tay quay OA q10 = 0.5; % rad omega = 1.0; % rad/s t = linspace(0, 5*pi/omega, 360); a(1) = 0.2; b(1) = 0.542; % chon xap xi ban dau cho q2, q3 for i=1:length(t) q1(i)= q10+omega*t(i); dq1(i) = omega; ddq1(i) = 0; f = ptlk_bkbl(q1(i), a(1), b(1)); k = 1; while (norm(f,2) > 1.0e-10 && k < 30) [J1, J23] = Jacob_bkbl(q1(i), a(1), b(1)); deltax = inv(J23)*(-f); a(2) = a(1)+deltax(1); 219 b(2) = b(1)+deltax(2); a(1) = a(2); b(1) = b(2); f = ptlk_bkbl(q1(i), a(1), b(1)); k = k+1; end q2(i)=a(1); q3(i)=b(1); % van toc [J1, J23] = Jacob_bkbl(q1(i), q2(i), q3(i)); dq23 = -inv(J23)*J1*dq1(i); dq2(i) = dq23(1); dq3(i) = dq23(2); % gia toc [dJ1, dJ23] = Jacdot_bkbl(q1(i),q2(i),q3(i),dq1(i),dq2(i),dq3(i)); ddq23 = -inv(J23)*(dJ23*dq23 + J1*ddq1(i) + dJ1*dq1(i)); ddq2(i) = ddq23(1); ddq3(i) = ddq23(2); end figure(1) plot(t,q2,'k-', t, dq2,'k--', 'linewidth',1), grid on xlabel('t [s]'), legend('q_2', 'dq_2') axis([0 max(t) -0.5 1.1]) figure(2) plot(t,q3,'k-', t, dq3,'k--', 'linewidth',1), grid on xlabel('t [s]'), legend('q_3', 'dq_3') axis([0 max(t) -0.8 2.2]) x0 = 0; y0 = 0; r = 0.005; figure(3), hold on plot(x0,y0,'o-','linewidth',1), plot(dx,dy,'o-','linewidth',1) for i = 1:5:length(t)/2 x1(i) = l1*cos(q1(i)); x2(i) = x1(i) + l2*cos(q2(i)); y1(i) = l1*sin(q1(i)); y2(i) = y1(i) + l2*sin(q2(i)); x4(i) = dx; x3(i) = dx + l3*cos(q3(i)); y4(i) = dy; y3(i) = dy + l3*sin(q3(i)); x_day = [x0, x1(i), x2(i), x3(i), x4(i), x0]; y_day = [y0, y1(i), y2(i), y3(i), y4(i), y0]; plot(x_day, y_day, 'k-','linewidth',1), grid on, axis equal tam1 = [x1(i), y1(i)]; vehinhtron(tam1,r); tam2 = [x2(i), y2(i)]; vehinhtron(tam2,r); end axis([-0.2 0.8 -0.2 0.6]) ð th ca các gĩc q2, q 3 cùng các đo hàm qɺ2, q ɺ 3 đưc th hin như trên hình 9 33. Cu hình ca cơ cu ti mt s v trí đưc ch ra trên hình 934. 220 1 2 0.5 1 0 0 -0.5 0 5 10 15 0 5 10 15 t [s] t [s] q dq q2 dq2 3 3 Hình 933. ð th các gĩc q2 và q3 theo thi gian 0.6 0.5 0.4 0.3 0.2 0.1 0 -0.1 -0.2 -0.2 0 0.2 0.4 0.6 0.8 Hình 934. Cu hình cơ cu bn khâu bn l 9.7 Bài tốn đng hc ngưc rơbt Trong bài tốn đng hc ngưc rơbt, ta cn xác đnh các gĩc khp – hay ta đ khp – khi mun bàn kp ca rơbt mt v trí mong mun. Trong phn này ta xét mt tay máy hai bc t do phng cĩ sơ đ như trên hình 935. Cn tìm các ta đ suy rng q1, q 2 đ bàn kp E mt v trí mong mun cho trưc xE, y E . Bit chiu dài các khâu là OA= L1 = 0.50 m , AE= L2 = 0.70 m . ði vi bài tốn này, rơbt cĩ cu trúc khá đơn gin, ta cĩ th gii bng phương pháp hình hc, bng cách gii tam giác OAE khi bit chiu dài ba cnh ca nĩ. 221 đây ta s s dng phương pháp gii h phương trình phi tuyn, đ gii bài tốn này, E q vì l phương pháp cĩ th đưc m rng d 2 yE dàng cho các rơbt cĩ s bc t do ln hơn và A cĩ cu trúc phc tp hơn. Mi liên h gia ta đ đim E vi các ta đ q1 suy rng đưc xác đnh bi O xE xE = l1cos q 1 + l 2 cos( q 1 − q 2 ) yE = l1sin q 1 + l 2 sin( q 1 − q 2 ) Hình 935 ðây là mt h hai phương trình phi tuyn đi vi hai n q1, q 2 . Phương pháp NewtonRaphson đưc s dng đ gii h này, vi ma trn Jacobi như sau −l1sin q 1 l 2 sin( q 1 − q 2 )  J =   q lcos q− l cos( q − q ) 1 1 2 1 2  Trin khai trong Matlab vi các mfile sau: 1. Phương trình liên kt function f = taymay2dof_bkbl(q1, q2, xE, yE) global l1 l2 f1 = l1*cos(q1) + l2*cos(q1 - q2) - xE; f2 = l1*sin(q1) + l2*sin(q1 - q2) - yE; f=[f1; f2]; 2. Ma trn Jacobi function Jq = taymay2dof_Jacob(q1, q2) global l1 l2 Jq = [-l1*sin(q1), l2*sin(q1 - q2); l1*cos(q1), -l2*cos(q1 - q2)]; 3. Chương trình chính function taymay2dof_main global l1 l2 l1 = 0.50; l2 = 0.70; n = 20; % so diem chia xE = linspace(1.1, 0.1, n); yE = linspace(0.1, 1.1, n); a(1) = 0.5; b(1) = 0.2; % chon xap xi ban dau cho q1, q2 for i=1:n f = taymay2dof_ptlk(a(1), b(1), xE(i), yE(i)); k = 1; while (norm(f,2) > 1.0e-10 && k < 30) 222 Jq = taymay2dof_Jacob(a(1), b(1)); dq = inv(Jq)*(-f); a(2) = a(1)+dq(1); b(2) = b(1)+dq(2); a(1) = a(2); b(1) = b(2); f = taymay2dof_ptlk(a(1), b(1), xE(i), yE(i)); k = k+1; end q1(i) = a(1); q2(i) = b(1); end k = 1:1:n; figure(1) plot(k, q1,'k-', k, q2,'k--', 'linewidth',1), grid on legend('q_1', 'q_2'), xlabel('k = 1..n'), ylabel('[rad]') axis([0 n+1 0 2.1]) figure(2) x0 = 0; y0 = 0; plot(x0,y0,'o-','linewidth',1), hold on for i = 1:n xA(i) = l1*cos(q1(i)); xE(i) = xA(i)+l2*cos(q1(i)-q2(i)); yA(i) = l1*sin(q1(i)); yE(i) = yA(i)+l2*sin(q1(i)-q2(i)); x_day = [x0, xA(i), xE(i)]; y_day = [y0, yA(i), yE(i)]; plot(x_day, y_day, 'k-','linewidth',1), grid on, axis equal end Chy chương trình cho ta các kt qu như trong bng 95. ð th ca các gĩc khp và cu hình tương ng ca rơbt đưc đưa ra như trên hình 936. 2 1 0.8 1.5 0.6 1 y [m] 0.4 [rad] 0.2 q 0.5 1 0 q2 0 -0.2 0 5 10 15 20 0 0.5 1 k = 1..n x [m] Hình 936. ð th các ta đ q1, q 2 và cu hình ca rơbt 223 Bng 95. Kt qu tính tốn bài tốn đng hc ngưc Cho trưc Gii tìm đưc k xE [m] yE [m] q1 [rad] q2 [rad] 1 1.1000 0.1000 0.5701 0.8152 2 1.0474 0.1526 0.7334 0.9965 3 0.9947 0.2053 0.8806 1.1411 4 0.9421 0.2579 1.0181 1.2600 5 0.8895 0.3105 1.1488 1.3584 6 0.8368 0.3632 1.2738 1.4387 7 0.7842 0.4158 1.3931 1.5024 8 0.7316 0.4684 1.5063 1.5499 9 0.6789 0.5211 1.6123 1.5816 10 0.6263 0.5737 1.7098 1.5974 11 0.5737 0.6263 1.7974 1.5974 12 0.5211 0.6789 1.8739 1.5816 13 0.4684 0.7316 1.9380 1.5499 14 0.4158 0.7842 1.9889 1.5024 15 0.3632 0.8368 2.0257 1.4387 16 0.3105 0.8895 2.0478 1.3584 17 0.2579 0.9421 2.0545 1.2600 18 0.2053 0.9947 2.0444 1.1411 19 0.1526 1.0474 2.0148 0.9965 20 0.1000 1.1000 1.9596 0.8152 224 Tài liu tham kho [1] Ottmar Beucher: Matlab und Simulink: Grundlegende Einführung für Studenten und Ingenierure in der Praxis. Pearson Studium Verlag; Auflage 3 (2006). [2] Helmut Bode: MATLABSimulink Analyse und Simulation dynamischer Systeme. Springer Verlag, 2006. [3] Jưrg Kahlert: Simulation technischer Systeme: Eine beispielorientierte Einführung. Vieweg & SohnVerlag, Wiesbaden, 2004. [4] Jaan Kiusalaas: Numerical Methods in Engineering with MATLAB. Cambridge University Press, 2005. [5] David McMahon: MATLAB Demystified A SelfTeaching Guide. McGraw Hill, 2007. [6] John H. Mathews, Kurtis D. Fink: Numerical methods using MATLAB. 3. Ed. Prentice Hall 1999. [7] Dieter Pietruszka: MATLAB in der Ingenieurpraxis. Modellbildung, Berechnung und Simulation. 2. Ed. Vieweg+Teubner Verlag, 2006. [8] L.F. Shampine, I. Gladwell, S. Thomson: Solving ODEs with Matlab. Cambridge University Press, 2003. [9] Won Young Yang, Wenwu Cao, TaeSang Chung, John Morris: Applied Numerical Methods Using Matlab. John Wiley & Sons, Inc., 2005. [10] D. Gross, W. Hauger, W. Schnell, J. Schrưder: Technische Mechanik, Band 1: Statik. (8. Auflage). Springer Verlag, Berlin 2004. [11] D. Gross, W. Hauger, W. Schnell, J. Schrưder: Technische Mechanik, Band 2: Elastostatik (8. Auflage). Springer Verlag, Berlin 2005. [12] D. Gross, W. Hauger, W. Schnell, J. Schrưder: Technische Mechanik, Band 3: Kinetik (8. Auflage). Springer Verlag, Berlin 2004. [13] Dankert/Dankert: Technische Mechanik: Statik, Festigkeitslehre, Kinematik/Kinetik. (3. Auflage). Teubner B.G. 2004. [14] R.C. Hibbele: Engineering Mechanics –Strength of Materials (10.Edition). Prentice Hall, Upper Saddle River , New Jersey 2004. [15] Nguyn Hồng Hi. Nguyn Vit Anh: Lp trình Matlab và ng dng. NXB Khoa hc và K thut, Hà Ni 2006. [16] Nguyn Văn Khang: Dao đng k thut (in ln 4). NXB Khoa hc và K thut, Hà Ni 2005. [17] ðinh Văn Phong: Phương pháp s trong cơ hc. Nhà xut bn Khoa hc và K thut, 1999. 225 [18] Nguyn Phùng Quang: Matlab & Simulink dành cho k sư điu khin t đng. NXB Khoa hc và K thut, Hà Ni 2006. [19] ð Sanh: Cơ hc k thut, tp 1 và tp 2. Nhà xut bn Giáo dc, Hà ni 2008. [20] Phm Th Ngc Yn, Ngơ Hu Tình, Lê Tn Hùng và Nguyn Th Lan Hương: Cơ s Matlab & ng dng. NXB Khoa hc và K thut, Hà Ni, 2007. [21] Lê Quang Minh, Nguyn Văn Vưng: Sc bn vt liu, Tp 1,2,3. NXB Giáo dc, Hà ni 1997. [22] ðng Vit Cương, Nguyn Nht Thăng, Nh Phương Mai: Sc bn vt liu, tp 1 và 2. NXB Khoa Hc và K Thut, 2002. [23] Phm Minh Hồng: Maple và các bài tốn ng dng. NXB Khoa hc và K thut, 2008. 226